homework 1 glm (lesson 6).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0524a136-e54d-4cd7-80e0-43d3a75a76fa",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:14:42.959839Z",
     "shell.execute_reply.started": "2024-08-06T09:14:42.955263Z",
     "to_execute": "2024-08-06T09:14:42.952Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [],
   "source": [
    "import ssl\n",
    "ssl._create_default_https_context=ssl._create_unverified_context"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7d6e9d16-f767-431f-a7b9-a882c5f78c65",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:14:43.489726Z",
     "shell.execute_reply.started": "2024-08-06T09:14:42.962338Z",
     "to_execute": "2024-08-06T09:14:42.952Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from functools import partial\n",
    "from sklearn.datasets import fetch_openml\n",
    "from sklearn.metrics import (\n",
    "    mean_absolute_error,\n",
    "    mean_squared_error,\n",
    "    mean_tweedie_deviance,\n",
    ")\n",
    "\n",
    "\n",
    "def load_mtpl2(n_samples=None):\n",
    "    df_freq = fetch_openml(data_id=41214, as_frame=True).data\n",
    "    df_freq[\"IDpol\"] = df_freq[\"IDpol\"].astype(int)\n",
    "    df_freq.set_index(\"IDpol\", inplace=True)\n",
    "\n",
    "    df_sev = fetch_openml(data_id=41215, as_frame=True).data\n",
    "\n",
    "    df_sev = df_sev.groupby(\"IDpol\").sum()\n",
    "\n",
    "    df = df_freq.join(df_sev, how=\"left\")\n",
    "    df[\"ClaimAmount\"].fillna(0, inplace=True)\n",
    "\n",
    "    for column_name in df.columns[df.dtypes.values == object]:\n",
    "        df[column_name] = df[column_name].str.strip(\"'\")\n",
    "    return df.iloc[:n_samples]\n",
    "\n",
    "def plot_obs_pred(\n",
    "    df,\n",
    "    feature,\n",
    "    weight,\n",
    "    observed,\n",
    "    predicted,\n",
    "    y_label=None,\n",
    "    title=None,\n",
    "    ax=None,\n",
    "    fill_legend=False,\n",
    "):\n",
    "    df_ = df.loc[:, [feature, weight]].copy()\n",
    "    df_[\"observed\"] = df[observed] * df[weight]\n",
    "    df_[\"predicted\"] = predicted * df[weight]\n",
    "    df_ = (\n",
    "        df_.groupby([feature])[[weight, \"observed\", \"predicted\"]]\n",
    "        .sum()\n",
    "        .assign(observed=lambda x: x[\"observed\"] / x[weight])\n",
    "        .assign(predicted=lambda x: x[\"predicted\"] / x[weight])\n",
    "    )\n",
    "\n",
    "    ax = df_.loc[:, [\"observed\", \"predicted\"]].plot(style=\".\", ax=ax)\n",
    "    y_max = df_.loc[:, [\"observed\", \"predicted\"]].values.max() * 0.8\n",
    "    p2 = ax.fill_between(\n",
    "        df_.index,\n",
    "        0,\n",
    "        y_max * df_[weight] / df_[weight].values.max(),\n",
    "        color=\"g\",\n",
    "        alpha=0.1,\n",
    "    )\n",
    "    if fill_legend:\n",
    "        ax.legend([p2], [\"{} distribution\".format(feature)])\n",
    "    ax.set(\n",
    "        ylabel=y_label if y_label is not None else None,\n",
    "        title=title if title is not None else \"Train: Observed vs Predicted\",\n",
    "    )\n",
    "\n",
    "\n",
    "def score_estimator(\n",
    "    estimator,\n",
    "    X_train,\n",
    "    X_test,\n",
    "    df_train,\n",
    "    df_test,\n",
    "    target,\n",
    "    weights,\n",
    "    tweedie_powers=None,\n",
    "):\n",
    "    \"\"\"Evaluate an estimator on train and test sets with different metrics\"\"\"\n",
    "\n",
    "    metrics = [\n",
    "        (\"D² explained\", None),  # Use default scorer if it exists\n",
    "        (\"mean abs. error\", mean_absolute_error),\n",
    "        (\"mean squared error\", mean_squared_error),\n",
    "    ]\n",
    "    if tweedie_powers:\n",
    "        metrics += [\n",
    "            (\n",
    "                \"mean Tweedie dev p={:.4f}\".format(power),\n",
    "                partial(mean_tweedie_deviance, power=power),\n",
    "            )\n",
    "            for power in tweedie_powers\n",
    "        ]\n",
    "\n",
    "    res = []\n",
    "    for subset_label, X, df in [\n",
    "        (\"train\", X_train, df_train),\n",
    "        (\"test\", X_test, df_test),\n",
    "    ]:\n",
    "        y, _weights = df[target], df[weights]\n",
    "        for score_label, metric in metrics:\n",
    "            if isinstance(estimator, tuple) and len(estimator) == 2:\n",
    "                # Score the model consisting of the product of frequency and\n",
    "                # severity models.\n",
    "                est_freq, est_sev = estimator\n",
    "                y_pred = est_freq.predict(X) * est_sev.predict(X)\n",
    "            else:\n",
    "                y_pred = estimator.predict(X)\n",
    "\n",
    "            if metric is None:\n",
    "                if not hasattr(estimator, \"score\"):\n",
    "                    continue\n",
    "                score = estimator.score(X, y, sample_weight=_weights)\n",
    "            else:\n",
    "                score = metric(y, y_pred, sample_weight=_weights)\n",
    "\n",
    "            res.append({\"subset\": subset_label, \"metric\": score_label, \"score\": score})\n",
    "\n",
    "    res = (\n",
    "        pd.DataFrame(res)\n",
    "        .set_index([\"metric\", \"subset\"])\n",
    "        .score.unstack(-1)\n",
    "        .round(4)\n",
    "        .loc[:, [\"train\", \"test\"]]\n",
    "    )\n",
    "    return res\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b26f5a69-218a-447a-adda-2751690ea8c7",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:18.080469Z",
     "shell.execute_reply.started": "2024-08-06T09:14:43.491153Z",
     "to_execute": "2024-08-06T09:14:42.952Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.10/site-packages/sklearn/datasets/_openml.py:1022: FutureWarning: The default value of `parser` will change from `'liac-arff'` to `'auto'` in 1.4. You can set `parser='auto'` to silence this warning. Therefore, an `ImportError` will be raised from 1.4 if the dataset is dense and pandas is not installed. Note that the pandas parser may return different data types. See the Notes Section in fetch_openml's API doc for details.\n",
      "  warn(\n",
      "/usr/local/lib/python3.10/site-packages/sklearn/datasets/_openml.py:1022: FutureWarning: The default value of `parser` will change from `'liac-arff'` to `'auto'` in 1.4. You can set `parser='auto'` to silence this warning. Therefore, an `ImportError` will be raised from 1.4 if the dataset is dense and pandas is not installed. Note that the pandas parser may return different data types. See the Notes Section in fetch_openml's API doc for details.\n",
      "  warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       ClaimNb  Exposure Area  VehPower  VehAge  DrivAge  BonusMalus VehBrand  \\\n",
      "IDpol                                                                           \n",
      "139.0      1.0      0.75    F       7.0     1.0     61.0        50.0      B12   \n",
      "190.0      1.0      0.14    B      12.0     5.0     50.0        60.0      B12   \n",
      "414.0      1.0      0.14    E       4.0     0.0     36.0        85.0      B12   \n",
      "424.0      2.0      0.62    F      10.0     0.0     51.0       100.0      B12   \n",
      "463.0      1.0      0.31    A       5.0     0.0     45.0        50.0      B12   \n",
      "\n",
      "        VehGas  Density Region  ClaimAmount   PurePremium  Frequency  \\\n",
      "IDpol                                                                  \n",
      "139.0  Regular  27000.0    R11       303.00    404.000000   1.333333   \n",
      "190.0   Diesel     56.0    R25      1981.84  14156.000000   7.142857   \n",
      "414.0  Regular   4792.0    R11      1456.55  10403.928571   7.142857   \n",
      "424.0  Regular  27000.0    R11     10834.00  17474.193548   3.225806   \n",
      "463.0  Regular     12.0    R73      3986.67  12860.225806   3.225806   \n",
      "\n",
      "       AvgClaimAmount  \n",
      "IDpol                  \n",
      "139.0          303.00  \n",
      "190.0         1981.84  \n",
      "414.0         1456.55  \n",
      "424.0         5417.00  \n",
      "463.0         3986.67  \n"
     ]
    }
   ],
   "source": [
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import (\n",
    "    FunctionTransformer,\n",
    "    KBinsDiscretizer,\n",
    "    OneHotEncoder,\n",
    "    StandardScaler,\n",
    ")\n",
    "\n",
    "df = load_mtpl2()\n",
    "\n",
    "df.loc[(df[\"ClaimAmount\"] == 0) & (df[\"ClaimNb\"] >= 1), \"ClaimNb\"] = 0\n",
    "\n",
    "df[\"ClaimNb\"] = df[\"ClaimNb\"].clip(upper=4)\n",
    "df[\"Exposure\"] = df[\"Exposure\"].clip(upper=1)\n",
    "df[\"ClaimAmount\"] = df[\"ClaimAmount\"].clip(upper=200000)\n",
    "\n",
    "log_scale_transformer = make_pipeline(\n",
    "    FunctionTransformer(func=np.log), StandardScaler()\n",
    ")\n",
    "\n",
    "column_trans = ColumnTransformer(\n",
    "    [\n",
    "        (\n",
    "            \"binned_numeric\",\n",
    "            KBinsDiscretizer(n_bins=10, subsample=int(2e5), random_state=0),\n",
    "            [\"VehAge\", \"DrivAge\"],\n",
    "        ),\n",
    "        (\n",
    "            \"onehot_categorical\",\n",
    "            OneHotEncoder(),\n",
    "            [\"VehBrand\", \"VehPower\", \"VehGas\", \"Region\", \"Area\"],\n",
    "        ),\n",
    "        (\"passthrough_numeric\", \"passthrough\", [\"BonusMalus\"]),\n",
    "        (\"log_scaled_numeric\", log_scale_transformer, [\"Density\"]),\n",
    "    ],\n",
    "    remainder=\"drop\",\n",
    ")\n",
    "X = column_trans.fit_transform(df)\n",
    "\n",
    "df[\"PurePremium\"] = df[\"ClaimAmount\"] / df[\"Exposure\"]\n",
    "\n",
    "df[\"Frequency\"] = df[\"ClaimNb\"] / df[\"Exposure\"]\n",
    "df[\"AvgClaimAmount\"] = df[\"ClaimAmount\"] / np.fmax(df[\"ClaimNb\"], 1)\n",
    "\n",
    "with pd.option_context(\"display.max_columns\", 15):\n",
    "    print(df[df.ClaimAmount > 0].head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e90e600f-b598-4db3-860d-18297d3be98b",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:18.208843Z",
     "shell.execute_reply.started": "2024-08-06T09:15:18.082285Z",
     "to_execute": "2024-08-06T09:14:42.952Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "df_train, df_test, X_train, X_test = train_test_split(df, X, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7295219f-9a00-4aad-8e4c-9fc777caf58e",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:20.045514Z",
     "shell.execute_reply.started": "2024-08-06T09:15:18.210905Z",
     "to_execute": "2024-08-06T09:14:42.952Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluation of the Product Model and the Tweedie Regressor on target PurePremium\n",
      "                                      T              \n",
      "subset                            train          test\n",
      "metric                                               \n",
      "D² explained               1.690000e-02  1.420000e-02\n",
      "mean Tweedie dev p=1.5000  7.640770e+01  7.640880e+01\n",
      "mean Tweedie dev p=1.7000  3.682880e+01  3.692270e+01\n",
      "mean Tweedie dev p=1.8000  3.037600e+01  3.045390e+01\n",
      "mean Tweedie dev p=1.9000  3.382120e+01  3.387830e+01\n",
      "mean Tweedie dev p=1.9900  2.015347e+02  2.015587e+02\n",
      "mean Tweedie dev p=1.9990  1.914538e+03  1.914387e+03\n",
      "mean Tweedie dev p=1.9999  1.904747e+04  1.904558e+04\n",
      "mean abs. error            2.739865e+02  2.731249e+02\n",
      "mean squared error         3.295505e+07  3.213056e+07\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import TweedieRegressor\n",
    "\n",
    "glm_pure_premium = TweedieRegressor(power=1.9, alpha=0.1, solver='newton-cholesky')\n",
    "glm_pure_premium.fit(\n",
    "    X_train, df_train[\"PurePremium\"], sample_weight=df_train[\"Exposure\"]\n",
    ")\n",
    "\n",
    "tweedie_powers = [1.5, 1.7, 1.8, 1.9, 1.99, 1.999, 1.9999]\n",
    "\n",
    "scores_glm_pure_premium = score_estimator(\n",
    "    glm_pure_premium,\n",
    "    X_train,\n",
    "    X_test,\n",
    "    df_train,\n",
    "    df_test,\n",
    "    target=\"PurePremium\",\n",
    "    weights=\"Exposure\",\n",
    "    tweedie_powers=tweedie_powers,\n",
    ")\n",
    "\n",
    "scores = pd.concat(\n",
    "    [scores_glm_pure_premium],\n",
    "    axis=1,\n",
    "    sort=True,\n",
    "    keys=(\"TweedieRegressor\"),\n",
    ")\n",
    "print(\"Evaluation of the Product Model and the Tweedie Regressor on target PurePremium\")\n",
    "with pd.option_context(\"display.expand_frame_repr\", False):\n",
    "    print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "016fd47e-1120-4ba3-912e-edfa664495e2",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:20.049012Z",
     "shell.execute_reply.started": "2024-08-06T09:15:20.046778Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [],
   "source": [
    "n_iter = glm_pure_premium.n_iter_\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e940c2bd-cf38-49b3-a61f-3f1a83e9f133",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:22.972444Z",
     "shell.execute_reply.started": "2024-08-06T09:15:20.050440Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The version of SecretFlow: 1.6.1b0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.10/subprocess.py:1796: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n",
      "  self.pid = _posixsubprocess.fork_exec(\n",
      "2024-08-06 09:15:22,016\tWARNING services.py:1996 -- WARNING: The object store is using /tmp instead of /dev/shm because /dev/shm has only 67104768 bytes available. This will harm performance! You may be able to free up space by deleting files in /dev/shm. If you are inside a Docker container, you can increase /dev/shm size by passing '--shm-size=2.26gb' to 'docker run' (or add it to the run_options list in a Ray cluster config). Make sure to set this to more than 30% of available RAM.\n",
      "2024-08-06 09:15:22,144\tINFO worker.py:1724 -- Started a local Ray instance.\n"
     ]
    }
   ],
   "source": [
    "import secretflow as sf\n",
    "\n",
    "print('The version of SecretFlow: {}'.format(sf.__version__))\n",
    "\n",
    "sf.shutdown()\n",
    "\n",
    "sf.init(['alice', 'bob'], address='local')\n",
    "\n",
    "alice, bob = sf.PYU('alice'), sf.PYU('bob')\n",
    "spu = sf.SPU(\n",
    "    sf.utils.testing.cluster_def(\n",
    "        ['alice', 'bob'],\n",
    "        {\"protocol\": \"REF2K\", \"field\": \"FM128\", \"fxp_fraction_bits\": 40},\n",
    "    ),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b8ec6ca9-a5c0-4d2f-a46d-833b530c669a",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:24.192229Z",
     "shell.execute_reply.started": "2024-08-06T09:15:22.974469Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [],
   "source": [
    "from secretflow.data import FedNdarray, PartitionWay\n",
    "\n",
    "x, y = X_train, df_train[\"PurePremium\"]\n",
    "w = df_train[\"Exposure\"]\n",
    "\n",
    "\n",
    "def x_to_vdata(x):\n",
    "    x = x.todense()\n",
    "    v_data = FedNdarray(\n",
    "        partitions={\n",
    "            alice: alice(lambda: x[:, :15])(),\n",
    "            bob: bob(lambda: x[:, 15:])(),\n",
    "        },\n",
    "        partition_way=PartitionWay.VERTICAL,\n",
    "    )\n",
    "    return v_data\n",
    "\n",
    "\n",
    "v_data = x_to_vdata(x)\n",
    "\n",
    "label_data = FedNdarray(\n",
    "    partitions={alice: alice(lambda: y.values)()},\n",
    "    partition_way=PartitionWay.VERTICAL,\n",
    ")\n",
    "\n",
    "sample_weight = FedNdarray(\n",
    "    partitions={alice: alice(lambda: w.values)()},\n",
    "    partition_way=PartitionWay.VERTICAL,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c06a4f96-33ae-48ae-82aa-fc30f52953bc",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:24.227640Z",
     "shell.execute_reply.started": "2024-08-06T09:15:24.196275Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [],
   "source": [
    "from secretflow.device.driver import reveal\n",
    "from secretflow.ml.linear.ss_glm.core import get_dist\n",
    "\n",
    "dist = 'Tweedie'\n",
    "ss_glm_power = 1.9\n",
    "\n",
    "\n",
    "class DirectRevealModel:\n",
    "    def __init__(self, model) -> None:\n",
    "        self.model = model\n",
    "\n",
    "    def predict(self, X):\n",
    "        vdata = x_to_vdata(X)\n",
    "        y = self.model.predict(vdata)\n",
    "        return reveal(y).reshape((-1,))\n",
    "\n",
    "    def score(self, X, y, sample_weight=None):\n",
    "        y = y.values\n",
    "        y_pred = self.predict(X)\n",
    "\n",
    "        constant = np.mean(y)\n",
    "        if sample_weight is not None:\n",
    "            constant *= sample_weight.shape[0] / np.sum(sample_weight)\n",
    "\n",
    "        # Missing factor of 2 in deviance cancels out.\n",
    "        deviance = get_dist(dist, 1, ss_glm_power).deviance(y_pred, y, None)\n",
    "        deviance_null = get_dist(dist, 1, ss_glm_power).deviance(\n",
    "            np.average(y, weights=sample_weight) + np.zeros(y.shape), y, None\n",
    "        )\n",
    "        return 1 - (deviance + constant) / (deviance_null + constant)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8eea5348-f070-4bbd-a4ce-20d2cb093396",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:45.226758Z",
     "shell.execute_reply.started": "2024-08-06T09:15:24.230529Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:jax._src.xla_bridge:Unable to initialize backend 'cuda': \n",
      "INFO:jax._src.xla_bridge:Unable to initialize backend 'rocm': module 'jaxlib.xla_extension' has no attribute 'GpuAllocatorConfig'\n",
      "INFO:jax._src.xla_bridge:Unable to initialize backend 'tpu': INTERNAL: Failed to open libtpu.so: libtpu.so: cannot open shared object file: No such file or directory\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls updating weights...\n",
      "INFO:root:epoch 1 train times: 12.741058826446533s\n",
      "INFO:root:epoch 1 validation time cost: 1.1556816101074219\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls calculating partials...\n",
      "INFO:root:irls updating weights...\n",
      "INFO:root:epoch 2 train times: 3.0125348567962646s\n",
      "INFO:root:epoch 2 validation time cost: 0.19248199462890625\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "from secretflow.ml.linear.ss_glm import SSGLM\n",
    "\n",
    "model = SSGLM(spu)\n",
    "\n",
    "ss_glm_power = 1.9\n",
    "start = time.time()\n",
    "model.fit_irls(\n",
    "    v_data,\n",
    "    label_data,\n",
    "    None,\n",
    "    sample_weight,\n",
    "    2,\n",
    "    'Log',\n",
    "    'Tweedie',\n",
    "    ss_glm_power,\n",
    "    l2_lambda=0.1,\n",
    "    infeed_batch_size_limit=10000000,\n",
    "    fraction_of_validation_set=0.2,\n",
    "    stopping_rounds=2,\n",
    "    stopping_metric='deviance',\n",
    "    stopping_tolerance=0.001,\n",
    ")\n",
    "\n",
    "wrapped_model = DirectRevealModel(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a65669dc-902b-4e1b-8323-5985996670dc",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:15:45.252950Z",
     "shell.execute_reply.started": "2024-08-06T09:15:45.229018Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.41737620e-03],\n",
       "       [ 1.89247951e-01],\n",
       "       [ 5.89404225e-01],\n",
       "       [-1.12057433e-01],\n",
       "       [-1.91406030e-02],\n",
       "       [ 3.76760662e-01],\n",
       "       [-1.78643540e-01],\n",
       "       [ 7.97630921e-02],\n",
       "       [ 6.48671240e-02],\n",
       "       [-6.64113700e-01],\n",
       "       [ 3.58013153e-01],\n",
       "       [-9.31122601e-01],\n",
       "       [-4.33228731e-01],\n",
       "       [-4.25371081e-02],\n",
       "       [ 2.17337251e-01],\n",
       "       [ 7.76762515e-02],\n",
       "       [ 4.42767352e-01],\n",
       "       [ 3.43533568e-02],\n",
       "       [ 3.02862942e-01],\n",
       "       [ 5.23277819e-01],\n",
       "       [-2.29991794e-01],\n",
       "       [ 7.61393070e-01],\n",
       "       [ 6.59996927e-01],\n",
       "       [-4.33898538e-01],\n",
       "       [ 2.84868866e-01],\n",
       "       [-5.70980132e-01],\n",
       "       [ 1.28583670e-01],\n",
       "       [ 1.76516414e-01],\n",
       "       [ 1.00819051e-01],\n",
       "       [-7.33127445e-02],\n",
       "       [-6.76781118e-01],\n",
       "       [-2.44344637e-01],\n",
       "       [-5.11223078e-01],\n",
       "       [-5.94665557e-02],\n",
       "       [-1.31659880e-01],\n",
       "       [-4.47464734e-01],\n",
       "       [ 2.72815317e-01],\n",
       "       [ 4.01874661e-01],\n",
       "       [ 4.99068677e-01],\n",
       "       [ 2.64340192e-01],\n",
       "       [ 9.79155481e-01],\n",
       "       [-2.00659722e-01],\n",
       "       [-6.43545151e-01],\n",
       "       [ 5.36597788e-01],\n",
       "       [ 4.61995959e-01],\n",
       "       [-2.98327297e-01],\n",
       "       [ 4.34670162e+00],\n",
       "       [-2.82204330e-01],\n",
       "       [-5.53125978e-01],\n",
       "       [-7.98071548e-02],\n",
       "       [ 5.19124866e-01],\n",
       "       [-2.48122420e-02],\n",
       "       [-4.89941418e-01],\n",
       "       [-5.83527446e-01],\n",
       "       [-1.33512449e+00],\n",
       "       [-1.00906539e+00],\n",
       "       [-4.37948197e-01],\n",
       "       [-1.15588337e-01],\n",
       "       [-3.31586003e-01],\n",
       "       [ 1.25630423e-01],\n",
       "       [-8.83829653e-01],\n",
       "       [ 4.85097356e-02],\n",
       "       [ 3.30506265e-01],\n",
       "       [-1.05094850e+00],\n",
       "       [-2.47250766e-01],\n",
       "       [ 4.02382225e-01],\n",
       "       [ 2.50227308e+00],\n",
       "       [ 8.51446509e-01],\n",
       "       [ 4.72118527e-01],\n",
       "       [ 2.13906214e-01],\n",
       "       [ 1.14421956e-01],\n",
       "       [-1.06497288e-01],\n",
       "       [-1.41345167e+00],\n",
       "       [ 8.99265930e-02],\n",
       "       [ 5.71904898e-01],\n",
       "       [-1.63847847e+01]], dtype=float32)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reveal(model.spu_w)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d7033237-f6cb-4518-8b64-147eb95b294a",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:17:21.461152Z",
     "shell.execute_reply.started": "2024-08-06T09:15:45.254836Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2024-08-06 09:15:52.084] [info] [thread_pool.cc:30] Create a fixed thread pool with size 19\n"
     ]
    }
   ],
   "source": [
    "tweedie_powers = [1.5, 1.7, 1.8, 1.9, 1.99, 1.999, 1.9999]\n",
    "\n",
    "scores_ss_glm_pure_premium = score_estimator(\n",
    "    wrapped_model,\n",
    "    X_train,\n",
    "    X_test,\n",
    "    df_train,\n",
    "    df_test,\n",
    "    target=\"PurePremium\",\n",
    "    weights=\"Exposure\",\n",
    "    tweedie_powers=tweedie_powers,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "34e55522-886f-42f5-9785-d59c0ca090c0",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:17:21.492446Z",
     "shell.execute_reply.started": "2024-08-06T09:17:21.462803Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluation of the Tweedie Regressor and SS GLM on target PurePremium\n",
      "                          TweedieRegressor                     SSGLMRegressor                       \n",
      "subset                               train          test                train                   test\n",
      "metric                                                                                              \n",
      "D² explained                  1.690000e-02  1.420000e-02         -0.008159399          -0.0010857582\n",
      "mean Tweedie dev p=1.5000     7.640770e+01  7.640880e+01           104.685472             107.990646\n",
      "mean Tweedie dev p=1.7000     3.682880e+01  3.692270e+01            42.253085              42.932537\n",
      "mean Tweedie dev p=1.8000     3.037600e+01  3.045390e+01            32.950962              33.348731\n",
      "mean Tweedie dev p=1.9000     3.382120e+01  3.387830e+01            35.103795              35.341118\n",
      "mean Tweedie dev p=1.9900     2.015347e+02  2.015587e+02           202.247704             202.379453\n",
      "mean Tweedie dev p=1.9990     1.914538e+03  1.914387e+03          1915.211589             1915.16268\n",
      "mean Tweedie dev p=1.9999     1.904747e+04  1.904558e+04          19048.14074           19046.352948\n",
      "mean abs. error               2.739865e+02  2.731249e+02          4743.568604            32958.44976\n",
      "mean squared error            3.295505e+07  3.213056e+07  781756613944.103027  276791909504836.03125\n",
      "MultiIndex([('TweedieRegressor', 'train'),\n",
      "            ('TweedieRegressor',  'test'),\n",
      "            (  'SSGLMRegressor', 'train'),\n",
      "            (  'SSGLMRegressor',  'test')],\n",
      "           names=[None, 'subset'])\n",
      "明密文训练数据比较\n",
      "metric\n",
      "D² explained                   -1.4828047\n",
      "mean Tweedie dev p=1.5000        0.370091\n",
      "mean Tweedie dev p=1.7000        0.147284\n",
      "mean Tweedie dev p=1.8000         0.08477\n",
      "mean Tweedie dev p=1.9000        0.037923\n",
      "mean Tweedie dev p=1.9900        0.003538\n",
      "mean Tweedie dev p=1.9990        0.000352\n",
      "mean Tweedie dev p=1.9999        0.000035\n",
      "mean abs. error                 16.313147\n",
      "mean squared error           23720.903742\n",
      "dtype: object\n",
      "明密文测试数据比较\n",
      "metric\n",
      "D² explained                     -1.0764618\n",
      "mean Tweedie dev p=1.5000          0.413327\n",
      "mean Tweedie dev p=1.7000          0.162768\n",
      "mean Tweedie dev p=1.8000          0.095056\n",
      "mean Tweedie dev p=1.9000          0.043179\n",
      "mean Tweedie dev p=1.9900          0.004072\n",
      "mean Tweedie dev p=1.9990          0.000405\n",
      "mean Tweedie dev p=1.9999          0.000041\n",
      "mean abs. error                  119.671714\n",
      "mean squared error           8614599.208325\n",
      "dtype: object\n"
     ]
    }
   ],
   "source": [
    "scores = pd.concat(\n",
    "    [scores_glm_pure_premium, scores_ss_glm_pure_premium],\n",
    "    axis=1,\n",
    "    sort=True,\n",
    "    keys=(\"TweedieRegressor\", \"SSGLMRegressor\"),\n",
    ")\n",
    "print(\"Evaluation of the Tweedie Regressor and SS GLM on target PurePremium\")\n",
    "with pd.option_context(\"display.expand_frame_repr\", False):\n",
    "    print(scores)\n",
    "    # print(scores.shape) # (10, 4)\n",
    "    print(scores.columns)\n",
    "    # print(scores['TweedieRegressor', 'train'])\n",
    "    train_plain = scores['TweedieRegressor', 'train']\n",
    "    train_cipher = scores['SSGLMRegressor', 'train']\n",
    "    train_cmp = (train_cipher - train_plain) / train_plain\n",
    "    print(f\"明密文训练数据比较\\n{train_cmp}\")\n",
    "    test_plain = scores['TweedieRegressor', 'test']\n",
    "    test_cipher = scores['SSGLMRegressor', 'test']\n",
    "    train_cmp = (test_cipher - test_plain) / test_plain\n",
    "    print(f\"明密文测试数据比较\\n{train_cmp}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b3c51360-d227-48a8-8d90-d55adc7af529",
   "metadata": {
    "execution": {
     "parties": "[\"alice\"]",
     "shell.execute_reply.end": "2024-08-06T09:17:29.808023Z",
     "shell.execute_reply.started": "2024-08-06T09:17:21.495008Z",
     "to_execute": "2024-08-06T09:14:42.953Z"
    },
    "libroFormatter": "formatter-string"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "subset                                 train          test\n",
      "observed                        3.917618e+07  1.299546e+07\n",
      "predicted, tweedie, power=1.90  3.951751e+07  1.325198e+07\n",
      "predicted, ss glm, power=1.90   1.250159e+09  2.951692e+09\n"
     ]
    }
   ],
   "source": [
    "res = []\n",
    "for subset_label, x, df in [\n",
    "    (\"train\", X_train, df_train),\n",
    "    (\"test\", X_test, df_test),\n",
    "]:\n",
    "    exposure = df[\"Exposure\"].values\n",
    "    res.append(\n",
    "        {\n",
    "            \"subset\": subset_label,\n",
    "            \"observed\": df[\"ClaimAmount\"].values.sum(),\n",
    "            \"predicted, tweedie, power=%.2f\" % glm_pure_premium.power: np.sum(\n",
    "                exposure * glm_pure_premium.predict(x)\n",
    "            ),\n",
    "            \"predicted, ss glm, power=%.2f\" % ss_glm_power: np.sum(\n",
    "                exposure * wrapped_model.predict(x)\n",
    "            ),\n",
    "        }\n",
    "    )\n",
    "print(pd.DataFrame(res).set_index(\"subset\").T)\n"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}