{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cf0c31f6-2c6c-4021-a7f7-974b37bbd632",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc as pm\n",
    "import numpy as np\n",
    "import arviz as az\n",
    "\n",
    "%load_ext lab_black"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "671c3bbc-f380-4830-a3b8-3f68bf359987",
   "metadata": {},
   "source": [
    "# 4. Priors as Hidden Mixtures*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd8abdb6-bebe-40bf-8fc7-626d5db0dead",
   "metadata": {},
   "source": [
    "Adapted from [Unit 7: tasmixture.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit7/tasmixture.odc) and [Jeremymus.odc](https://raw.githubusercontent.com/areding/6420-pymc/main/original_examples/Codes4Unit7/Jeremymus.odc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "800d9106-eec5-4e4c-9bc3-cf3271621145",
   "metadata": {},
   "source": [
    "## Student T likelihood with Normal prior vs Normal likelihood with Gamma prior\n",
    "This example demonstrates the equivalence of a couple methods. Variable names from original problem are retained right now, but they are kind of a mess. Will fix later.\n",
    "\n",
    "Confusingly, the Student's T distribution parameter sigma in PyMC is equivalent to the tau parameter in the BUGS parameterization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0cb127c8-7de7-4bef-aa33-20acefd71516",
   "metadata": {
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [mu1, prec, mu2]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:01&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 1 seconds.\n"
     ]
    }
   ],
   "source": [
    "df = 6\n",
    "x = y = 10\n",
    "mu0 = 6\n",
    "tau1 = 10\n",
    "tau = 0.4\n",
    "a = df / 2\n",
    "b = df / (2 * tau)\n",
    "\n",
    "with pm.Model() as m1:\n",
    "    mu1 = pm.StudentT(\"mu1\", nu=df, mu=mu0, sigma=tau)\n",
    "    pm.Normal(\"X\", mu1, tau=tau1, observed=x)\n",
    "\n",
    "    prec = pm.Gamma(\"prec\", a, b)\n",
    "    mu2 = pm.Normal(\"mu2\", mu0, tau=prec)\n",
    "    pm.Normal(\"Y\", mu2, tau=tau1, observed=y)\n",
    "\n",
    "    trace = pm.sample(3000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5b4da943-bf70-479b-b0c1-11518364598b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>mu1</th>\n",
       "      <td>9.824</td>\n",
       "      <td>0.322</td>\n",
       "      <td>9.204</td>\n",
       "      <td>10.408</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>15555.0</td>\n",
       "      <td>9439.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mu2</th>\n",
       "      <td>9.909</td>\n",
       "      <td>0.318</td>\n",
       "      <td>9.312</td>\n",
       "      <td>10.502</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>12622.0</td>\n",
       "      <td>9347.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>prec</th>\n",
       "      <td>0.233</td>\n",
       "      <td>0.126</td>\n",
       "      <td>0.036</td>\n",
       "      <td>0.460</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>14735.0</td>\n",
       "      <td>8701.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "mu1   9.824  0.322   9.204   10.408      0.003    0.002   15555.0    9439.0   \n",
       "mu2   9.909  0.318   9.312   10.502      0.003    0.002   12622.0    9347.0   \n",
       "prec  0.233  0.126   0.036    0.460      0.001    0.001   14735.0    8701.0   \n",
       "\n",
       "      r_hat  \n",
       "mu1     1.0  \n",
       "mu2     1.0  \n",
       "prec    1.0  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef87067b-0c44-4a72-b149-74a4a93ad021",
   "metadata": {},
   "source": [
    "## Jeremy's IQ - Normal prior as a scale mixture of uniforms\n",
    "\n",
    "Any symmetric unimodal distribution is a scale mixture of uniforms.\n",
    "\n",
    "$$y|\\mu, \\sigma^2  \\sim \\mathcal{N}(\\mu, \\sigma^2)$$\n",
    "\n",
    "is the same as\n",
    "\n",
    "$$y|\\mu, \\xi  \\sim U(\\mu - \\sigma d^{1/2}, \\mu + \\sigma d^{1/2})$$\n",
    "\n",
    "\n",
    "$$d \\sim Ga(3/2, 1/2)$$\n",
    "\n",
    "If\n",
    "$d \\sim Ga(3/2, s/2)$, $s \\lt 1$ heavy tails, $s\\gt 1$ light tails\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cff3ee14-5535-40f0-a7f3-5fb07248f65e",
   "metadata": {
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [d, theta]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:01&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 1 seconds.\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [d, theta]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:01&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 1 seconds.\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [d, theta]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "    /* Turns off some styling */\n",
       "    progress {\n",
       "        /* gets rid of default border in Firefox and Opera. */\n",
       "        border: none;\n",
       "        /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "        background-size: auto;\n",
       "    }\n",
       "    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:01&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 3_000 draw iterations (4_000 + 12_000 draws total) took 1 seconds.\n"
     ]
    }
   ],
   "source": [
    "precy = 0.0125\n",
    "precth = 0.0083333\n",
    "mu = 110\n",
    "s_list = [0.5, 1, 2]\n",
    "y = 98\n",
    "\n",
    "summaries = []\n",
    "\n",
    "for s in s_list:\n",
    "    beta = s / 2\n",
    "    with pm.Model() as m2:\n",
    "        d = pm.Gamma(\"d\", 1.5, beta)\n",
    "        lb = mu - (d / precth) ** 0.5\n",
    "        ub = mu + (d / precth) ** 0.5\n",
    "\n",
    "        theta = pm.Uniform(\"theta\", lb, ub)\n",
    "\n",
    "        pm.Normal(\"y\", mu=theta, tau=precy, observed=y)\n",
    "\n",
    "        trace = pm.sample(3000)\n",
    "\n",
    "    summaries.append(az.summary(trace, kind=\"stats\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f6d183ac-08b0-4825-8c2e-1abd8a30226d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "s=0.5\n",
      "          mean     sd  hdi_3%  hdi_97%\n",
      "d        5.214  4.249   0.043   12.921\n",
      "theta  101.093  7.655  87.282  115.913\n",
      "\n",
      "\n",
      "s=1\n",
      "          mean     sd  hdi_3%  hdi_97%\n",
      "d        2.838  2.287   0.008    6.886\n",
      "theta  102.950  6.902  90.222  116.055\n",
      "\n",
      "\n",
      "s=2\n",
      "          mean     sd  hdi_3%  hdi_97%\n",
      "d        1.512  1.192   0.014    3.715\n",
      "theta  104.845  5.870  93.745  115.580\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for s, summary in zip(s_list, summaries):\n",
    "    print(f\"{s=}\")\n",
    "    print(summary)\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "616d0040-668a-4d79-9fab-44a98c90ac76",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: Wed Mar 22 2023\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.11.0\n",
      "IPython version      : 8.9.0\n",
      "\n",
      "pytensor: 2.10.1\n",
      "\n",
      "pymc : 5.1.2\n",
      "numpy: 1.24.2\n",
      "arviz: 0.14.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -p pytensor"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}