2025-03-12 14:10:26 -04:00
{
"cells": [
2025-03-22 14:57:50 -04:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 1\n",
"- 2 False: The variance of a Wiener process with scale coefficient sigma = 1 and time t is t2\n",
"- 3 False: \n",
"- 4 False: \n",
"- 5 False: \n",
"- 6 False: \n",
"- 8 False: "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 3\n",
"Show that the following identity holds for any given prior and posterior:\n",
" Var [θ] = E[Var [θ |y]] + Var [E[θ |y]] \n",
"\n",
"Clarification of terms:\n",
"1. Var [θ] –
"2. E[Var [θ |y]] –
"3. Var [E[θ |y]] –
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem 4\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.stats import norm\n",
"\n",
"import matplotlib.pyplot as plt\n",
"from scipy.stats import norm\n",
"import pandas as pd\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sample Mean (ybar): 4.84 hours\n",
"Posterior Mean (µ_n): 4.85\n",
"Posterior Standard Deviation (σ
]
}
],
"source": [
"np.random.seed(42)\n",
"n = 100 #sample size\n",
"mu_true = 5 #avg sleep time\n",
"sigma = 1.5 #standard deviation\n",
"data = np.random.normal(mu_true, sigma, n)\n",
"ybar = np.mean(data) \n",
"print(f\"Sample Mean (ybar): {ybar:.2f} hours\")\n",
"\n",
"# Prior parameters (initial belief about sleep time)\n",
"mu0 = 5 \n",
"sigma0 = 1.0\n",
"\n",
"# Posterior parameters\n",
"sigma_n = np.sqrt(1 / (1 / sigma0**2 + n / sigma**2))\n",
"mu_n = (mu0 / sigma0**2 + n * ybar / sigma**2) / (1 / sigma0**2 + n / sigma**2)\n",
"\n",
"print(f\"Posterior Mean (µ_n): {mu_n:.2f}\")\n",
"print(f\"Posterior Standard Deviation (σ
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAIjCAYAAAA0vUuxAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjZ1JREFUeJzs3Xl8VPW9//HXzGQySSbbTHaSQELYFBAQBVkUXKlLldrWpbdFbPV6b9XqtXahv6q1y+Vq69Jaq7WtWrfaaq1W64Z7BdxYBQGBkA2yzmSdLJPMnN8fkwwEErKQ5GR5Px+Pecic+c45n0niJO/5nvP5WgzDMBAREREREZFuWc0uQEREREREZLhTcBIREREREemBgpOIiIiIiEgPFJxERERERER6oOAkIiIiIiLSAwUnERERERGRHig4iYiIiIiI9EDBSUREREREpAcKTiIiIiIiIj1QcBIRGUJLly5l6dKl4fsFBQVYLBYeffRR02oaKXJycli5cqXZZQyan/zkJ1gsFrPLEBGRbig4iYgcxd69e7nmmmuYOHEiUVFRxMfHs2jRIn7961/T1NRkdnlDpqCggCuvvJK8vDyioqJIT0/ntNNO47bbbus07ne/+92IDoEvv/wyP/nJTwZ0nw0NDdx2223MmDEDp9NJUlISs2fP5oYbbuDAgQMDeqzBtnLlSmJjY7t93GKxcN111w1hRSIiQyfC7AJERIarf/3rX3z1q1/F4XCwYsUKZsyYgd/v5/333+d73/se27dv56GHHjK7zEG3Z88eTj75ZKKjo/nmN79JTk4OpaWlbNy4kTvuuIPbb789PPZ3v/sdycnJI3Zm6OWXX+b+++8fsPDU2trKaaedxs6dO7niiiu4/vrraWhoYPv27Tz11FN86UtfYty4cQNyLBERGVwKTiIiXdi3bx+XXXYZEyZM4K233iIjIyP82LXXXsuePXv417/+ZWKFQ+eee+6hoaGBzZs3M2HChE6PVVRUmFTVyPD888+zadMmnnzySb72ta91eqy5uRm/329SZaOLz+fD6XSaXYaIjHI6VU9EpAt33nknDQ0N/OlPf+oUmjpMmjSJG264IXy/ra2Nn/3sZ+Tl5eFwOMjJyeFHP/oRLS0t/Tr+zp07+cpXvoLb7SYqKoqTTjqJf/7zn0eM27p1K0uWLCE6OpqsrCx+/vOf88gjj2CxWCgoKOg09pVXXuHUU0/F6XQSFxfH+eefz/bt23usZe/evWRlZR0RmgBSU1PD/87JyWH79u28++67WCwWLBZL+Hqu7q7fefTRR4+o1TAMfv7zn5OVlUVMTAynn356t3XW1NRw4403kp2djcPhYNKkSdxxxx0Eg8HwmI7ryH71q1/x0EMPhb9HJ598Mh9//HF43MqVK7n//vsBwvUfWvPTTz/N3LlziYuLIz4+npkzZ/LrX/+6x68dwKJFi454rOPUz5488cQTzJ07l+joaNxuN5dddhnFxcVHjPvwww/5whe+QEJCAjExMSxZsoS1a9d2GtPxfdi5cyeXXHIJ8fHxJCUlccMNN9Dc3NxjLf1RUVHBt771LdLS0oiKimLWrFn8+c9/7jTmnXfewWKx8M4773Ta3tU1gB2nC+7du5fzzjuPuLg4/uM//gOA3bt38+Uvf5n09HSioqLIysrisssuo7a2dlBem4iMLZpxEhHpwosvvsjEiRNZuHBhr8ZfddVV/PnPf+YrX/kK3/3ud/nwww9ZvXo1O3bs4B//+Eefjr19+3YWLVpEZmYmP/zhD3E6nfztb39j+fLl/P3vf+dLX/oSAPv37+f000/HYrGwatUqnE4nf/zjH3E4HEfs8/HHH+eKK65g2bJl3HHHHTQ2NvLAAw+wePFiNm3aRE5OTrf1TJgwgTfeeIO33nqLM844o9tx9957L9dffz2xsbH8v//3/wBIS0vr02sHuPXWW/n5z3/Oeeedx3nnncfGjRs555xzjpidaWxsZMmSJezfv59rrrmG8ePHs27dOlatWkVpaSn33ntvp/FPPfUU9fX1XHPNNVgsFu68804uvvhi8vPzsdvtXHPNNRw4cIA1a9bw+OOPd3rumjVruPzyyznzzDO54447ANixYwdr167tFKAP1xE2H3vsMX784x/3ufnDL37xC2655RYuueQSrrrqKiorK7nvvvs47bTT2LRpE4mJiQC89dZbnHvuucydO5fbbrsNq9XKI488whlnnMG///1v5s2b12m/l1xyCTk5OaxevZoPPviA3/zmN1RXV/PYY4/1qq6qqqpejWtqamLp0qXs2bOH6667jtzcXJ555hlWrlxJTU3NUb92R9PW1sayZctYvHgxv/rVr4iJicHv97Ns2TJaWlq4/vrrSU9PZ//+/bz00kvU1NSQkJDQr2OJiIQZIiLSSW1trQEYF110Ua/Gb9682QCMq666qtP2m2++2QCMt956K7xtyZIlxpIlS8L39+3bZwDGI488Et525plnGjNnzjSam5vD24LBoLFw4UJj8uTJ4W3XX3+9YbFYjE2bNoW3eTwew+12G4Cxb98+wzAMo76+3khMTDSuvvrqTvWVlZUZCQkJR2w/3LZt24zo6GgDMGbPnm3ccMMNxvPPP2/4fL4jxk6fPr3T6+tw2223GV39ynnkkUc61VpRUWFERkYa559/vhEMBsPjfvSjHxmAccUVV4S3/exnPzOcTqfx+eefd9rnD3/4Q8NmsxlFRUWGYRz8GiclJRlerzc87oUXXjAA48UXXwxvu/baa7us84YbbjDi4+ONtra2rr9I3WhsbDSmTp1qAMaECROMlStXGn/605+M8vLyI8Ye/jUqKCgwbDab8Ytf/KLTuE8//dSIiIgIbw8Gg8bkyZONZcuWdfqaNTY2Grm5ucbZZ599xDEuvPDCTvv89re/bQDGli1bjvp6rrjiCgM46u3aa68Nj7/33nsNwHjiiSfC2/x+v7FgwQIjNjbWqKurMwzDMN5++20DMN5+++1Ox+vq/4+OGn74wx92Grtp0yYDMJ555pmjvgYRkf7SqXoiIoepq6sDIC4urlfjX375ZQBuuummTtu/+93vAvTpWiiv18tbb73FJZdcQn19PVVVVVRVVeHxeFi2bBm7d+9m//79ALz66qssWLCA2bNnh5/vdrvDpy11WLNmDTU1NVx++eXh/VVVVWGz2Zg/fz5vv/32UWuaPn06mzdv5utf/zoFBQX8+te/Zvny5aSlpfGHP/yh16+tN9544w38fj/XX399p9mZG2+88YixzzzzDKeeeioul6vT6zrrrLMIBAK89957ncZfeumluFyu8P1TTz0VgPz8/B7rSkxMxOfzsWbNmj69nujoaD788EO+973vAaFTE7/1rW+RkZHB9ddff9RTOZ977jmCwSCXXHJJp9eXnp7O5MmTw9+3zZs3s3v3br72ta/h8XjC43w+H2eeeSbvvfdep1MXIXSd3qGuv/564ODP8tFERUWxZs2aLm+He/nll0lPT+fyyy8Pb7Pb7XznO9+hoaGBd999t8fjdee///u/O93vmFF67bXXaGxs7Pd+RUS6o1P1REQO03HdSX19fa/GFxYWYrVamTRpUqft6enpJCYmUlhY2Otj79mzB8MwuOWWW7jlllu6HFNRUUFmZiaFhYUsWLDgiMcPr2P37t0A3Z5m15vrbKZMmcLjjz9OIBDgs88+46WXXuLOO+/kP//zP8nNzeWss87qcR+90fG1mjx5cqftKSkpnUIPhF7X1q1bSUlJ6XJfhzeuGD9+fKf7Hfurrq7usa5vf/vb/O1vf+Pcc88lMzOTc845h0suuYQvfOELPT43ISGBO++8kzvvvJPCwkLefPNNfvWrX/Hb3/6WhIQEfv7zn3f5vN27d2MYxhFfiw52uz08DuCKK67otoba2tpOX7/D95mXl4fVaj3iuriu2Gy2Xn+/CwsLmTx5MlZr589pjzvuuPDj/REREUFWVlanbbm5udx0003cfffdPPnkk5x66qlceOGFfP3rX9dpeiIyIBScREQOEx8fz7hx49i2bVufnjcQi5d2zAzcfPPNLFu2rMsxhwej3u7z8ccfJz09/Yj
"text/plain": [
"<Figure size 1000x600 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Define x-axis range\n",
"x = np.linspace(3, 8, 1000)\n",
"\n",
"# Prior distribution\n",
"prior = norm.pdf(x, mu0, sigma0)\n",
"\n",
"# Likelihood distribution (centered at sample mean)\n",
"likelihood = norm.pdf(x, ybar, sigma / np.sqrt(n))\n",
"\n",
"# Posterior distribution\n",
"posterior = norm.pdf(x, mu_n, sigma_n)\n",
"\n",
"# Plot\n",
"plt.figure(figsize=(10, 6))\n",
"plt.plot(x, prior, label=\"Prior\", linestyle=\"--\")\n",
"plt.plot(x, likelihood, label=\"Likelihood\", linestyle=\"-.\")\n",
"plt.plot(x, posterior, label=\"Posterior\", linestyle=\"-\")\n",
"plt.title(\"College Students Sleep Hours\")\n",
"plt.xlabel(\"Sleep Time (hours)\")\n",
"plt.ylabel(\"Density\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Prior Data (Sample Mean) Posterior\n",
"0 Precision 1.0 33.333 34.333\n",
"1 SD 1.0 0.173 0.171\n",
"2 Mean 100.0 97.450 97.524\n"
]
}
],
"source": [
"# Create a DataFrame\n",
"df = pd.DataFrame({\n",
" \"\": [\"Precision\", \"SD\", \"Mean\"],\n",
" \"Prior\": [1 / sigma0**2, sigma0, mu0],\n",
" \"Data (Sample Mean)\": [n / sigma**2, sigma / np.sqrt(n), ybar],\n",
" \"Posterior\": [1 / sigma_n**2, sigma_n, mu_n]\n",
"})\n",
"\n",
"# Display the DataFrame\n",
"print(df.round(3))"
]
},
2025-03-12 14:10:26 -04:00
{
"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-22 14:57:50 -04:00
"execution_count": null,
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'>"
]
},
2025-03-22 14:57:50 -04:00
"execution_count": 6,
2025-03-15 16:20:08 -04:00
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
2025-03-22 14:57:50 -04:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKiJJREFUeJzt3X9QlPeBx/HP4o9Fc4g/QH4kKGoiolGMRKnWNFqpyPU8TVLPMKH+qDEdR5IYatKQ83cyR35c1KZSbTNVkvGsJjPG5BrPnmLUOqCpGC4hg4x6IhpZFBJYgQIKe3/03GYjoCDw7PJ9v2a+M3me5/s8+3mys+GTZ5/dtblcLpcAAAAM4md1AAAAgM5GAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMA4FCAAAGKe71QG8UWNjoy5duqSAgADZbDar4wAAgNvgcrl09epVhYeHy8+v5Ws8FKAmXLp0SREREVbHAAAAbXDhwgXdc889Lc6hADUhICBA0t/+Bfbp08fiNAAA4HY4nU5FRES4/463hALUhBtve/Xp04cCBACAj7md21e4CRoAABiHAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMA4FCAAAGIcCBAAAjGNpAUpPT9f48eMVEBCggQMHavbs2SosLPSYU1tbq6VLl2rAgAH6h3/4Bz322GMqLS1t8bgul0urVq1SWFiYevXqpfj4eJ0+fbojTwUAAPgQSwvQ4cOHtXTpUh07dkz79+/XtWvXNH36dFVXV7vnPPfcc/rP//xPvf/++zp8+LAuXbqkRx99tMXjvv7663rrrbe0ZcsWHT9+XHfddZcSEhJUW1vb0acEAAB8gM3lcrmsDnHDlStXNHDgQB0+fFg/+MEPVFlZqeDgYO3YsUM/+clPJEmnTp1SdHS0cnJy9L3vfe+mY7hcLoWHh+sXv/iFli9fLkmqrKxUSEiIMjMz9fjjj98yh9PpVGBgoCorK/kxVAAAfERr/n571a/BV1ZWSpL69+8vScrNzdW1a9cUHx/vnjNixAgNGjSo2QJ07tw5ORwOj30CAwMVFxennJycJgtQXV2d6urq3MtOp7PdzglojeLiYpWVlVkdo1WCgoI0aNAgq2MAQKt4TQFqbGzUsmXL9P3vf1/333+/JMnhcKhnz57q27evx9yQkBA5HI4mj3NjfUhIyG3vk56errVr197hGQB3pri4WNFRUarxsbdqe/v7q6CwkBIEwKd4TQFaunSp8vPzdfTo0U5/7LS0NKWmprqXnU6nIiIiOj0HzFZWVqaa2lptj45WdO/eVse5LQU1NUouKFBZWRkFCIBP8YoClJKSoj/+8Y86cuSI7rnnHvf60NBQ1dfXq6KiwuMqUGlpqUJDQ5s81o31paWlCgsL89hn7NixTe5jt9tlt9vv/ESAdhDdu7fGBQRYHQMAujRLPwXmcrmUkpKiDz74QAcPHtSQIUM8tsfGxqpHjx7KyspyryssLFRxcbEmTpzY5DGHDBmi0NBQj32cTqeOHz/e7D4AAMAslhagpUuXavv27dqxY4cCAgLkcDjkcDj017/+VdLfbl5etGiRUlNT9cknnyg3N1cLFy7UxIkTPW6AHjFihD744ANJks1m07Jly/TKK6/oo48+0hdffKF58+YpPDxcs2fPtuI0AQCAl7H0LbDNmzdLkqZMmeKxftu2bVqwYIEkacOGDfLz89Njjz2muro6JSQk6De/+Y3H/MLCQvcnyCTphRdeUHV1tZ566ilVVFRo8uTJ2rdvn/z9/Tv0fAAAgG/wqu8B8hZ8DxCscPLkScXGxio3NtZn7gE6efWqYnNzlZubq3HjxlkdB4DhWvP3m98CAwAAxqEAAQAA41CAAACAcShAAADAOBQgAABgHAoQAAAwDgUIAAAYhwIEAACMQwECAADGoQABAADjUIAAAIBxKEAAAMA4FCAAAGAcChAAADAOBQgAABiHAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMA4FCAAAGIcCBAAAjEMBAgAAxqEAAQAA41CAAACAcbpbHQAAOltxcbHKysqsjtEqQUFBGjRokNUxgC6DAgTAKMXFxYqOilJNba3VUVqlt7+/CgoLKUFAO6EAATBKWVmZamprtT06WtG9e1sd57YU1NQouaBAZWVlFCCgnVCAABgpundvjQsIsDoGAItwEzQAADAOBQgAABjH0gJ05MgRzZw5U+Hh4bLZbNqzZ4/HdpvN1uR44403mj3mmjVrbpo/YsSIDj4TAADgSywtQNXV1YqJiVFGRkaT20tKSjzG1q1bZbPZ9Nhjj7V43FGjRnnsd/To0Y6IDwAAfJSlN0EnJiYqMTGx2e2hoaEeyx9++KGmTp2qoUOHtnjc7t2737QvAADADT5zD1Bpaak+/vhjLVq06JZzT58+rfDwcA0dOlRPPPGEiouLW5xfV1cnp9PpMQAAQNflMwXonXfeUUBAgB599NEW58XFxSkzM1P79u3T5s2bde7cOT300EO6evVqs/ukp6crMDDQPSIiIto7PgAA8CI+U4C2bt2qJ554Qv7+/i3OS0xM1Jw5czRmzBglJCRo7969qqio0HvvvdfsPmlpaaqsrHSPCxcutHd8AADgRXziixD//Oc/q7CwULt27Wr1vn379tXw4cN15syZZufY7XbZ7fY7iQgAAHyIT1wB+v3vf6/Y2FjFxMS0et+qqiqdPXtWYWFhHZAMAAD4IksLUFVVlfLy8pSXlydJOnfunPLy8jxuWnY6nXr//ff15JNPNnmMadOmadOmTe7l5cuX6/DhwyoqKlJ2drYeeeQRdevWTUlJSR16LgAAwHdY+hbYiRMnNHXqVPdyamqqJGn+/PnKzMyUJO3cuVMul6vZAnP27FmVlZW5ly9evKikpCSVl5crODhYkydP1rFjxxQcHNxxJwIAAHyKpQVoypQpcrlcLc556qmn9NRTTzW7vaioyGN5586d7RENAAB0YT5xDxAAAEB78olPgQHwbgUFBVZHuG2+lBVAx6EAAWizkvp6+UlKTk62Okqr1dXXWx0BgIUoQADarOL6dTVKejsyUuMGDLA6zm3ZW16ulUVFun79utVRAFiIAgTgjkX16qVxAQFWx7gtBTU1VkcA4AW4CRoAABiHAgQAAIxDAQIAAMahAAEAAONQgAAAgHEoQAAAwDgUIAAAYBwKEAAAMA4FCAAAGIcCBAAAjEMBAgAAxqEAAQAA41CAAACAcShAAADAOBQgAABgHAoQAAAwDgUIAAAYhwIEAACMQwECAADGoQABAADjUIAAAIBxKEAAAMA4FCAAAGAcChAAADAOBQgAABiHAgQAAIxDAQIAAMahAAEAAONYWoCOHDmimTNnKjw8XDabTXv27PHYvmDBAtlsNo8xY8aMWx43IyNDkZGR8vf3V1xcnD799NMOOgMAAOCLLC1A1dXViomJUUZGRrNzZsyYoZKSEvf4wx/+0OIxd+3apdTUVK1evVonT55UTEyMEhISdPny5faODwAAfFR3Kx88MTFRiYmJLc6x2+0KDQ297WOuX79eixcv1sKFCyVJW7Zs0ccff6ytW7fqxRdfvKO8AACga/D6e4AOHTqkgQMHKioqSkuWLFF5eXmzc+vr65Wbm6v4+Hj3Oj8/P8XHxysnJ6fZ/erq6uR0Oj0GAADoury6AM2YMUPvvvuusrKy9Nprr+nw4cNKTExUQ0NDk/PLysrU0NCgkJAQj/UhISFyOBzNPk56eroCAwPdIyIiol3PAwAAeBdL3wK7lccff9z9z6NHj9aYMWM0bNgwHTp0SNOmTWu3x0lLS1Nqaqp72el0UoIAAOjCvPoK0HcNHTpUQUFBOnPmTJPbg4KC1K1bN5WWlnqsLy0tbfE+Irvdrj59+ngMAADQdflUAbp48aLKy8sVFhbW5PaePXsqNjZWWVlZ7nWNjY3KysrSxIkTOysmAADwcpYWoKqqKuXl5SkvL0+SdO7cOeXl5am4uFhVVVV6/vnndezYMRUVFSkrK0uzZs3Svffeq4SEBPc
2025-03-15 16:20:08 -04:00
"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-22 14:57:50 -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",
2025-03-22 14:57:50 -04:00
"execution_count": null,
2025-03-15 16:20:08 -04:00
"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": [
2025-03-15 15:01:32 -04:00
"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
}
