ensemble fitting worked with gumbel and normalgit add .1

ihmeuw-msca · Aug 13, 2024 · 825f375 · 825f375
1 parent 1ac70f1
commit 825f375
Showing 1 changed file with 58 additions and 27 deletions.
diff --git a/simulations.ipynb b/simulations.ipynb
@@ -6,41 +6,72 @@
    "metadata": {},
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[-1.88951825 -1.67562088 -1.65214474 -1.62064284 -1.562736   -1.50039701\n",
-      " -1.46558166 -1.41107123 -1.33930541 -1.33348901 -1.22021545 -1.16896446\n",
-      " -1.11823493 -1.11364725 -1.0355334  -0.99155624 -0.9322802  -0.91862877\n",
-      " -0.89289741 -0.87632799 -0.85337447 -0.82952583 -0.80436763 -0.77718609\n",
-      " -0.77606566 -0.74560126 -0.6891231  -0.66738905 -0.65431026 -0.60114386\n",
-      " -0.51794973 -0.51688887 -0.4953134  -0.44496799 -0.40446697 -0.39829545\n",
-      " -0.39799426 -0.31320473 -0.284331   -0.28041655 -0.22160632 -0.17987713\n",
-      " -0.17326124 -0.11955245 -0.08257138 -0.00575042  0.01167684  0.03075504\n",
-      "  0.05233829  0.06073081  0.08723611  0.09763968  0.10284089  0.13299994\n",
-      "  0.1459378   0.25129287  0.2637708   0.26529311  0.28513058  0.3164859\n",
-      "  0.35035412  0.35802117  0.36623514  0.3697894   0.40023227  0.41085433\n",
-      "  0.48258617  0.48964917  0.50151693  0.50449019  0.51650857  0.53792918\n",
-      "  0.56122631  0.57075674  0.57988314  0.5942923   0.62610387  0.6539129\n",
-      "  0.66528933  0.73613501  0.74289003  0.80039489  0.80555165  0.81712376\n",
-      "  0.94335029  0.99278054  1.01699842  1.02001017  1.05444052  1.07084284\n",
-      "  1.09992702  1.16661601  1.17901527  1.20899573  1.2415302   1.3356225\n",
-      "  1.4081054   1.61840104  1.83143626  1.85685823]\n",
-      "<ensemble.ensemble_model.EnsembleResult object at 0x16223f020>\n"
-     ]
+     "data": {
+      "text/plain": [
+       "array([-0.13232203,  0.43740831, -1.17082002, -0.23752026, -0.29514267,\n",
+       "       -0.12570246, -0.0199769 ,  0.584793  , -1.23747854, -0.87275655])"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
     "import numpy as np\n",
-    "import scipy.stats\n",
+    "import scipy\n",
     "import matplotlib.pyplot as plt\n",
     "from ensemble.distributions import distribution_dict\n",
     "from ensemble.ensemble_model import EnsembleFitter\n",
     "\n",
+    "def ensemble_cdf(x, distributions, weights, mean, variance):\n",
+    "  my_objs = []\n",
+    "  for distribution in distributions:\n",
+    "    my_objs.append(distribution_dict[distribution](mean, variance))\n",
+    "  return sum(weight * distribution.cdf(x) for distribution, weight in zip(my_objs, weights))\n",
     "\n",
-    "std_norm = distribution_dict[\"normal\"](0, 1).rvs(100)\n",
+    "def ppf_to_solve(x, q):\n",
+    "  return ensemble_cdf(x, [\"normal\", \"gumbel\"], [0.7, 0.3], 0, 1) - q\n",
+    "\n",
+    "def ppf_single(q):\n",
+    "  factor = 10.\n",
+    "  left = -factor\n",
+    "  right = factor\n",
+    "\n",
+    "  while ppf_to_solve(left, q) > 0:\n",
+    "    left, right = left * factor, left\n",
+    "\n",
+    "  while ppf_to_solve(right, q) < 0:\n",
+    "    left, right = right, right * factor\n",
+    "\n",
+    "  return scipy.optimize.brentq(ppf_to_solve, left, right, args=q)\n",
+    "\n",
+    "def ensemble_rvs(size):\n",
+    "  ppf_vec = np.vectorize(ppf_single, otypes=\"d\")\n",
+    "  unif_samp = scipy.stats.uniform.rvs(size=size)\n",
+    "  return ppf_vec(unif_samp)\n",
+    "\n",
+    "ensemble_rvs(10)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[0.77402318 0.2793297 ]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# std_norm = distribution_dict[\"normal\"](0, 1).rvs(100)\n",
+    "temp = ensemble_rvs(10)\n",
     "model = EnsembleFitter([\"normal\", \"gumbel\"], None)\n",
-    "res = model.fit(std_norm)\n",
+    "res = model.fit(temp)\n",
     "print(res.weights)\n",
     "\n",
     "\n",
@@ -80,7 +111,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -91,7 +122,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {