COGMOD-HWI/HW3/Q5/Q5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import stan\n",
    "import arviz as az\n",
    "\n",
    "# stan problems\n",
    "import nest_asyncio\n",
    "nest_asyncio.apply()\n",
    "\n",
    "np.random.seed(42)\n",
    "N = 100\n",
    "alpha = 2.3\n",
    "beta = 4.0\n",
    "sigma = 2.0\n",
    "x = np.random.normal(size=N)\n",
    "y = alpha + beta * x + sigma * np.random.normal(size=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_0 = 3\n",
    "sigma_0 = 3\n",
    "sample_mean = np.mean(y)\n",
    "n = N\n",
    "\n",
    "# posterior\n",
    "sigma_sq = sigma ** 2\n",
    "sigma_0_sq = sigma_0 ** 2\n",
    "\n",
    "mu_post = ((n / sigma_sq) * sample_mean + (1 / sigma_0_sq) * mu_0) / (n / sigma_sq + 1 / sigma_0_sq)\n",
    "sigma_post_sq = 1 / (n / sigma_sq + 1 / sigma_0_sq)\n",
    "sigma_post = np.sqrt(sigma_post_sq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHHCAYAAABXx+fLAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAATHdJREFUeJzt3XlcFWX/P/7XAeSwH0CBA4aAopQmYBqEa9ZRJDNpMTTvWO409XbJkEr85lrdlKZixS1WKlLuZVpqJJFoKu7RokmiIBiLS8IRVBDO9fvDn/PpxCL7Aeb1fDzm8bnPNddc856jeV6fmWtmFEIIASIiIiIZMTJ0AUREREQtjQGIiIiIZIcBiIiIiGSHAYiIiIhkhwGIiIiIZIcBiIiIiGSHAYiIiIhkhwGIiIiIZIcBiIiIiGSHAYiIahUeHg53d3dDl0HVcHd3R3h4uKHLIGqTGICI2rGEhAQoFAppMTMzQ48ePTBt2jQUFhYaurwW5+7urvd9ODo6YtCgQfjqq6+aZX+7d+/GggULmmVsImocBd8FRtR+JSQkICIiAosWLYKHhwdu3bqFAwcO4LPPPoObmxt+++03WFhY1DrG7du3odPpoFQqW6jq5uPu7g47OzvMmjULAJCXl4dVq1bh/PnzWLlyJSZPntyk+5s2bRri4uLQXP/MlpWVwcjICB06dGiW8YnaMxNDF0BEzS8oKAj9+vUDAEyYMAEdO3bEsmXLsGPHDowbN67abUpLS2FpadmkP646nQ7l5eUwMzNrsjHrq3PnzvjXv/4lfQ4NDYWnpyeWL1/e5AGoOQghcOvWLZibmzdpKL116xZMTU1hZMQLAyQP/JtOJEOPPfYYACArKwvAnXk+VlZWOHfuHJ544glYW1tj/Pjx0rp/zgEqLS3FrFmz4OrqCqVSCS8vL7z//vtVznQoFApMmzYN69evR69evaBUKpGUlFRtTU8++SS6du1a7bqAgAApwAFAcnIyBg4cCFtbW1hZWcHLywtz5sxp0HehVqvxwAMPSN8FAPz0008ICgqCjY0NrKys8Pjjj+Pw4cN6292+fRsLFy5E9+7dYWZmho4dO2LgwIFITk4GcOd7i4uLk76Hu8tdOp0OsbGx6NWrF8zMzODk5IRJkybh2rVrevtxd3fHk08+ie+++w79+vWDubk5Vq1aJa375xyg8+fPY8yYMbC3t4eFhQUeeeQR7Nq1S69PamoqFAoFNm3ahDfffBOdO3eGhYUFtFptg75DoraIZ4CIZOjcuXMAgI4dO0ptFRUVCAwMxMCBA/H+++/XeGlMCIGnnnoKe/fuxUsvvQRfX1989913eO211/Dnn39i+fLlev1/+OEHbNmyBdOmTUOnTp1qnFAdEhKC0NBQHDt2DA8//LDUfuHCBRw+fBhLliwBAJw6dQpPPvkkvL29sWjRIiiVSmRmZuLgwYMN+i5u376N3Nxc6bs4deoUBg0aBBsbG7z++uvo0KEDVq1ahUcffRT79u2Dv78/AGDBggWIiYnBhAkT4OfnB61Wi+PHj+PkyZMYNmwYJk2ahLy8PCQnJ+Ozzz6rst9JkyZJlyhnzJiBrKwsfPTRR/jpp59w8OBBvTNvGRkZGDduHCZNmoSJEyfCy8ur2mMpLCxE//79cePGDcyYMQMdO3bEunXr8NRTT+GLL77A008/rdf/rbfegqmpKaKiolBWVgZTU9MGfYdEbZIgonZr7dq1AoD4/vvvxeXLl0Vubq7YtGmT6NixozA3NxcXL14UQggRFhYmAIjZs2dXGSMsLEy4ublJn7dv3y4AiLfffluv33PPPScUCoXIzMyU2gAIIyMjcerUqXvWWlxcLJRKpZg1a5Ze++LFi4VCoRAXLlwQQgixfPlyAUBcvny5zt/DXW5ubmL48OHi8uXL4vLly+Lnn38WY8eOFQDE9OnThRBCBAcHC1NTU3Hu3Dlpu7y8PGFtbS0GDx4stfn4+IiRI0fWur+pU6eK6v6Z/fHHHwUAsX79er32pKSkKu1ubm4CgEhKSqr2eMLCwqTPM2fOFADEjz/+KLVdv35deHh4CHd3d1FZWSmEEGLv3r0CgOjatau4ceNGrcdA1F7xEhiRDGg0Gjg4OMDV1RVjx46FlZUVvvrqK3Tu3Fmv35QpU+451u7du2FsbIwZM2botc+aNQtCCHz77bd67UOGDEHPnj3vOa6NjQ2CgoKwZcsWvUtpmzdvxiOPPIIuXboAAGxtbQEAO3bsgE6nu+e4/7Rnzx44ODjAwcEBPj4+2Lp1K1588UW89957qKysxJ49exAcHKx3Oc7Z2RkvvPACDhw4IF0msrW1xalTp3D27Nl617B161aoVCoMGzYMV65ckZa+ffvCysoKe/fu1evv4eGBwMDAe467e/du+Pn5YeDAgVKblZUVXn75ZWRnZ+P06dN6/cPCwmBubl7v+onaAwYgIhmIi4tDcnIy9u7di9OnT+P8+fNVflBNTExw33333XOsCxcuwMXFBdbW1nrtDzzwgLT+7zw8POpcZ0hICHJzc5GWlgbgzqW6EydOICQkRK/PgAEDMGHCBDg5OWHs2LHYsmVLncOQv78/kpOT8f333+PQoUO4cuUKEhMTYW5ujsuXL+PGjRvVXmJ64IEHoNPpkJubCwBYtGgRioqK0KNHD/Tu3RuvvfYafvnllzrVcPbsWRQXF8PR0VEKY3eXkpISXLp0Sa9/Xb/DCxcu1Fj73fUNGZeoPeIcICIZ8PPz05tEXB2lUtksdwDV5wzDqFGjYGFhgS1btqB///7YsmULjIyMMGbMGL3x9u/fj71792LXrl1ISkrC5s2b8dhjj2HPnj0wNjaudR+dOnWCRqNp8PHcNXjwYJw7dw47duzAnj178Omnn2L58uWIj4/HhAkTat1Wp9PB0dER69evr3a9g4OD3ufmOkvDsz8kZzwDRET14ubmhry8PFy/fl2v/cyZM9L6hrK0tMSTTz6JrVu3QqfTYfPmzRg0aBBcXFz0+hkZGeHxxx/HsmXLcPr0abzzzjv44Ycfqlw6qi8HBwdYWFggIyOjyrozZ87AyMgIrq6uUpu9vT0iIiKwceNG5ObmwtvbW+/Bh3+/6+vvunXrhqtXr2LAgAHQaDRVFh8fnwbV7+bmVmPtd9cT0R0MQERUL0888QQqKyvx0Ucf6bUvX74cCoUCQUFBjRo/JCQEeXl5+PTTT/Hzzz/rXf4CgL/++qvKNr6+vgDuPBiwMYyNjTF8+HDs2LED2dnZUnthYSE2bNiAgQMHwsbGBgBw9epVvW2trKzg6empV4OlpSUAoKioSK/v888/j8rKSrz11ltVaqioqKjSv66eeOIJHD16VLqECNx5ZMHHH38Md3f3Os3FIpILXgIjonoZNWoUhg4div/3//4fsrOz4ePjgz179mDHjh2YOXMmunXr1qjx7z6HKCoqCsbGxnj22Wf11i9atAj79+/HyJEj4ebmhkuXLuF///sf7rvvPr3Jvw319ttvS88Z+s9//gMTExOsWrUKZWVlWLx4sdSvZ8+eePTRR9G3b1/Y29vj+PHj+OKLLzBt2jSpT9++fQEAM2bMQGBgIIyNjTF27FgMGTIEkyZNQkxMDNLT0zF8+HB06NABZ8+exdatW7FixQo899xz9a599uzZ2LhxI4KCgjBjxgzY29tj3bp1yMrKwpdffsmHHBL9naFvQyOi5nP3Nvhjx47V2i8sLExYWlrWuO7vt8ELcefW6ldffVW4uLiIDh06iO7du4slS5YInU6n1w+AmDp1ar3rHj9+vAAgNBpNlXUpKSli9OjRwsXFRZiamgoXFxcxbtw48ccff9xzXDc3t3veui6EECdPnhSBgYHCyspKWFhYiKFDh4pDhw7p9Xn77beFn5+fsLW1Febm5uL+++8X77zzjigvL5f6VFRUiOnTpwsHBwehUCiq3BL/8ccfi759+wpzc3NhbW0tevfuLV5//XWRl5dXp5r/eRu8EEKcO3dOPPfcc8LW1laYmZkJPz8/sXPnTr0+d2+D37p16z2/C6L2iu8CIyIiItnh+VAiIiKSHQYgIiIikh0GICIiIpIdBiAiIiKSHQY
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prior_samples = np.random.normal(mu_0, sigma_0, 10000)\n",
    "posterior_samples = np.random.normal(mu_post, sigma_post, 10000)\n",
    "\n",
    "plt.hist(prior_samples, bins=50, density=True, alpha=0.5, label='Prior')\n",
    "plt.hist(posterior_samples, bins=50, density=True, alpha=0.5, label='Posterior')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('Average Height (meters)')\n",
    "plt.ylabel('Density')\n",
    "plt.title('Prior vs Posterior')\n",
    "plt.savefig('Q5.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stanCode = \"\"\"\n",
    "data {\n",
    "  int<lower=0> N;  // observations\n",
    "  vector[N] x;     // predictor\n",
    "  vector[N] y;     // response\n",
    "}\n",
    "parameters {\n",
    "  real alpha;      // intercept\n",
    "  real beta;       // slope\n",
    "  real<lower=0> sigma;  // noise standard deviation\n",
    "}\n",
    "model {\n",
    "  sigma ~ inv_gamma(1, 1);   // noise\n",
    "  alpha ~ normal(0, 10);     // intercept\n",
    "  beta ~ normal(0, 10);      // slope\n",
    "  y ~ normal(alpha + beta * x, sigma);\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "data = {\n",
    "    'N': N,\n",
    "    'x': x,\n",
    "    'y': y\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Building: found in cache, done.Sampling:   0%\n",
      "Sampling:  25% (3000/12000)\n",
      "Sampling:  50% (6000/12000)\n",
      "Sampling:  75% (9000/12000)\n",
      "Sampling: 100% (12000/12000)\n",
      "Sampling: 100% (12000/12000), done.\n",
      "Messages received during sampling:\n",
      "  Gradient evaluation took 4.2e-05 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:\n",
      "  Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/httpstan_lipc1ba5/model_oj7ef663.stan', line 16, column 2 to column 38)\n",
      "  If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,\n",
      "  but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.\n",
      "  Gradient evaluation took 2.8e-05 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Gradient evaluation took 2.2e-05 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.\n",
      "  Adjust your expectations accordingly!\n",
      "  Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:\n",
      "  Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/httpstan_lipc1ba5/model_oj7ef663.stan', line 16, column 2 to column 38)\n",
      "  If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,\n",
      "  but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.\n",
      "  Gradient evaluation took 3.8e-05 seconds\n",
      "  1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.\n",
      "  Adjust your expectations accordingly!\n"
     ]
    },
    {
     "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>alpha</th>\n",
       "      <td>2.315</td>\n",
       "      <td>0.193</td>\n",
       "      <td>1.946</td>\n",
       "      <td>2.672</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>7672.0</td>\n",
       "      <td>5912.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>3.716</td>\n",
       "      <td>0.214</td>\n",
       "      <td>3.305</td>\n",
       "      <td>4.108</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>6876.0</td>\n",
       "      <td>5564.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.910</td>\n",
       "      <td>0.140</td>\n",
       "      <td>1.662</td>\n",
       "      <td>2.175</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>7423.0</td>\n",
       "      <td>5615.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",
       "alpha  2.315  0.193   1.946    2.672      0.002    0.002    7672.0    5912.0   \n",
       "beta   3.716  0.214   3.305    4.108      0.003    0.002    6876.0    5564.0   \n",
       "sigma  1.910  0.140   1.662    2.175      0.002    0.002    7423.0    5615.0   \n",
       "\n",
       "       r_hat  \n",
       "alpha    1.0  \n",
       "beta     1.0  \n",
       "sigma    1.0  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = stan.build(stanCode, data=data)\n",
    "fit = model.sample(num_samples=2000, num_chains=4)\n",
    "\n",
    "az.summary(az.from_pystan(fit))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}