2025-03-12 14:10:26 -04:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
2025-03-15 16:20:08 -04:00
"# Problem 5\n"
2025-03-12 14:10:26 -04:00
]
},
{
"cell_type": "code",
2025-03-15 16:20:08 -04:00
"execution_count": 2,
2025-03-12 14:10:26 -04:00
"metadata": {},
"outputs": [
{
2025-03-13 13:47:02 -04:00
"data": {
"text/plain": [
2025-03-15 16:20:08 -04:00
"<Axes: ylabel='Count'>"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAK1hJREFUeJzt3X1wVPWhxvFnw0tCbBLEhLxoePMlBIQgKGmoVigpIdfLBfVylREJFLHDEKumahsv79rG6gjUkgttR4gOF1FnFHtbJr0QBeokoITmapyYITRxQdjgomFJ0iSQnPtHh223JCEJm5zd/L6fmTPDOed3zj7nzE545pyzuw7LsiwBAAAYJMTuAAAAAH2NAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMM5AuwMEora2Np06dUoRERFyOBx2xwEAAF1gWZbOnz+vhIQEhYR0fo2HAtSOU6dOKTEx0e4YAACgB06cOKEbbrih0zEUoHZERERI+tsJjIyMtDkNAADoCo/Ho8TERO//452hALXj0m2vyMhIChAAAEGmK4+v8BA0AAAwDgUIAAAYhwIEAACMQwECAADGoQABAADjUIAAAIBxKEAAAMA4FCAAAGAcChAAADAOBQgAABiHAgQAAIxDAQIAAMahAAEAAOPwa/AArorT6ZTb7bY7RrdER0drxIgRdscAYCMKEIAeczqdSk5KUmNTk91RuiU8LEwVlZWUIMBgFCAAPeZ2u9XY1KQdyclKDg+3O06XVDQ2amFFhdxuNwUIMBgFCMBVSw4P1+SICLtjAECX8RA0AAAwDgUIAAAYhwIEAACMQwECAADGoQABAADjUIAAAIBxKEAAAMA4FCAAAGAcChAAADAOBQgAABjH1gKUl5enO+64QxERERo+fLjmzZunyspKnzFNTU1asWKFrrvuOn3rW9/S/fffr9ra2k73a1mWVq9erfj4eA0ZMkTp6ek6duxYbx4KAAAIIrYWoAMHDmjFihU6dOiQ9u7dqwsXLmjWrFlqaGjwjnnyySf1P//zP3r77bd14MABnTp1Svfdd1+n+33xxRf1yiuvaOvWrTp8+LCuueYaZWRkqCnIfrEaAAD0Dlt/DLWwsNBnvqCgQMOHD1dpaam++93v6ty5c3r11Ve1c+dOfe9735Mkbd++XcnJyTp06JC+/e1vX7ZPy7K0adMmrVy5UnPnzpUkvf7664qNjdXu3bv14IMP9v6BAQCAgBZQvwZ/7tw5SdKwYcMkSaWlpbpw4YLS09O9Y8aOHasRI0aopKSk3QJUXV0tl8vls01UVJRSU1NVUlLSbgFqbm5Wc3Ozd97j8fjtmIDucDqdcrvddsfosoqKCrsjAECPBEwBamtr0xNPPKHvfOc7uvXWWyVJLpdLgwcP1tChQ33GxsbGyuVytbufS8tjY2O7vE1eXp7WrVt3lUcAXB2n06nkpCQ1BuGt2uaWFrsjAEC3BEwBWrFihcrLy/Xhhx/2+Wvn5uYqJyfHO+/xeJSYmNjnOWA2t9utxqYm7UhOVnJ4uN1xumTP2bNaVVOjixcv2h0FALolIApQdna2fv/73+vgwYO64YYbvMvj4uLU0tKiuro6n6tAtbW1iouLa3dfl5bX1tYqPj7eZ5tJkya1u01oaKhCQ0Ov/kAAP0gOD9fkiAi7Y3RJRWOj3REAoEds/RSYZVnKzs7Wu+++q/fff1+jR4/2WT9lyhQNGjRIRUVF3mWVlZVyOp1KS0trd5+jR49WXFyczzYej0eHDx/ucBsAAGAWWwvQihUrtGPHDu3cuVMRERFyuVxyuVz661//KulvDy8vXbpUOTk5+uCDD1RaWqolS5YoLS3N5wHosWPH6t1335UkORwOPfHEE3r++ef1u9/9Tp9++qkWLVqkhIQEzZs3z47DBAAAAcbWW2BbtmyRJE2fPt1n+fbt27V48WJJ0saNGxUSEqL7779fzc3NysjI0H/913/5jK+srPR+gkySnnnmGTU0NOjRRx9VXV2d7rzzThUWFiosLKxXjwcAAAQHWwuQZVlXHBMWFqb8/Hzl5+d3eT8Oh0Pr16/X+vXrrzojAADof/gtMAAAYBwKEAAAMA4FCAAAGIcCBAAAjEMBAgAAxqEAAQAA41CAAACAcShAAADAOBQgAABgHAoQAAAwDgUIAAAYhwIEAACMQwECAADGoQABAADjUIAAAIBxKEAAAMA4FCAAAGAcChAAADAOBQgAABiHAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMM5AuwMAAK7M6XTK7XbbHaNboqOjNWLECLtjAO2iAAFAgHM6nUpOSlJjU5PdUbolPCxMFZWVlCAEJAoQAAQ4t9utxqYm7UhOVnJ4uN1xuqSisVELKyrkdrspQAhIFCAACBLJ4eGaHBFhdwygX+AhaAAAYBxbC9DBgwc1Z84cJSQkyOFwaPfu3T7rHQ5Hu9NLL73U4T7Xrl172fixY8f28pEAAIBgYmsBamhoUEpKivLz89tdf/r0aZ9p27Ztcjgcuv/++zvd7/jx4322+/DDD3sjPgAACFK2PgOUmZmpzMzMDtfHxcX5zL/33nuaMWOGxowZ0+l+Bw4ceNm2AAAAlwTNM0C1tbX6wx/+oKVLl15x7LFjx5SQkKAxY8booYcektPp7HR8c3OzPB6PzwQAAPqvoClAr732miIiInTfffd1Oi41NVUFBQUqLCzUli1bVF1drbvuukvnz5/vcJu8vDxFRUV5p8TERH/HBwAAASRoCtC2bdv00EMPKSwsrNNxmZmZmj9/viZOnKiMjAzt2bNHdXV1euuttzrcJjc3V+fOnfNOJ06c8Hd8AAAQQILie4D+9Kc/qbKyUm+++Wa3tx06dKhuueUWVVVVdTgmNDRUoaGhVxMRAAAEkaC4AvTqq69qypQpSklJ6fa29fX1On78uOLj43shGQAACEa2FqD6+nqVlZWprKxMklRdXa2ysjKfh5Y9Ho/efvttPfLII+3uY+bMmdq8ebN3/qmnntKBAwdUU1Oj4uJi3XvvvRowYIAWLFjQq8cCAACCh623wI4cOaIZM2Z453NyciRJWVlZKigokCTt2rVLlmV1WGCOHz/u8wvJJ0+e1IIFC3T27FnFxMTozjvv1KFDhxQTE9N7BwIAAIKKrQVo+vTpsiyr0zGPPvqoHn300Q7X19TU+Mzv2rXLH9EAAEA/FhTPAAEAAPgTBQgAABiHAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMA4FCAAAGIcCBAAAjEMBAgAAxqEAAQAA41CAAACAcShAAADAOBQgAABgHAoQAAAwzkC7AwCAHSoqKuyO0GXBlBUIFhQgAEY53dKiEEkLFy60O0q3Nbe02B0B6DcoQACMUnfxotok/XbUKE2+7jq743TJnrNntaqmRhcvXrQ7CtBvUIAAGClpyBBNjoiwO0aXVDQ22h0B6Hd4CBoAABiHAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMA4FCAAAGIcCBAAAjGNrATp48KDmzJmjhIQEORwO7d6922f94sWL5XA4fKbZs2dfcb/5+fkaNWqUwsLClJqaqo8++qiXjgAAAAQjWwtQQ0ODUlJSlJ+f3+GY2bNn6/Tp097pjTfe6HSfb775pnJycrRmzRodPXpUKSkpysjI0JkzZ/wdHwAABClbfww1MzNTmZmZnY4JDQ1VXFxcl/e5YcMGLVu2TEuWLJEkbd26VX/4wx+0bds2/fSnP72qvAAAoH8I+GeA9u/fr+HDhyspKUnLly/X2bNnOxzb0tKi0tJSpaene5eFhIQoPT1dJSUlHW7X3Nwsj8fjMwEAgP4roAvQ7Nmz9frrr6uoqEi/+MUvdODAAWVmZqq1tbXd8W63W62trYqNjfVZHhsbK5fL1eHr5OXlKSoqyjslJib69TgAAEBgsfUW2JU8+OCD3n9PmDBBEydO1I033qj9+/dr5syZfnud3Nxc5eTkeOc9Hg8lCACAfiygrwD9szFjxig
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
2025-03-13 13:47:02 -04:00
]
},
"metadata": {},
"output_type": "display_data"
2025-03-12 14:10:26 -04:00
}
],
"source": [
"import numpy as np\n",
2025-03-15 16:20:08 -04:00
"import stan\n",
"import arviz as az\n",
"import nest_asyncio\n",
"import seaborn as sns\n",
2025-03-12 14:10:26 -04:00
"import matplotlib.pyplot as plt\n",
2025-03-13 13:47:02 -04:00
"\n",
2025-03-15 16:20:08 -04:00
"nest_asyncio.apply() # Allow nested event loops for Jupyter notebooks\n",
2025-03-15 15:01:32 -04:00
"# Simulate data\n",
"N = 100\n",
"alpha = 2.3\n",
"beta = 4.0\n",
"sigma = 2.0\n",
"\n",
"x = np.random.normal(size=N)\n",
"y = alpha + beta * x + sigma * np.random.normal(size=N)\n",
2025-03-15 16:20:08 -04:00
"sns.histplot(y, color=\"red\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"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% (1500/6000)\n",
"Sampling: 50% (3000/6000)\n",
"Sampling: 75% (4500/6000)\n",
"Sampling: 100% (6000/6000)\n",
"Sampling: 100% (6000/6000), done.\n",
"Messages received during sampling:\n",
" Gradient evaluation took 4.7e-05 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 0.47 seconds.\n",
" Adjust your expectations accordingly!\n",
" Gradient evaluation took 4.9e-05 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.\n",
" Adjust your expectations accordingly!\n",
" Gradient evaluation took 1e-05 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.\n",
" Adjust your expectations accordingly!\n",
" Gradient evaluation took 1.5e-05 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.\n",
" Adjust your expectations accordingly!\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \\\n",
"alpha 2.701 0.197 2.344 3.082 0.003 0.003 3785.0 2755.0 \n",
"beta 3.838 0.176 3.488 4.142 0.003 0.003 3854.0 3215.0 \n",
"sigma 1.987 0.143 1.734 2.265 0.002 0.002 3842.0 2617.0 \n",
"\n",
" r_hat \n",
"alpha 1.0 \n",
"beta 1.0 \n",
"sigma 1.0 \n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAJOCAYAAABm7rQwAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXecJGWd/99V3T1xd2cDG0hLDkoSUYJyCiKKgOJ54olizulU9OQQ9aenggmVA0/UMwMqoKIiQTISdhd22V02x8l5pnN3xef5/fFU6p6ws3EIz9sX7sx01VNPPfXU0/X91DcYUkqJRqPRaDQajUaj0Wg0Go1Gsw8xp7sDGo1Go9FoNBqNRqPRaDSaFx5alNJoNBqNRqPRaDQajUaj0exztCil0Wg0Go1Go9FoNBqNRqPZ52hRSqPRaDQajUaj0Wg0Go1Gs8/RopRGo9FoNBqNRqPRaDQajWafo0UpjUaj0Wg0Go1Go9FoNBrNPkeLUhqNRqPRaDQajUaj0Wg0mn2OFqU0Go1Go9FoNBqNRqPRaDT7HC1KaTQajUaj0Wg0Go1Go9Fo9jlalNJoNHuN9773vZx11lk7vV97ezuGYfCrX/1qj/dJo9FoNBqN5vmMfv7SaDTPJbQopdFoNBqNRqPRaDQajUaj2edoUUqj0Wg0Go1Go9FoNBqNRrPP0aKURqPZadasWcMll1zCIYccQnNzM0cccQQf//jHyWazk+730EMPYRgGt99+O5deeiltbW3Mnj2bD3/4w1QqlTHbe57HF7/4RRYuXMjcuXN5+9vfzujoaM023/nOdzjttNOYM2cOs2fP5hWveAV33nnnHj1fjUaj0Wg0mulGP39pNJrnI+np7oBGo3nu0dXVxVFHHcXb3/525syZw5YtW7jqqqt4+umneeKJJ3a4/yc/+UkuvPBCbrnlFlatWsWXv/xlqtUqv/3tb2u2+8Y3vsGZZ57Jr371KwYHB/nc5z7Hf/zHf3DjjTdG23R0dPDRj36UQw45BMdx+Otf/8oFF1zAXXfdxXnnnbfHz12j0Wg0Go1mOtDPXxqN5vmIIaWU090JjUbz3MbzPJ544gle9apXsWLFCk4++WRAJdpsb2/noYceAtSburPPPps3velN/OUvf4n2/973vsfll1/OunXrOOaYY2hvb+ewww7jNa95Dffff3/NdldeeSWWZWEYxph+CCEQQvCGN7yBlpaWmmNoNBqNRqPRPJ/Qz18ajeb5gA7f02g0O43jOFx11VUce+yxNDc3k8lkeNWrXgXAxo0bd7j/W9/61prf3/a2tyGEYNmyZTV/v+CCC2p+P+GEE3Ach4GBgehvy5cv58ILL2ThwoWk02kymQz33XfflPqh0Wg0Go1G81xBP39pNJrnI1qU0mg0O80VV1zB1VdfzYc+9CHuvPNOnnzySf70pz8BYFnWDvdfsGBBze8LFy4EoLe3t+bvc+fOrfm9sbGx5hhdXV2cc845eJ7Hj3/8Yx5//HGefPJJzjvvvCn1Q6PRaDQajea5gn7+0mg0z0d0TimNRrPT/P73v+cLX/gCn/vc56K/5fP5Ke8/ODhY83v45u2AAw7YqX7cfffdUR6DhoaG6O+lUmmn2tFoNBqNRqN5tqOfvzQazfMR7Sml0Wh2mkqlQiaTqfnbL37xiynvf9ttt9X8fsstt2CaJqeddtpO9yOVSmGa8VK2YcOGKSX71Gg0Go1Go3kuoZ+/NBrN8xHtKaXRaHaa8847j+9973ssXLiQAw44gFtuuYWlS5dOef8VK1bw0Y9+lLe85S2sXLmSr3zlK7zjHe/g6KOP3ql+vPa1r+Xzn/8873rXu/jABz5AZ2cn/+///T8WL16MEGJnT0uj0Wg0Go3mWYt+/tJoNM9HtKeURqPZaa677jpe+9rXctlll/H2t78dy7L43e9+t1P7F4tFLr74Yq666ire/e53c8MNN+x0P4477jhuuukmnn76aS688EKuueYavvvd70ZJPzUajUaj0WieL+jnL41G83zEkFLK6e6ERqN5YRCWJH7wwQc566yzprs7Go1Go9FoNM979POXRqN5NqM9pTQajUaj0Wg0Go1Go9FoNPscLUppNBqNRqPRaDQajUaj0Wj2OTp8T6PRaDQajUaj0Wg0Go1Gs8/RnlIajUaj0Wg0Go1Go9FoNJp9jhalNBqNRqPRaDQajUaj0Wg0+xwtSmk0Go1Go9FoNBqNRqPRaPY56alsJISgt7eXmTNnYhjG3u6TRqPRaDQazR5DSkmxWOSAAw7ANJ877+P085dGo9FoNJrnKlN9/pqSKNXb28vBBx+8xzqn0Wg0Go1Gs6/p6urioIMOmu5uTBn9/KXRaDQajea5zo6ev6YkSs2cOTNqbNasWXumZxqNRqPRaDT7gEKhwMEHHxw9zzxX0M9fGo1Go9FonqtM9flrSqJU6DI+a9Ys/VCk0Wg0Go3mOclzLQROP39pNBqNRqN5rrOj56/nTmIFjUaj0Wg0Go1Go9FoNBrN8wYtSmk0Go1Go9FoNBqNRqPRaPY5Uwrf02g0ml3Gd6HYD24FPBt8R/29cab6r3UBpPRSpNFoNBqN5rlP1fFpbkhNdzc0Go3mOYO2BDUazZ6hNAh9q6BvJQxugHw35Lug2AdSTLyfmYG5h8P8Y2Dx6XDIK2HRCWDqBzqNRqPRaDTPHUq2x/3rBzjugDaOXDBjuruj0WieBfi+j+u6092NvUImkyGV2n2bTYtSGo1m1ygNwfaHof1R9d/I5uADA+YcCrMXwxFnQ9vB0HYQNMyAdCOkGkBKcIpg5SHXCcObYGANrP+raqJ1AZz07/CSS2HBsdN1hhqNRqPRaDRTpuJ4AAyXbC1KaTQvcKSU9Pf3k8vlprsre5XZs2ezaNGi3Somo0UpjUYzdUa3w5o/wqa7ofspQEJTGxxyJrzsfXDgKbDweGjcxQexYr8SuNb9BZbcAI9fB4e9Gl71n3DomfAcq5yl0Wg0Go1Go9FoXniEgtSCBQtoaWl5zlUA3hFSSiqVCoODgwDsv//+u9yWFqU0Gs3k2CUlEq28GToeBQwlPr3mSjjy3D0bajdzEZzwVvVfeQRW3gRL/hd+fSEcfDq84dtwwEv2zLE0Go1Go9FoNBqNZg/j+34kSM2bN2+6u7PXaG5uBmBwcJAFCxbsciifFqU0Gs345Lth6Q2w/NdgF2DeUfDar8KJ/w6zDtj7x2+dB6/8DzjtI0qceujb8LOz4dSPKEGscebe74NGo9G8gPGLRapPP03rGWdgZDLT3R2N5jmDlNPdA41GM52EOaRaWlqmuSe7jytcKm6FmQ0zMQ1zzOfhObquq0UpjUazh8j3wMPfUp5RUsLxb1FC0EEvm57wuXQjvOz9cPxb4YFvKKFsw9/h4l+qPmk0Go1mr2Bv2oRfKOKNjpJZuHC6u6PRaDQazXOK50PInu3ZAAgpxhWl9sQ5alFKo9EorDw8/B1Y9jP1+6kfhjM+oZKUPxtomgXnfwdOfBv88QPwi9fDOf8PzvgkmGMXSI1Go9FoNBqNRqPR7D4Ge09g05acRvNCR0p45ja4/uUqf9OJb4P/WAHnXf3sEaSSHPQy+Mgj8KI3wr1fhtveC251unul0Wg0Go1Go9FoNHuc9vZ2DMNg5cqVu9XOWWedxWc+85md2key9+ORtaeURvNCpjgAt38Mtt4PB74M3nkr7H/SdPdqxzS1wVt/CQedCvd8UeW/uuT3MGPBdPdMo9FonndUnnyK1lecQXru3Onuikaj0Wg0ml3kT3/6E5mdzBEp90GSPO0ppdG8UNlyH9zwSuh8Ai64Bj5w73NDkAoxDDjj43DJ72BwA/zsHBhYN9290mg0muclXn//dHdBo9FMIyPVEdrz7dPdDY1GsxvMnTuXmTOffcWitCil0bzQ8F34x5fgxn+DGQvhww/Byz/4rMnLJISkY6TM3Wv6+MnDW/nqX9fysRuXc+n/LeWtP36cd/18KZ+8eQVfuv0Zvv+PjfypfAKbLrgFIVz4xXnQuXS6T0Gj0Wg0muc1d267ky3ZLdPdDc0+ZGnfUtaN6Jd/Gs3eRAjBd77zHY488kgaGxt
"text/plain": [
"<Figure size 1200x600 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABsMAAAHxCAYAAAA4D8L7AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAA6pxJREFUeJzs3Xd8pGW5//HvM5Pee7KbZJNszfbeK+wCC0sHEQWkqShgAUWxYTt6RP3B8Sh6VBBQQGnS67K9sL33nk3ZTbLpPZmZ5/fHzESWbenPlM/79eJ1jjPJM9+sAnfu676vyzBN0xQAAAAAAAAAAAAQgGxWBwAAAAAAAAAAAAB6C8UwAAAAAAAAAAAABCyKYQAAAAAAAAAAAAhYFMMAAAAAAAAAAAAQsCiGAQAAAAAAAAAAIGBRDAMAAAAAAAAAAEDAohgGAAAAAAAAAACAgEUxDAAAAAAAAAAAAAGLYhgAAAAAAAAAAAACFsUwIMAsX75chmFo3rx5PfZMwzBkGEaPPQ8AAMAfsAYCAADofb2xlwUAn0YxDAAAAAAsduzYMRmGodzcXKujAAAAAEDACbE6AAAAAAAAAAAgOE2ZMkV79+5VVFSU1VEABDCKYQAAAAAAAAAAS0RFRSk/P9/qGAACHG0SAR+3YcMGfec739GUKVOUkZGhsLAwpaen66qrrtJHH33U4ed8svWOw+HQr3/9a40cOVKRkZFKSUnRTTfdpH379l3wOa+++qpmzZqluLg4RUdHa+bMmXr33XfP+rV79uzRj3/8Y82cOVOZmZkKCwtTcnKyFixYoJdeeqnD2QEAAKz217/+VRMnTlR0dLQSEhJ0xRVXaN26def8eofDoSeffFLz5s1TUlKSwsPDlZeXp69+9asqLCw87WvvuOMO5eXlSZIKCgraZ5V9emZZXV2d/vrXv+r666/XkCFDFB0drejoaI0ePVo/+MEPVF1d3Ss/OwAAQFccPHhQd911l/Ly8hQeHq6YmBjl5ORo0aJFevrpp9u/7kIzw1avXq2FCxcqISFBMTExmjx5sv7+979LOveM10++/txzz2nKlCmKiYlRamqqPve5z+n48eOSJNM09Yc//EHjxo1TdHS0UlJSdMcdd6isrOyMZ7a1tem5557TLbfcovz8fMXFxSkyMlLDhg3T17/+dZWUlHT3jwxALzJM0zStDgHg3BYsWKBly5Zp5MiRys7OVnR0tA4fPqwtW7ZIkv7nf/5H3/jGN9q/fvny5brooos0d+5cLV++vP31Y8eOKS8vTzk5OZo4caLeeustzZ07VykpKdqwYYOOHDmimJgYffjhh5o+ffppGbyLh0ceeUQ///nPNWPGDGVlZWnfvn3avn27DMPQq6++quuuu+607/viF7+op556Svn5+crJyVFCQoKOHz+u9evXy+Vy6YEHHtBjjz3WS39yAAAA3eNdAz3wwAP6n//5H82cOVPZ2dnauXOndu3apZCQEL300ktnrIHq6up09dVXa/ny5YqJidHEiROVmpqqnTt3av/+/UpOTtbixYs1fvx4SdKTTz6p999/X6+++qqio6N14403nva8Z555RpJ7I2j27NlKTU3VsGHDlJmZqaqqKm3evFkVFRUaPHiw1q1bp+Tk5N7/wwEAADiPXbt2aebMmaqtrdWwYcM0cuRI2e12FRUVaefOnRo0aJC2bdsm6dx7WZL0r3/9S7fccotcLpdGjx6tUaNGqbi4WKtXr9ZDDz2kRx99VJK7qPVJ3nXcww8/rN/+9reaM2eOkpKStGHDBh0/flzZ2dnavn27vvKVr+jNN9/UvHnzFBkZqTVr1qisrExjxozRxo0bFRYW1v7MoqIiZWdnKz4+XsOHD1d2drYaGhq0bds2lZSUKDU1VWvXrtXgwYN77w8WQNeZAHzau+++a5aUlJzx+tq1a824uDgzNDTULCoqan992bJlpiRz7ty5p3390aNHTUmmJDMlJcXcvn17+3sOh8P82te+Zkoyc3JyzObm5tO+1/t9CQkJ5rp1605778c//rEpyRw6dOgZGZcvX24ePnz4jNf37dtnZmVlmZLM9evXd+jPAQAAoK9510CRkZHmkiVLTnvv17/+tSnJjI+PN0tLS0977/Of/7wpybzyyivPeO/xxx83JZlDhgwxHQ5H++vetVpOTs458xQWFpofffSR6XQ6T3u9oaHB/MIXvmBKMu+9994u/rQAAAA958477zQlmf/1X/91xnuNjY3mihUr2v/zufayiouLzZiYGFOS+bvf/e6091asWGFGR0e3r9c+zft6cnKyuW3bttM+e9asWaYkc/To0eagQYPMY8eOtb9fXl5uDh482JRkPvfcc6c9s7a21nzjjTfMlpaW015vbW01v/e975mSzCuuuOLCfzgALEGbRMDHXX755erXr98Zr0+fPl333Xef2tra9MYbb3TqmT/84Q81ZsyY9v9st9v1m9/8RpmZmSooKNCrr7561u/72c9+pqlTp5722ve+9z3Fx8frwIEDZ7T8mTt3rgYOHHjGc4YNG6Yf/ehHkqRXXnmlU9kBAAD62j333KOLL774tNceeughTZo0STU1NXryySfbX9+7d6/++c9/qn///nrhhReUlpZ22vd985vf1BVXXKGDBw/qvffe61SOrKwszZ8/Xzbb6b/GRUVF6U9/+pNCQkL08ssvd/KnAwAA6HmlpaWSpCuuuOKM9yIjIzVnzpwLPuOpp55SfX29pk+frq9//eunvTdnzhx99atfveAzfvazn2ns2LGnffaDDz4oSdq5c6f+93//Vzk5Oe3vp6SktD93yZIlpz0rNjZWV1999Wm3xSQpNDRUv/zlL9W/f3+9//77qquru2AuAH0vxOoAAC6soqJC77zzjnbt2qWqqiq1tbVJcvdelqT9+/d36nm33377Ga+Fh4frs5/9rB577DEtX75cn//858/4mquuuuqs3zdw4EBt3bpVxcXFys7OPu39+vp6vffee9q6datOnTql1tZWSdKJEye6lB0AAKCvnW3tJElf+MIXtGnTJi1fvlzf//73JUnvvvuuTNPU5ZdfrtjY2LN+37x58/Tuu+9q7dq1uvLKKzudZ+3atVq1apWOHz+uxsbG9rZAYWFhKi8vV1VVlRITEzv9XAAAgJ4yZcoUvfvuu/rqV7+qn/70p5o7d64iIiI69YwVK1ZIkm655Zazvn/LLbfot7/97XmfcbZi3JAhQyRJISEhuvTSS8/5/rlmgG3fvl1LlizR0aNH1dDQIJfLJck9M9blcunQoUPt7bAB+A6KYYCP++tf/6oHHnhADQ0N5/ya2traDj8vISFBCQkJZ33PO7i9qKjorO8PGDDgrK/HxcVJkpqbm097/a233tKdd96pioqKc+bpTHYAAAAreNdI53r9k2unI0eOSHKfZH7qqafO+9zy8vJO5SgrK9MNN9yg1atXn/framtrKYYBAABLPfTQQ1q9erU++ugjLVy4UKGhoRo7dqzmzJmjm2++WZMnT77gM7xrrNzc3LO+f67XP+lse1kxMTGSpH79+ikk5Mztce+Bpk/vczU0NOi2227Ta6+9dt7PZK8L8E0UwwAftnnzZt1zzz2y2+169NFHddVVV2nAgAGKioqSYRj6y1/+onvuueeMIaHdda7nfbolz/kUFxfrs5/9rJqamvSd73xHt9xyi3JzcxUTEyObzaYPP/xQl112WY9nBwAA6GufXM94TwaPGzfutJY8Z/Pp9tMX8sUvflGrV6/W9OnT9dOf/lRjx45VYmKiQkNDJUn9+/fXiRMnWF8BAADLRUVFafHixdq4caPef/99rV27VmvXrtWmTZv02GOP6d5779UTTzzRoWcZhtGp1z/pfHtZndnnktyjQl577TXl5+frV7/6lSZPnqyUlJT2tokzZszQxx9/zFoM8FEUwwAf9vLLL8s0TX3ta1/Td77znTPe97ZJ7Izq6mpVV1ef9XbYsWPHJLnnUXTXW2+9paamJl133XV69NFHz3i/K9kBAACscPToUY0bN+6M18+2dvK2jJ45c6b+8Ic/9FiGhoYGvfvuu7LZbHr33XfPWMs1NDT
"text/plain": [
"<Figure size 2208x552 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"\n",
2025-03-15 15:01:32 -04:00
"\n",
"stan_code = \"\"\"\n",
"data {\n",
" int<lower=1> N; // Number of observations\n",
" vector[N] x; // Covariate\n",
" vector[N] y; // Outcome\n",
"}\n",
"\n",
"parameters {\n",
" real alpha; // Intercept\n",
" real beta; // Slope\n",
" real<lower=0> sigma2; // Variance (sigma squared)\n",
"}\n",
"\n",
"transformed parameters {\n",
" real<lower=0> sigma; // Standard deviation\n",
" sigma = sqrt(sigma2);\n",
"}\n",
"\n",
"model {\n",
" // Priors\n",
" alpha ~ normal(0, 10);\n",
" beta ~ normal(0, 10);\n",
" sigma2 ~ inv_gamma(1, 1);\n",
"\n",
" // Likelihood\n",
" y ~ normal(alpha + beta * x, sigma);\n",
"}\n",
"\"\"\"\n",
"# Prepare data dictionary\n",
"data_dict = {\n",
" 'N': N,\n",
" 'x': x,\n",
" 'y': y\n",
"}\n",
"\n",
"# Build and fit the model\n",
"model = stan.build(stan_code, data=data_dict)\n",
"fit = model.sample(num_chains=4, num_samples=1000, num_warmup=500)\n",
"summary = az.summary(fit, var_names=['alpha', 'beta', 'sigma'])\n",
"print(summary)\n",
"\n",
"# Trace plots for convergence diagnostics\n",
"az.plot_trace(fit, var_names=['alpha', 'beta', 'sigma'], compact=False, legend=True)\n",
"plt.tight_layout()\n",
"plt.show()\n",
"\n",
"# Plot posterior distributions\n",
"az.plot_posterior(fit, var_names=['alpha', 'beta', 'sigma'])\n",
2025-03-12 14:10:26 -04:00
"plt.show()"
]
2025-03-15 15:01:32 -04:00
},
{
"cell_type": "code",
2025-03-15 16:20:08 -04:00
"execution_count": 5,
2025-03-15 15:01:32 -04:00
"metadata": {},
2025-03-15 16:20:08 -04:00
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Building...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"Building: found in cache, done.Sampling: 0%"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"Sampling: 25% (1500/6000)\n",
"Sampling: 50% (3000/6000)\n",
"Sampling: 75% (4500/6000)\n",
"Sampling: 100% (6000/6000)\n",
"Sampling: 100% (6000/6000), done.\n",
"Messages received during sampling:\n",
" Gradient evaluation took 0.000227 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 2.27 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_tzmle_c1/model_hxisyy5r.stan', line 26, column 4 to column 40)\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 0.00388 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 38.8 seconds.\n",
" Adjust your expectations accordingly!\n",
" Gradient evaluation took 0.000115 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 1.15 seconds.\n",
" Adjust your expectations accordingly!\n",
" Gradient evaluation took 4.5e-05 seconds\n",
" 1000 transitions using 10 leapfrog steps per transition would take 0.45 seconds.\n",
" Adjust your expectations accordingly!\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \\\n",
"alpha 2.220 0.064 2.102 2.344 0.001 0.001 3953.0 3010.0 \n",
"beta 3.980 0.065 3.851 4.097 0.001 0.001 3775.0 2697.0 \n",
"sigma 1.993 0.044 1.912 2.076 0.001 0.001 3210.0 2842.0 \n",
"\n",
" r_hat \n",
"alpha 1.0 \n",
"beta 1.0 \n",
"sigma 1.0 \n",
"Uncertainty (N=100):\n",
" mean sd hdi_3% hdi_97%\n",
"alpha 2.701 0.197 2.344 3.082\n",
"beta 3.838 0.176 3.488 4.142\n",
"sigma 1.987 0.143 1.734 2.265\n",
"\n",
"Uncertainty (N=1000):\n",
" mean sd hdi_3% hdi_97%\n",
"alpha 2.220 0.064 2.102 2.344\n",
"beta 3.980 0.065 3.851 4.097\n",
"sigma 1.993 0.044 1.912 2.076\n"
]
}
],
2025-03-15 15:01:32 -04:00
"source": [
"# Simulate larger dataset\n",
"N_large = 1000\n",
"x_large = np.random.normal(size=N_large)\n",
"y_large = alpha + beta * x_large + sigma * np.random.normal(size=N_large)\n",
"\n",
"# Prepare data dictionary for larger dataset\n",
"data_dict_large = {\n",
" 'N': N_large,\n",
" 'x': x_large,\n",
" 'y': y_large\n",
"}\n",
"\n",
2025-03-15 16:20:08 -04:00
"# Build and fit the model with larger dataset\n",
"model_large = stan.build(stan_code, data=data_dict_large)\n",
"fit_large = model_large.sample(num_chains=4, num_samples=1000, num_warmup=500)\n",
2025-03-15 15:01:32 -04:00
"\n",
"# Summarize results\n",
"summary_large = az.summary(fit_large, var_names=['alpha', 'beta', 'sigma'])\n",
"print(summary_large)\n",
"\n",
"# Compare uncertainty\n",
"print(\"Uncertainty (N=100):\")\n",
"print(summary[['mean', 'sd', 'hdi_3%', 'hdi_97%']])\n",
"print(\"\\nUncertainty (N=1000):\")\n",
"print(summary_large[['mean', 'sd', 'hdi_3%', 'hdi_97%']])"
]
2025-03-12 14:10:26 -04:00
}
],
"metadata": {
"kernelspec": {
2025-03-13 13:47:02 -04:00
"display_name": "venv",
2025-03-12 14:10:26 -04:00
"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",
2025-03-13 13:47:02 -04:00
"version": "3.12.3"
2025-03-12 14:10:26 -04:00
}
},
"nbformat": 4,
"nbformat_minor": 2
}
