COGMOD-HWI/HW03/hw3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample Mean (ybar): 4.66 hours\n",
      "Posterior Mean (µ_n): 4.70\n",
      "Posterior Standard Deviation (σ_n): 0.21\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "n = 50 #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.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 (σ_n): {sigma_n:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAIjCAYAAADLH25TAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAATBhJREFUeJzt3XlcFuX+//H3DbLJ6gIiiohg7opbhuZWGq5pWprZwb1OkkuWJZ7cUiNL006p6am0XH62uhzTzH3LOm64Zy64VO4KCCoozO8Pv951CygoeMPwej4e9+PBzFwz85l7LHl7zXWNxTAMQwAAAABgIg72LgAAAAAAchtBBwAAAIDpEHQAAAAAmA5BBwAAAIDpEHQAAAAAmA5BBwAAAIDpEHQAAAAAmA5BBwAAAIDpEHQAAAAAmA5BBwAKoWbNmqlZs2bW5WPHjslisWj27Nl2q6mgKF++vHr27GnvMvLM6NGjZbFY7F0GANw3gg4AU5g9e7YsFou2bduW6fZmzZqpevXqD7iq3HPkyBG9+OKLqlChglxdXeXl5aVGjRrpgw8+0NWrV+1d3gNz7Ngx9erVSyEhIXJ1dZW/v7+aNGmiUaNG2bSbNm1agQ5ty5Yt0+jRo3P1mElJSRo1apSqV68ud3d3lShRQmFhYRo0aJD+/PPPXD0XAOQHRexdAADgzr7//ns988wzcnFxUWRkpKpXr67U1FRt2rRJQ4cO1b59+zRz5kx7l5nnDh8+rPr168vNzU29e/dW+fLlderUKe3YsUMTJkzQmDFjrG2nTZumkiVLFtiel2XLlmnq1Km5FnauX7+uJk2a6Ndff1WPHj00YMAAJSUlad++fZo/f76eeuopBQQE5Mq5ACC/IOgAQB5KTk6Wu7v7Pe8fFxenZ599VkFBQVqzZo1Kly5t3RYVFaXDhw/r+++/z41S873JkycrKSlJsbGxCgoKstl29uxZO1VVMCxatEg7d+7UvHnz9Nxzz9lsu3btmlJTU+1UGQDkHR5dA1Bo3bhxQ2PHjlVISIhcXFxUvnx5DR8+XCkpKTbtLBZLpv+yfvtYjVuPz61fv179+/eXn5+fypYtK0m6fPmyBg8erPLly8vFxUV+fn5q2bKlduzYccca3333XSUlJenTTz+1CTm3hIaGatCgQTm+puz69ddf9fTTT6t48eJydXVVvXr1tGTJkgztdu/eraZNm8rNzU1ly5bVuHHjNGvWLFksFh07dsym7fLly9W4cWO5u7vL09NTbdu21b59++5ay5EjR1S2bNkMIUeS/Pz8rD+XL19e+/bt0/r162WxWGSxWKzjkbIaf3Lr3v29VsMwNG7cOJUtW1ZFixZV8+bNs6wzPj5egwcPVmBgoFxcXBQaGqoJEyYoPT3d2ubWOKiJEydq5syZ1ntUv359bd261dquZ8+emjp1qiRZ6/97zQsWLFDdunXl6ekpLy8v1ahRQx988MFdvztJatSoUYZttx6FvJu5c+eqbt26cnNzU/HixfXss8/q5MmTGdr98ssvatWqlby9vVW0aFE1bdpUmzdvtmlz6z78+uuv6tKli7y8vFSiRAkNGjRI165du2stAJAd9OgAMJWEhASdP38+w/rr169nWNe3b199/vnnevrpp/Xqq6/ql19+UUxMjA4cOKCFCxfecw39+/eXr6+vRo4cqeTkZEnSP//5T33zzTd6+eWXVbVqVV24cEGbNm3SgQMHVKdOnSyP9d///lcVKlRQw4YNs3Xu3Lymffv2qVGjRipTpoyGDRsmd3d3ffXVV+rYsaO+/fZbPfXUU5KkP/74Q82bN5fFYlF0dLTc3d31ySefyMXFJcMx58yZox49eigiIkITJkzQlStXNH36dD366KPauXOnypcvn2U9QUFBWrVqldasWaPHHnssy3ZTpkzRgAED5OHhoX/961+SpFKlSuXo2iVp5MiRGjdunNq0aaM2bdpox44deuKJJzL0fly5ckVNmzbVH3/8oRdffFHlypXTTz/9pOjoaJ06dUpTpkyxaT9//nxdvnxZL774oiwWi95991116tRJR48elZOTk1588UX9+eefWrlypebMmWOz78qVK9WtWzc9/vjjmjBhgiTpwIED2rx5s03gvd2tcPjFF1/ozTffzPFkA+PHj9eIESPUpUsX9e3bV+fOndOHH36oJk2aaOfOnfLx8ZEkrVmzRq1bt1bdunU1atQoOTg4aNasWXrssce0ceNGPfzwwzbH7dKli8qXL6+YmBj9/PPP+ve//61Lly7piy++yFF9AJApAwBMYNasWYakO36qVatmbR8bG2tIMvr27WtznNdee82QZKxZs8a6TpIxatSoDOcMCgoyevTokaGGRx991Lhx44ZNW29vbyMqKipH15SQkGBIMjp06JCt9jm5pqZNmxpNmza1LsfFxRmSjFmzZlnXPf7440aNGjWMa9euWdelp6cbDRs2NCpWrGhdN2DAAMNisRg7d+60rrtw4YJRvHhxQ5IRFxdnGIZhXL582fDx8TH69etnU9/p06cNb2/vDOtvt3fvXsPNzc2QZISFhRmDBg0yFi1aZCQnJ2doW61aNZvru2XUqFFGZn/13bp3t2o9e/as4ezsbLRt29ZIT0+3ths+fLghyea+jx071nB3dzd+++03m2MOGzbMcHR0NE6cOGEYxl/fcYkSJYyLFy9a2y1evNiQZPz3v/+1rouKisq0zkGDBhleXl4Z/nzdzZUrV4xKlSoZkoygoCCjZ8+exqeffmqcOXMmQ9vbv6Njx44Zjo6Oxvjx423a7dmzxyhSpIh1fXp6ulGxYkUjIiLC5ju7cuWKERwcbLRs2TLDOZ588kmbY/bv39+QZOzatStH1wcAmeHRNQCmMnXqVK1cuTLDp2bNmjbtli1bJkkaMmSIzfpXX31Vku5r3Eu/fv3k6Ohos87Hx0e//PJLjma3SkxMlCR5enpmq31uXtPFixe1Zs0adenSRZcvX9b58+d1/vx5XbhwQRERETp06JD++OMPSdIPP/yg8PBwhYWFWfcvXry4unfvbnPMlStXKj4+Xt26dbMe7/z583J0dFSDBg20du3aO9ZUrVo1xcbG6vnnn9exY8f0wQcfqGPHjipVqpT+85//ZPvasmPVqlVKTU3VgAEDbHo/Bg8enKHt119/rcaNG6tYsWI219WiRQulpaVpw4YNNu27du2qYsWKWZcbN24sSTp69Ohd6/Lx8VFycrJWrlyZo+txc3PTL7/8oqFDh0q6+ahenz59VLp0aQ0YMOCOjzZ+9913Sk9PV5cuXWyuz9/fXxUrVrTet9jYWB06dEjPPfecLly4YG2XnJysxx9/XBs2bLB5lE+6Oc7s7wYMGCDprz/LAHA/eHQNgKk8/PDDqlevXob1t34JveX48eNycHBQaGioTTt/f3/5+Pjo+PHj91xDcHBwhnXvvvuuevToocDAQNWtW1dt2rRRZGSkKlSokOVxbo2buHz5crbOm5vXdPjwYRmGoREjRmjEiBGZtjl79qzKlCmj48ePKzw8PMP22+s4dOiQJGX52Fl2xok89NBDmjNnjtLS0rR//34tXbpU7777rl544QUFBwerRYsWdz1Gdtz6ripWrGiz3tfX1yakSDeva/fu3fL19c30WLdPlFCuXDmb5VvHu3Tp0l3r6t+/v7766iu1bt1aZcqU0RNPPKEuXbqoVatWd93X29tb7777rt59910dP35cq1ev1sSJE/XRRx/J29tb48aNy3S/Q4cOyTCMDN/FLU5OTtZ2ktSjR48sa0hISLD5/m4/ZkhIiBwcHDKM6wKAe0HQAVCo3c+LEdPS0jJd7+bmlmFdly5d1LhxYy1cuFA//vij3nvvPU2YMEHfffedWrdunelxvLy8FBAQoL179+aortx42eOtf3l/7bXXFBERkWmb24NMdo85Z84c+fv7Z9hepEj2/0pydHRUjRo1VKNGDYWHh6t58+aaN2/eXYNOVt9NVvcyO9LT09WyZUu9/vrrmW5/6KGHbJZv7+27xTCMu57Lz89PsbGxWrFihZYvX67ly5dr1qxZioyM1Oeff57tmoOCgtS7d2899dR
      "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",
    "# Generate random samples for histograms\n",
    "prior_samples = norm.rvs(loc=mu0, scale=sigma0, size=1000)\n",
    "likelihood_samples = norm.rvs(loc=ybar, scale=sigma / np.sqrt(n), size=1000)\n",
    "posterior_samples = norm.rvs(loc=mu_n, scale=sigma_n, size=1000)\n",
    "\n",
    "# Plot histograms\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.hist(prior_samples, bins=30, density=True, alpha=0.5, label=\"Prior\")\n",
    "plt.hist(likelihood_samples, bins=30, density=True, alpha=0.5, label=\"Likelihood\")\n",
    "plt.hist(posterior_samples, bins=30, density=True, alpha=0.5, label=\"Posterior\")\n",
    "plt.title(\"Hours College Students Sleep\")\n",
    "plt.xlabel(\"Sleep Time (hours)\")\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              Prior  Data (Sample Mean)  Posterior\n",
      "0  Precision    1.0              44.444     45.444\n",
      "1         SD    1.0               0.150      0.148\n",
      "2       Mean    5.5               4.844      4.859\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))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}