cs146-3.1-pre-class-work.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normal likelihoods and normal-inverse-gamma priors\n",
    "\n",
    "Today we explore how samples from a prior distribution can be interpreted as instances of the likelihood function. Specifically, we look at how samples from a normal-inverse-gamma (NIG) distribution can be interpreted as normal distributions.\n",
    "\n",
    "**In short:** Each sample from the NIG distribution is a pair $(x, \\sigma^2)$. These values specify the mean and variance of a normal distribution and so we can think of the sample (the pair of values) as an instance of the normal distribution (which will be our likelihood function). More below.\n",
    "\n",
    "## Normal-inverse-gamma in SciPy\n",
    "\n",
    "Even though SciPy does have classes defined for the normal distribution (`scipy.stats.norm`) and the inverse-gamma distribution (`scipy.stats.invgamma`), it does not have one defined for the normal-inverse-gamma distribution. To help you, the functions below implement the probability density function and a sampler for the normal-inverse-gamma distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sts\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Function definitions for the normal-inverse-gamma distribution. The parameters\n",
    "of the distribution, namely mu (μ), either lambda (λ) or nu (ν), alpha (α),\n",
    "beta (β), are used as defined here:\n",
    "\n",
    "  https://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution\n",
    "\n",
    "Note that we use the symbol nu (ν) rather than lambda (λ) for the third\n",
    "parameter. This is to match the notation used in the conjugate priors table on\n",
    "Wikipedia:\n",
    "\n",
    "  https://en.wikipedia.org/wiki/Conjugate_prior#Table_of_conjugate_distributions\n",
    "'''\n",
    "\n",
    "def norminvgamma_pdf(x, sigma2, mu, nu, alpha, beta):\n",
    "    '''\n",
    "    The probability density function of the normal-inverse-gamma distribution at\n",
    "    x (mean) and sigma2 (variance).\n",
    "    '''\n",
    "    return (\n",
    "        sts.norm.pdf(x, loc=mu, scale=np.sqrt(sigma2 / nu)) *\n",
    "        sts.invgamma.pdf(sigma2, a=alpha, scale=beta))\n",
    "\n",
    "def norminvgamma_rvs(mu, nu, alpha, beta, size=1):\n",
    "    '''\n",
    "    Generate n samples from the normal-inverse-gamma distribution. This function\n",
    "    returns a (size x 2) matrix where each row contains a sample, (x, sigma2).\n",
    "    '''\n",
    "    # Sample sigma^2 from the inverse-gamma distribution\n",
    "    sigma2 = sts.invgamma.rvs(a=alpha, scale=beta, size=size)\n",
    "    # Sample x from the normal distribution\n",
    "    x = sts.norm.rvs(loc=mu, scale=np.sqrt(sigma2 / nu), size=size)\n",
    "    return np.vstack((x, sigma2)).transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1\n",
    "\n",
    "1. Generate 10 samples from the normal-inverse-gamma (NIG) distribution with parameters as\n",
    "   provided below.\n",
    "   \n",
    "   Each sample corresponds to the mean and variance of a normal\n",
    "   distribution.\n",
    "   \n",
    "   With these NIG parameters, the prior 95% confidence interval for\n",
    "   the mean is about [-10, 10] and for the variance [0.1, 10].\n",
    "   \n",
    "   In practice you would\n",
    "   work the other way around: use confidence intervals (or other information) to determine values for the\n",
    "   prior hyperparameters.\n",
    "\n",
    "\n",
    "2. Plot the 10 normal distributions corresponding to your 10 samples. To see the functions\n",
    "   clearly, plot your graphs on the domain [-15, 15].\n",
    "   \n",
    "   You should see that the 10 samples\n",
    "   (normal distributions) are all quite different. This means the prior is quite broad\n",
    "   (uncertain) over the mean and variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normal-inverse-gamma prior hyperparameters\n",
    "mu_0 = 0        # The prior mean is centered around 0.\n",
    "nu_0 = 0.054    # The smaller ν₀ is, the more uncertain we are about the prior mean.\n",
    "alpha_0 = 1.12  # α₀ and β₀ control the marginal prior over the variance.\n",
    "beta_0 = 0.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 9.5444441   3.67203472]\n",
      " [-4.03010843  0.70490904]\n",
      " [-4.06012632  0.32924424]\n",
      " [ 2.68216204  0.58756206]\n",
      " [ 4.62908863 11.66975558]\n",
      " [-0.31906779  0.85431638]\n",
      " [ 5.19282216  1.49825748]\n",
      " [ 3.74791633  1.3569993 ]\n",
      " [-3.86371762  0.66381644]\n",
      " [-1.55328481  0.16585573]]\n"
     ]
    }
   ],
   "source": [
    "# Generate the samples of loc and sigma\n",
    "samples = norminvgamma_rvs(mu_0, nu_0, alpha_0, beta_0, size=10)  \n",
    "\n",
    "# Plot the graph with the ten normal distributions\n",
    "x = np.linspace(-15, 15, 300)\n",
    "distributions = [sts.norm.pdf(x, loc=samples[i][0], scale=1/samples[i][1]) for i in range(len(samples))]\n",
    "for normal_distribution in distributions:\n",
    "    plt.plot(x, normal_distribution)\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Probability')\n",
    "plt.title('Normal distributions')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd5xb9Z3v/9fnqIw0murxuDdsHMCmYwiJaQESwBDYZZMlhWw2jc0jyd5symaT3fTsbvLb3OQmuWmXSwhcUiBtCSUkNBtTg20cqjHYxr3PePqM2vn8/jhHY814iqbIso4+z8dDD0mn6Xs00nu++p7v+R5RVYwxxgSPU+oCGGOMKQ4LeGOMCSgLeGOMCSgLeGOMCSgLeGOMCSgLeGOMCSgLeBNIIvJlEfnZGJZXETnef/xjEfnCJJVjnoh0iUjIf75KRD44Gdv2t3efiLx3srZngsUC3oyLiGwVkX0iksib9kERWVXCYk0KVf2wqn5ttOX89+DSUba1XVVrVDU70XIN9U9LVa9Q1Vsnum0TTBbwZiLCwMcnuhHxBO6zKCLhUpfBVLbAfanMUfVN4NMi0jDUTBF5o4isEZF2//6NefNWich/iMjjQA+w0J/27yLyhN+scbeINInIz0Wkw9/GgrxtfFdEdvjz1onI+YUWXET+WUT2iMhuEXn/oHm3iMi/+4+nisg9ItImIq0i8qiIOCJyGzAPuNsv62dEZIHf1PMBEdkOPJw3LT/sF4nI0/778nsRmeK/1kUisnNQWbaKyKUicjnwr8B1/us9m/c+ftB/7IjI50Vkm4jsF5H/JyL1/rxcOd4rIttF5KCI/Fve65wjImv993KfiHy70PfSHLss4M1ErAVWAZ8ePMMPrXuB7wFNwLeBe0WkKW+x9wA3ALXANn/aO/zps4FFwJPAT4EpwAbgS3nrrwFO9+f9Avi1iMRGK7Qflp8G3gwsBkZqZvkUsBNoBqbjhayq6nuA7cBb/SaY/8pb50LgJOCyYbb5d8D7gVlABu89GpGq/hH4T+AO//VOG2Kxv/dvbwIWAjXA9wctcx5wAnAJ8EUROcmf/l3gu6pah/e+/2q0MpljnwW8magvAv8oIs2Dpl8JvKqqt6lqRlV/CbwMvDVvmVtU9UV/ftqf9lNV3ayq7cB9wGZVfVBVM8CvgTNyK6vqz1S1xV//W0AVXniN5m/913lBVbuBL4+wbBqYCcxX1bSqPqqjD+D0ZVXtVtXeYebflvfaXwD+NncQdoLeDXxbVbeoahfwOeAdg349fEVVe1X1WeBZIPePIg0cLyJTVbVLVZ+ahPKYErOANxOiqi8A9wCfHTRrFodr5Tnb8GrmOTuG2OS+vMe9QzyvyT0RkU+JyAa/qaMNqAemFlDsWYNee3A5830T2ATcLyJbRGTwfg5lqP0abv42IEJh5R7N4Pd8G95xkul50/bmPe7h8Pv5AeB1wMt+U9hVk1AeU2IW8GYyfAn4EAPDezcwf9By84Bdec/HPZSp397+L3i18UZVbQDaASlg9T3A3EHlGpKqdqrqp1R1Id6vj0+KyCW52cOtNsrrD37tNHAQ6AaqczP8Wn3+L6PRtjv4PZ+H1wS0b+jF8zas+qqqvhOYBvx/wG/ye0iZ8mQBbyZMVTcBdwD/I2/yH4DXici7RCQsItcBS/Bq+5OhFi+8DgBhEfkiUFfgur8C/l5ElohINQPb9QcQkatE5HgREaADyPo38IJz4TjKfn3ea38V+I3fjfIVICYiV4pIBPg8XrNTzj5gwQg9jn4JfEJEjhORGg632WdGK5CIXC8izarqAm3+5Al37TSlZQFvJstXgf4an6q2AFfhHaRsAT4DXKWqByfp9f6E10b/Cl5TRB+jN43kynYf8B3gYbzml4dHWHwx8CDQhXfA94equsqf93Xg834PmyMONI/gNuAWvOaSGP4/Rv+4w0eAm/B+6XTjHeDN+bV/3yIizwyx3Zv9ba8GXsN7T/6xwDJdDrwoIl14B1zfoap9he+SORaJXfDDGGOCyWrwxhgTUBbwxhgTUBbwxhgTUBbwxhgTUMfUYEhTp07VBQsWlLoYxhhTNtatW3dQVQefSQ4cYwG/YMEC1q5dW+piGGNM2RCRYc/EtiYaY4wJKAt4Y4wJKAt4Y4wJKAt4Y4wJKAt4Y4wJKAt4Y4wJKAt4Y4wJKAt4YwKstyvFpnX7S10MUyLH1IlOxpjJdd+PnmfP5nZmLT6P6rpoqYtjjjKrwRsTYG37e0pdBFNCFvDGBFg2413QZ9iL/JlAsz+7MQHmZlzvgV24rSJZwBsTYNmsl+x2Zc7KZAFvTICpmwt4S/hKZAFvTCWwfK9IFvDGVACrwVcmC3hjKoDle2UqesCLSEhE1ovIPcV+LWPM0HJt8aayHI0a/MeBDUfhdYwxw7AafGUqasCLyBzgSuCmYr6OMWZk1gZfmYpdg/8O8BnAHW4BEblBRNaKyNoDBw4UuTjGVCjL94pUtIAXkauA/aq6bqTlVPVGVV2mqsuam5uLVRxjKprV4CtTMWvwy4GrRWQrcDtwsYj8rIivZ4wZhg77G9oEWdECXlU/p6pzVHUB8A7gYVW9vlivZ4wZntXgK5P1gzfGmIA6Khf8UNVVwKqj8VrGmCNZDb4yWQ3emIDKD3Vrg69MFvDGBJSbzQt4q8FXJAt4YwIqmz5cbbd8r0wW8MYEVGZAwFvCVyILeGMCKpvJa3i3fK9IFvDGBNSAJhobTbIiWcAbE1D5NXhroalMFvDGBNTAgLeEr0QW8MYEVMZ60VQ8C3hjAiq/Dd4SvjJZwBsTUNYGbyzgjQmorPWDr3gW8MYElLXBGwt4YwJqQBON9YOvSBbwxgTUwIOspSuHKR0LeGMCyvrBGwt4YwJq4HDBJSyIKRkLeGMCasAFPyzhK5IFvDEBlX8VJ7uiU2WygDcmoKzWbizgjQmo/Hy3sK9MFvDGBJS1wRsLeGOCKr8Gb23wFckC3piAGnD2qtXgK5IFvDEBpZbvFc8C3piAsjZ4YwFvTEBZDd5YwBsTUANq8DaaZEWygDcmqOyKfRXPAt6YgLJ2d2MBb0xADWiDtyaaimQBb0xADexFU8KCmJKxgDcmoFTBccR/bAlfiSzgjQkoVUVC4j8pbVlMaVjAGxNUrloNvsJZwBsTUKrg+DV4G2ysMlnAGxNQqopYDb6iFS3gRSQmIk+LyLMi8qKIfKVYr2WMOZIdZDXhIm47CVysql0iEgEeE5H7VPWpIr6mMcY3sAZf4sKYkihawKtXZejyn0b8m33MjDlK1D1cg7eEr0xFbYMXkZCI/AXYDzygqn8eYpkbRGStiKw9cOBAMYtjTGXJ6yZp+V6ZihrwqppV1dOBOcA5InLyEMvcqKrLVHVZc3NzMYtjTEUZ2AZf4sKYkjgqvWhUtQ1YBVx+NF7PGGO9aExxe9E0i0iD/zgOXAq8XKzXM8YMpAoi3s2OflWmYvaimQncKiIhvH8kv1LVe4r4esaYPOoqIoKI2GiSFaqYvWieA84o1vaNMSNTBXEAx9rgK5WdyWpMUGleDd4SviJZwBsTUPlt8JbvlckC3piAUj/hrQZfuSzgjQmoXA0eYcAFuE3lsIA3JqAG9KKxfpIVyQLemIDK9aKxNvjKZQFvTFDl9aLB+sFXJAt4YwKqvxeN9YOvWBbwxgSU5vWTtF40lckC3piAUtf6wVc6C3hjAio3mqT1g69cBQW8iPxWRK4UEfuHYEyZsDNZTaGB/SPgXcCrIvINETmxiGUyxkwG60VT8QoKeFV9UFXfDZwJbAUeEJEnROR9/gW1jTHHmPwzWa0GX5kKbnIRkSbg74EPAuuB7+IF/gNFKZkxZkIOj0WDnclaoQoaD15EfgecCNwGvFVV9/iz7hCRtcUqnDFm/A73ohHUxqKpSIVe8OMmVf1D/gQRqVLVpKouK0K5jDET1N+LxrE2mkpVaBPNvw8x7cnJLIgxZnJZLxozYg1eRGYAs4G4iJyBN/AoQB1QXeSyGWMmwu9FY2eyVq7RmmguwzuwOgf4dt70TuBfi1QmY8wksBq8GTHgVfVW4FYR+RtV/e1RKpMxZhKom3dFJ+sHX5FGa6K5XlV/BiwQkU8Onq+q3x5iNWPMMcA7yGqjSVay0ZpoEv59TbELYoyZXF4TjX8mqyV8RRqtieb/+PdfOTrFMcZMFlX1zmTF8r1SjdZE872R5qvq/5jc4hhjJk1/DR7rRVOhRmuiWXdUSmGMmXTeRbfxT3QqdWlMKRTSi8YYU4ZUAcdq8JVstCaa76jqP4nI3QxRB1DVq4tWMmPMhGjecMGW75VptCaa2/z7/1nsghhjJldusDGsBl+xRmuiWeffPyIiUbwRJRXYqKqpo1A+Y8w4DajB22iSFanQ4YKvBH4MbMYbj+Y4EfkHVb2vmIUzxoxf/1AFDrhZq8FXokKHC/4W8CZV3QQgIouAewELeGOOUfkX3bZG+MpU6HDB+3Ph7tsC7C9CeYwxk2TgYGMW8JVotF401/oPXxSRPwC/wmuDfzuwpshlM8ZMhGu9aCrdaE00b817vA+40H98AGgsSomMMZNC1Wt/R7DRJCvUaL1o3ne0CmKMmVyHL7otoy9sAqnQXjQx4APAUiCWm66q7y9SuYwxE2QX/DCFHmS9DZiBd4WnR/Cu8NRZrEIZYyZO7ZJ9Fa/QgD9eVb8AdPvj01wJnDLSCiIyV0RWisgGEXlRRD4+0cIaYwqjqv5okn4N3k50qkiF9oNP+/dtInIysBdYMMo6GeBTqvqMiNQC60TkAVV9aXxFNcYUzK+wiyOIYzX4SlVowN8oIo3AF4C78K7w9IWRVlDVPcAe/3GniGwAZgMW8MYUWS7QrQ2+shUU8Kp6k//wEWDhWF9ERBYAZwB/HmLeDcANAPPmzRvrpo0xQ+gPdLtkX0UrqA1eRJpE5H+LyDMisk5EviMiTQWuWwP8FvgnVe0YPF9Vb1TVZaq6rLm5eWylN8YMyWrwBgo/yHo73tAEfwO8DTgI3DHaSiISwQv3n6vq78ZbSGPM2OQCvb8XjZ3oVJEKbYOfoqpfy3v+7yLyVyOtIN7ZFT8BNqjqt8dbQGPM2OUC3a7JWtkKrcGvFJF3iIjj3/4WbzTJkSwH3gNcLCJ/8W8rJlRaY0xh+nvRYGeyVrDRBhvrxPuoCPBJ4Gf+LAfoAr403Lqq+pi/njHmKDvcBi/WD76CjTYWTe3RKogxZvIc7kUDWD/4ilVoGzwicjVwgf90lareU5wiGWMm6ogavOV7RSq0m+Q3gI/jnaT0EvBxf5ox5hiUa5LxOtFYDb5SFVqDXwGcrup9bETkVmA98NliFcwYM379NXjHavCVrNBeNAANeY/rJ7sgxphJ1N8PHjuTtYIVWoP/OrBeRFbiHba5APhc0UpljJkQ60VjoICA909Yegw4FzgbL+D/RVX3Frlsxphxsl40BgoIeFVVEblTVc/CG0nSGHOMG3AmK9ZCU6kKbYN/SkTOLmpJjDGTRgedyWo1+MpUaBv8m4APi8hWoBv8SoHqqcUqmDFm/PLb4HHoP+hqKkuhAX9FUUthjJlcg3rRWA2+Mo02Fk0M+DBwPPA88BNVzRyNghljxs/OZDUwehv8rcAyvHC/AvhW0UtkjJmw/m6RVoOvaKM10SxR1VMAROQnwNPFL5IxZqIG1+CxfvAVabQafDr3wJpmjCkfA3rRWD/4ijVaDf40EcldR1WAuP8814umrqilM8aMy4BeNNYGX7FGGw8+dLQKYoyZRIPHosELfbu6U2UZy2BjxpgycUQbPFaLr0QW8MYE0MBeNP60Mk74rmSG1w52l7oYZafgKzoZY8pH/njwhxO+hAWagK5khtO+cj9ZV3nt6yusmWkMrAZvTADpgDb43LTyTPgfrtxE1h88rbU7VeLSlBcLeGMCaGAbfO4gaylLNH6bD3T1P97W2lPCkpQfC3hjgqg/4PN60bjlmfB72/uY3RAHYIcF/JhYwBsTQIcvui1I7ltenvnOnvY+li1oBGBbiwX8WFjAGxNAhw+yDuwHX27SWZcDXUkWNCWYURezgB8jC3hjAujwJfu8M1kHTCsjBzqTqMKM+hjzplRbE80YWcAbE0BDHmQtwzb4Pe19gB/wTdVsa7W+8GNhAW9MAOV3k3RCXsC7ZRjwe/2An1kfY2Z9jAOdybLcj1KxgDcmgPIvut0f8NnyC8Y97b0AzKyL01gdxVXo6EuPspbJsYA3Jojyhgsu54Df19FHLOJQFw/TmIgAdrLTWFjAGxNA+W3wjpML+PK76kdLV4qpNVWICA3VUQAO9VgNvlAW8MYE0OFeNOCEvK95Odbg23rTNPrBPsW/b+uxGnyhLOCNCaABNfgyPsja1pOiodprmskFvTXRFM4C3pgAOnwmK3lNNGUY8L1p6uNewDf4bfBt1kRTMAt4YwJoyBp8GQZ8e0+6vwZfWxUm7AiHrImmYBbwxgTREL1o1C2vg6yqSltvmoa41zSTO9BqAV84C3hjAmioGny2zGrwXckMWVf7a/AAjdURDnVbE02hihbwInKziOwXkReK9RrGmKEFoRdNrq091wYP0JiwGvxYFLMGfwtweRG3b4wZxoCxaPyDrFpmAd/e6wV8rv87+DV4C/iCFS3gVXU10Fqs7Rtjhnd4qILyPZM1V4PPb6KZkojaiU5jUPI2eBG5QUTWisjaAwcOlLo4xgRC/2BjTn4bfHkdZG3r9WrqDXlNNA3VUQ51p8pybPtSKHnAq+qNqrpMVZc1NzeXujjGBMKAGrxTnsMF97fB59Xg6+MRMq7Sk8qWqlhlpeQBb4yZfIeHC5ayPciaa4PPP8iaq83n5pmRWcAbE0T9CV/ObfApqqMhqsKh/mn1FvBjUsxukr8EngROEJGdIvKBYr2WMWagXL47TvmORdPWkx7Q/g6HA96GKyhMuFgbVtV3FmvbxpiRDTzRyZtWbsMFt/Wmqc/rIgmH2+OtBl+YogW8MaZ0coONIeA4ZdoG35OmPj4wonI1+A4L+IJYG7wxAdRfg3ekbEeTbOtN9Y9Dk9PfRNNrJzsVwgLemAAa8qLb5RbweSNJ5tRUhQk5Yk00BbKANyaAyv2SfbmRJOsHBbyIUB+PWMAXyALemCDKq8GLI4iUVy+avrRLKuMe0UQDXjON9aIpjAW8MQGUX4MHb0TJcmqi6R+mYFANHqDOavAFs4A3JoD6hyXwv+FOSMor4HMDjcWPDPiGeMR60RTIAt6YAMpmdMA4NE5IyqqJZqhxaHLq4xHaLOALYgFvTABl0i6hiJPXRFNeNfj2/pEkh26DtyaawljAGxNAWT/gcxxHyqoXzVBjwec0VHtNNOX0i6RULOCNCaBsOks4fPjrLSEpqys6tfUOH/CN1VFcteEKCmEBb0wAZTKDavAhp6wuut3WkyYacohHQkfMm5Lwmm1a7dJ9o7KxaIwJIK+J5nA4hvw2eFUlubmdvpdbCU+Nk1g2HQkfe/W89t4U9dWR/mMI+XIBf6g7BXaNoBFZwBsTQNm0SzivBi+OoFmX9vteo2v1LnAEXKV77V6a338yzhBNIaXU2p2icZgy5QK+pdtq8KM59v51G2MmLJN2CYXzm2iExo4kXat3kXj9DGZ/5Y00XX8S6T3dtP7qlWPucn4tXSmaElVDzmvMr8GbEVkN3pgy09fXx8qVK3nxxRdRVZYuXcoll1xCVdXhQMwOaoOPiTCzrY+q4xtouOZ4xBHiJ0+l4cqFtN21mZ71+0mcNb0UuzOklu4US2fVDTlvSrW1wRfKavDGlJHu7m5uvvlmnn76aebPn89xxx3HmjVruOmmm+jp6elfbnATzdy+NKLQeO1ixDncrp04dyaROTW0/2kr7jF0IeuDXUmm1gxdg49HQ8QjIVq7LOBHYwFvTJnIZrPccccdtLa2cv311/P2t7+dt73tbVx//fW0trZy++23k816IZ3J6wefPthLcyrL/niY8JTYgG2KIzRcuRC3I0XP2n1HfZ+Gksq4dPZlaEoceZJTzpRE1GrwBbCAN6ZMPPbYY2zfvp2rr76aRYsW9U9ftGgR11xzDdu3b+eJJ54A/F40fht812O7UGB31dAtslXH1ROdV0vn47uOibb4Vr9tfUrN8AHfmIhYG3wBLOCNKQMtLS2sXr2apUuXcuqppx4x/9RTT2XJkiWsWrWKQ4cOkc14TTRuX4aeZ/ZzqDpMcoTsrjlvNtmWPvo2tBRxLwpzsCsJMOxBVoApiar+fwRmeBbwxpSB+++/n1AoxGWXXTbsMpdddhkiwsMPP9zfRNOzfj+aytJSWzXiqf3xpVMJNVTR+eiuYhR/THLdH6eOUIOfUh2xJpoCWMAbc4zbuXMnGzduZPny5dTVDd2zBKC+vp5zzz2X559/nl63nVBY6HpyD5E5NSSrIyMONiYhoWb5bFJbO0jt7CzGbhSstduvwQ9zkBW8rpKHum2ogtFYwBtzjFu5ciXxeJxzzz131GXPO+884vE4HbFNJHoyZPb3UHPuTH+wsZHb1xNnT0eiDl1P7Zmsoo9Li987pmmEGvzUmiq6khl6j6GeP8ciC3hjjmHbtm1j8+bNnHfeeQP6uQ8nFotx/vnnk4q20bVnOxIPU31asz9c8MijSTqxMNWnT6P32QO4Jbwk3sGuFJGQUDvMQWGAGXVeb6A97b1Hq1hlyU50MsbX2XKQ7S88y55XN9Lb2UEoEqF+2gzmn3o6sxafiBM6cuCrYnv44Yepqanh7LPPLnidM047i4fuW8ULPS9z8jlvRSKhgseDT7x+Jt1P76X7mf3Unjd7IkUft5auJE2JqiHHocmZ2ZAL+D4WNtccraKVHQt4U9FUlW3PrWfN3b9j+/N/AaCqOkF1fQOZdIqXH3uEp377S+K1dZx+2VWcueJqYomjEyivvfYa27Zt44orriAaHb65YjBRh/nJ+WypfpUDs1M04o8HX0AXyOjsGqLzaun+8x5qls8aMWSLZW9HH9PrRv61Mqs+DsDuNqvBj8QC3lSs1t27eOgnP2T7C8+SaJzCee/4OxaeeTZT585HHK/1MtnTzdZn17PhsZU8+ZtfsO7eOznvnX/H6W9e0b9MMagqK1eupLa2ljPPPHNM66b7Mrye2eyVHaxe9ziLTz9xTBfdTrx+Jod+/QrJze3Ejm8YT/EnZEdrD0tn14+4zIz6wzV4MzwLeHP0ZJKw9THYuQZ2/wU6d6PdB0lFHDQah7pZhJtOIjz3PFh4IVTVFqUYqsr6++5i9S9uIRyJcvH7/oFTLrmccOTI0QurqhOc8IbzOOEN57F/6xYe+dnNPHzzj3n58dWs+NgnqZ82oyhl3LJlC9u3b2fFihVEhijXSJIvt5JwQpw8+1Se3vk0mzZtGtMl+6pPbab93i10P7X7qAe86yq72nq5/OSZIy4Xi4RoSkStDX4UFvCmuFRh+5Ow9mZ04x/pDvfQ0hilo7mRjqmQDKVRUaAXaAVeILTjDhIbXepCM2mYvYKmJR8jHJsyKcVJ9fVy/4+/x8YnH2XRstfz5g99jERDY0HrTluwkLf929d4afXDrLz1Rm777Me54qOfZNFZr5+UsuWoKqtWraKurm7MtXdVJblmL91ZZdGpp7OxcyOrVq3ixIYLC77EnUQcqs+aTtfju8h2JAmN0lwymfZ19pHOKnMa46MuO7Mhxu42q8GPxALeFEc2Dc/9Cv78I/oOvcie2fXsOaeRXscLi1hsOnV1p1Adn0+0ajohJwYo6eRB+lqepSv0HHvc/exsuxV57FaawguZ+bp/ZOqMK3Cc8Y1d3rJrB3d96z85tHsX57/r7zn76r8ZcxuziLD0wkuYfeJS7v5fX+fO//oaZ1/zNs677j1jOgirqqR37KDvpZdI79xJev9+MvsP4Pb2sCscYUfzVM53ldYf/IDIrFnETjyRqsWLcWKxEbebeq0D3dfDpqTLkliUCy64gLvvvpvp0d1juuh2zbkz6XpsF52P7aZhxXEFrzdRO1q9GvncKdWjLjuzPs72lp5Rl6tkFvBmcrkuvPTf6MNf45C7k+0Lm2l5nVf7bmw8jXnTVjC16SJisVnDb8PPEzeTpH3DTRzcfAt7q17l4MufIPLyvzJ7zruZs+BDVEWnFlysV/78OH/84XcIR6O87fNfY97Jp01kL2mYPoN3fvWbrLz1Rtb8/jfs3fQKV338M1TXD92koaqkXnuN7kcfpevRx+h97jncjo7++U4iQbi5GWpqWLtoIYlkklmPPkZLa6v3ngI4DlUnnEDijW+gZvly4medhTOo62TnIzsgFmJ7W5pTww6nLz2dRx99lE17nyOaXVrw/oWb4sRPbab7qd3UXjiHUOLoXBBk5yEvsOcWUIOfVR/jqS2lH1rhWGYBbyaHKmx+GH3oyxxIv8zWRVPojNcTjdazYNZ1zJr5NuLxuWPapBOuovGUj9J4ykdZtO0JWtd8mV3uS2zVm9i+86fMnHEt8xZ8mOrqBcNuw81mefSXt7L27t8x8/gTuOoTn6Vu6uRc5y0cjfLmD32MWa87iQf/7w+47XP/xNWf+BwzF5/Qv0xq2zbaf38X7ffcQ3r7dgCiixZRd8UVxJYuIbZ0KVULFuAkEgCsXbuWlnvu4dprr+Wkr38dzWZJ79pF34aX6Xt5A73rnqH1/91G609uxqmupuaSS6i74goS5y0n25qmb+Mh9IxpuHt3EYo4hEIhLrroIu68805qY/tQVwcMFzySuovn0vvsAboe30X9WxZMyns2mlwNflbD6AE/p7Gazr4Mrd2p/qs8mYEs4M3E7VyH++AX2du7hm3z6+iJ1RGPz+bEeTcwY8ZfEwpNvA3Xmf9Gps6/n6k71tDz6BfYps+z2/01u/b+hmnNl7FgwUeprV0yYJ2u1hbu/d432bnhBU57y5Vc9HcfHPJA6kQtvfASps5bwF3f+k/u+PK/cOF172F+d4qOu+6md/16ECHxhnNpev/7SJx3PtE5Q/cv7+7u5qGHHmL+/PmccsopAEgoRHTePKLz5lF32VsAcLu76V6zhq6HHqbz/vvpuPtunLo6Ehd/GkIzyBxXDyt39Q8XfOqpp/LIg0/Q5m6hq7Ob2vrCunlGpieIn9xE1+O7qT1/Dk68+HGx41AP0+uqiA1xse3BlvgXBHKN2o0AABG8SURBVHlxdzvnL7aLsw7FAt6M376XyKz6Krs7VrJ9bg3JaC21NSdy8oKPMK35MkSKcGLQ3LOpftcfOem1R1n4yJfYwQZ2Zv/I/gN/pGnKBcxf8BEaG85m67PP8Ifvf4t0so8rPvYplpz/pskvS55ps+fyV1dcy/2338LDP/8pc1o7ODNWx7RPf4q6q64iMmP03jYPPvggfX19rFixYsRjA04iQe1FF1F70UXM+MLn6X7ySdr/8CTKDJIv/J49z/fArLfh7tsNC+pwHIfTFr2RVX/5PQ88eD/X/s21Be9X7cXz6H2hhY6V22lYsbDg9cbrxd0dLJ5WWO+ppf0B32EBPwwLeDN2rVtIrfoqW9rvZ8/sGG5zDVQtpqfuInYwncd2vEz7pj/Tne6mN9NLT7rHu8/00JfpI+tmcXFRVRTFVRdFCUuYqlAV0VCUWDhGNBSlKlRFIpygrqqOuqh/8x83XvoVmg9u5qQ1N9LtbGFn5lGeaV2NJOey6WGorjuJt37iGzTNGVvTUKFUlb5nn6X97nvouPdesm1tLGtqYutZJ/MCO+meNZvLzl9OUwHhvmHDBtavX8/y5cuZPr3wS+dJNEpi+fl0PVOD25Om6X0Xse8PzwKw+8M30Hf8DOqvvpq6+mXEu+fw3PPPsWTpEk488cSCth+dVUPi7Bl0PbabxJnTicxIFFy2sersS7Nxbwf/ePHigpZvqI4yuyHOi7s7Rl+4QlnAmyF1pbrY272XfT37Dt8f3EDvwaeZVdXKCXUuVQ1xXugN8WBLmK2pXcDPAQg7YRqrGklEEsTDceLhOI2xRmaFZxEPx4n4vWAccXDEa0YQhIxmSGVTJLNJktkkqWyKvkwfe7r3sPHQRjpSHXSnu48sbBTCMp0TXohwacph9kl7WHRFhqTTy5MHvkeDXs6MmvlMT0ynJlIz4bMzk5s3037PPXTccy/pHTuQqipqL7mY+muuIfHGN3JCJMIJz63n/v/zPW7/0mc4a8U1LL/ueiJVQ/eAOXDgAHfeeSczZ87kTW8a+y+Njge2kdnXQ9N7lxA/6Vwams+E216m+R8+ROZP/82+//gPWua+icSia6lpbuXOO+/kAx/4AM3NhdV66y5fQO+GFlpvf5lpHz0DiRTnBK9nd7TjKpw1v7Buq+DV4l/c1V6U8gSBBXwF6sv0sbd7L3t79rKnaw97e/Z6Id59OMy70l39y8dFOa06y/LqNHOnKVkVDoWOo6XhUubXnMg/x6bQUNVAY6yxP9iLdYp7xs3QmeqkPdnOoeQhDvQcYM/OLbT86Wn0lf0crIanQi7V8zs5q66FOe7vaG/5HQ/1hFnXHWK/W0Nz9XSmJ6Yzvdq7zUjM8B4npjOtehqNVY0Dyq+uS99LG+hatYrOhx4iuWEDOA6Jc89l6kc+Qu2bLyVUM7Bde8GpZ/Deb/6AR3/xU9bdeycbn1jNG97+bk6+6NIB3SlbW1v52c9+Rjgc5rrrriMcHttXsue5A3Q+spPEOTOIn9QEeFdzAmh+93XEP/we+ja+QtfPVyP7HU6+eyVPnH8at37/+1x3zjnMvuii/gO8wwklIjS+/XW0/PRFWn/7ClP+9oSCD9SOxbpthxCB0+cVfnLVKbPreWDDvhGv4VrJRLX0l+jKWbZsma5du7bUxShryWySlt4W9vfsZ2/3XvZ07xlwv7d7L4eSh45YrynW1B90M6qbmdV3kOmdz1OdeY10LeAIceqZMe99TJ3zLjRUT0qVlOvfVEm5LmlXSap6965LWvPne8ukXCWtSjL33J/nTfO2MXDb/jr+trOqXj/yVIq+7m5SyT4Qh6pEgmh1AgU01Y2b6iKjWTIhB3UEEFwRXA2RRcm4WbKaBRTURXABRVCqXIfqDMSSLrGeNKFMBsdVJF5FuKGe6qZmYrEqwiJERIg4A+/DeY97Ww+y69ln6Dmwn+rqOPNPWMLc151Ib3c3z6x5GslmuPSCC5k5tYmoI8Qch6gjVDkOVf6991yoEoewH649z+6n9Y5XiM6rpfmDpyD+JfrWP7CdJ367iQ995wKiMe8fRtv+Hn7+xad4w+vDhDf+gXscByebZfmfn2becQuoOf8Cqs86k9hJJyHDjHvTsXIHHX/aSvVZ02m89ngkNHk1eVXl6u8/TsZV7vv4+QWvt2l/J5d+ezX/tuIkPnRB8Y8RHItEZJ2qLhtyXjEDXkQuB74LhICbVPUbIy1vAe990DNKf6j1ZFO09nXSmuzgULKbtlQXrclOWpLtHOrroi3VxaFUNx3pHjrSvfRmMiBh1L8hYSKhaqoj9cQjtVRFaqkKJ4iEqgmHYjhShSD09R2iL9VJXzZFUiAlETKEyRAlKzEyUkVaYbJH386FYpV/H3WEqDj+vfc84gdbBCXb3UVf+yE6Dx4g09eLEw4zZcYsmmbNIRyNIgJelAuCIj0tSMsmtG0ryXCG3rhDssrBRUDChMNNqNaQ7RWynVmyHX04nb1I1sUVh3TYoT3h0JZw6KgWMiEHxAEJEXZi3nvoRHEkivjvuYu3/awKaVXSk/wVCwFRhWhaiTpCPBGhKuT0/0NIHkrSu7+XRSdPJRZyiIr3z2Hzk3tpmhJn8SlTSXd18tL69aSTSWa2tzHn1Vep6e4iCiRmzaJmwXyq58wmMWsWiblziNXWEhVIPr6H1OpdxGcmmHbVQmrm1uFMwq+1x149yPU/+TNfv/YU3nnOvDGte+0PH6e9N839n7iQUBF+WRzrShLw4nWheAV4M7ATWAO8U1VfGm6d8Qb8gVSarEJWlYwqrkJGlSzaP/3w/eBl8qeDm/fYW99b11vHD1/XpS+b8e7dLEk3S8rNknRdv4bq3dJ+TTYX1rlpGdf70mcUMn4QZBFcdXCL0PMkIkoElyhZwpohQoqwpghpkhApHEkTJkOENFFNEyNETXQKNXXHUR1rGhi6fq20yhEijkPV4BD2a5oREaLOwKD21vVCKLedXDhkMxnSfX2k+npJ9/XS3dZGV+tBOlsO0rJrBwe2bqF1907cbBYnFGbu0lNYcsHFLD77DURiMTSbRVMpNJXCTSa9++4e3I52sh0dZA8dIrv9ebJbniW5Zwud9d10zYfehZCZSf+VESQJ0e44cWcqicZFxKctIhqbTjJcQztRWjNJDva1s7+3nb29rezpPcC+7v3s69lHb2bQuCgKIQ0Rc+NM0xlMyTbT0D2VaE8CdQTp6SByYB+uCNmaOExrxGlsINTYQKi2HieewIlVQyiGZKK4PUK2zSV9KEtaID0rhju7mrQ4ZBTvF1FW2bG1g4wDtTOrSanS53/2OnszZFAyIe+3ymQJqxIF/x/04V8fEf+fzuDPx+HPhENEhM7eNI9s3E82q1x/zjyqQg4hgbD/Syh3C4n3GqFB09dvO8RNq7dw4eKpXH/OPGqrwnnreOs5IjiACDgIjkDIvxcOzx9wLwyYlls3lLcNgZKMuJmvVAH/BuDLqnqZ//xzAKr69eHWGW/Az3v4aVJy9E90CGmGEFnCHL4Pa4ZQ/nMyhMkOWDZ3C2lu3bxpZAmRwas/pwlrJm+ddP/0iL9ehLR3r2kikibq9BFx0t66/npDffwy6SqSvXX9t76eeno6m0n21oK/RrEb73RMrzDJXyI5vH+Ok6a65iDVNa3E4h1U5W6xLhxn5ChUBdcNoyqggqrjTxevsUf991IdVMEhRERDRAkj/ryR320hPz9cVTJumpTbh6tu/74MLBS44oIM3Kaog6iDAlnE+xRJhLRESTtR0qEYmf5fbmEykvfplPxPrvfpy0rEe5w3Lyt5n7z8dQasf/h5lhBZCeHieI/xHhejolMsoi7+X9v/7ea97/n3I01z1KXW7WLdm68c3+uPEPDFPMg6G9iR93wncMSoTCJyA3ADwLx5Y/tplvP2trtRB0KaRVDvI6KK4BLCRdQl5L/5jvo3XESVEN46jvrzcQm5rr+dLI7SPz23XkhzH0MFHerrOUQY6TDz86fLMNP7Jw01P4QQBmKo66BuiLQbJpkNoW7YC59MhGwqTjYdI5uKkU3HQI/8AlX5t6FKOt6wl2GCWfvnD36l3L0DIsioFx2TQRvK+3uI11iTe6xDvr8CLtA2i1QbpPwvnqCIujjhJKFIN6FoL6FID04ojThpJJTxH3v3oCCKiNv/2KslerXAEF5t1Rni/dD+e++RV+nyupCqumQ1S1YzpDRFRtPg1Zf719TBG3KyuM7AQbgUEBwkGwWk/59LGAiTJkYa6Pb+8Tnee3/EeyYD3ui8v5gg4uD4DWMD/ymN7x+zgh/6XthncciK43+jvcdZQrji9C/X/4/CX85b1/G/qXL4JuLHsDfdFRkw38VBJfd40DxxBj0fuJ6L41cecp/DoaI9N1/8+BCqMsW5glYxA36ov+wROaGqNwI3gleDH88Lfevar41nNWOMCbRiXpN1J5B/hskcYHcRX88YY0yeYgb8GmCxiBwnIlHgHcBdRXw9Y4wxeYrWRKOqGRH5GPAnvKbIm1X1xWK9njHGmIGKeiarqv4B+EMxX8MYY8zQitlEY4wxpoQs4I0xJqAs4I0xJqAs4I0xJqCOqdEkReQAsG2cq08FDk5icUopKPsSlP0A25djUVD2Aya2L/NVdcjB/Y+pgJ8IEVk73HgM5SYo+xKU/QDbl2NRUPYDircv1kRjjDEBZQFvjDEBFaSAv7HUBZhEQdmXoOwH2L4ci4KyH1CkfQlMG7wxxpiBglSDN8YYk8cC3hhjAqqsA15E3i4iL4qIKyLL8qYvEJFeEfmLf/txKctZiOH2xZ/3ORHZJCIbReSyUpVxPETkyyKyK+9vsaLUZRoLEbncf983ichnS12eiRCRrSLyvP93KKur24vIzSKyX0ReyJs2RUQeEJFX/fvGUpaxUMPsS1G+J2Ud8MALwLXA6iHmbVbV0/3bh49yucZjyH0RkSV4Y+kvBS4Hfuhf0Lyc/K+8v0XZjC7qv88/AK4AlgDv9P8e5exN/t+h3PqP34L3+c/3WeAhVV0MPOQ/Lwe3cOS+QBG+J2Ud8Kq6QVU3lrock2GEfbkGuF1Vk6r6GrAJOOfolq5inQNsUtUtqpoCbsf7e5ijTFVXA62DJl8D3Oo/vhX4q6NaqHEaZl+KoqwDfhTHich6EXlERM4vdWEmYKiLl88uUVnG62Mi8pz/07Qsfkb7gvDe51PgfhFZ51/svtxNV9U9AP79tBKXZ6Im/XtyzAe8iDwoIi8McRupJrUHmKeqZwCfBH4hInVHp8TDG+e+FHTx8lIaZb9+BCwCTsf7u3yrpIUdm2P+vR+j5ap6Jl6T00dF5IJSF8j0K8r3pKhXdJoMqnrpONZJAkn/8ToR2Qy8DijpgaXx7AtlcPHyQvdLRP4vcE+RizOZjvn3fixUdbd/v19E/huvCWqo41flYp+IzFTVPSIyE9hf6gKNl6ruyz2ezO/JMV+DHw8Rac4diBSRhcBiYEtpSzVudwHvEJEqETkOb1+eLnGZCuZ/8XL+Gu9gcrkIzIXjRSQhIrW5x8BbKK+/xVDuAt7rP34v8PsSlmVCivU9OeZr8CMRkb8G/jfQDNwrIn9R1cuAC4CvikgGyAIfVtWjclBjvIbbF1V9UUR+BbwEZICPqmq2lGUdo/8SkdPxmja2Av9Q2uIULmAXjp8O/LeIgPe9/4Wq/rG0RSqciPwSuAiYKiI7gS8B3wB+JSIfALYDby9dCQs3zL5cVIzviQ1VYIwxARXIJhpjjDEW8MYYE1gW8MYYE1AW8MYYE1AW8MYYE1AW8MYMQ0TmishrIjLFf97oP59f6rIZUwgLeGOGoao78E4h/4Y/6RvAjaq6rXSlMqZw1g/emBGISARYB9wMfAg4wx9Z0phjXlmfyWpMsalqWkT+Gfgj8BYLd1NOrInGmNFdgTfC38mlLogxY2EBb8wI/PFB3gycC3xi0KBQxhzTLOCNGYZ4I3P9CPgnVd0OfBP4n6UtlTGFs4A3ZngfArar6gP+8x8CJ4rIhSUskzEFs140xhgTUFaDN8aYgLKAN8aYgLKAN8aYgLKAN8aYgLKAN8aYgLKAN8aYgLKAN8aYgPr/Aauxzph9AP8tAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2\n",
    "\n",
    "Generate 1,000,000 samples from the normal-inverse-gamma prior above and calculate\n",
    "approximate 95% confidence intervals over the mean and the variance using the\n",
    "samples. You can use the `numpy.percentile` function for this.\n",
    "\n",
    "Your confidence intervals should approximately match the intervals [-10, 10] and [0.1, 10]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000000\n"
     ]
    }
   ],
   "source": [
    "sample_distributions = norminvgamma_rvs(mu_0, nu_0, alpha_0, beta_0, size=1000) \n",
    "samples_from_distributions = []\n",
    "\n",
    "for distribution in sample_distributions:\n",
    "    samples_from_distributions.extend(sts.norm.rvs(loc=distribution[0], scale=1/distribution[1], size=10) )\n",
    "\n",
    "lower_bound = np.percentile(sample_distributions, 2.5, interopolation='lower')\n",
    "upper_bound = np.percentile(sample_distributions, 97.5, interopolation='lower')\n",
    "print('The 95% confidence interval is from [{},{}]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-4.32875082982137, -3.7092704830926744, -3.822941067124121, -4.5938025771850635, -4.652843446433886, -3.586970797263276, -4.622715407643007, -4.415226151869573, -3.4745905003711, -4.220846121838704]\n"
     ]
    }
   ],
   "source": [
    "lst = []\n",
    "lst.extend(samples_from_distributions[0])\n",
    "print(lst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3\n",
    "Code the equations for calculating the posterior normal-inverse-gamma hyperparameters\n",
    "from the prior hyperparameters and data.\n",
    "\n",
    "The equations are found on found [Wikipedia](https://en.wikipedia.org/wiki/Conjugate_prior#Table_of_conjugate_distributions) and reproduced below, using $d_i$ for a datum.\n",
    "Note that $n$ is the number of data, and $\\overline{d}$ is the mean of the data.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  \\mu_{\\text{post}} &= \\frac{\\nu_0\\mu_0 + n\\overline{d}}{\\nu_0 + n} \\\\\n",
    "  \\nu_{\\text{post}} &= \\nu_0 + n \\\\\n",
    "  \\alpha_{\\text{post}} &= \\alpha_0 + \\frac{n}{2} \\\\\n",
    "  \\beta_{\\text{post}} &=\n",
    "     \\beta_0 +\n",
    "     \\frac{1}{2}\\sum_{i=1}^n (d_i-\\overline{d})^2 +\n",
    "     \\frac{n\\nu_0}{\\nu_0+n}\\frac{(\\overline{d}-\\mu_0)^2}{2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Once you have the update equations implemented, Bayesian inference is done!\n",
    "\n",
    " * The likelihood function is the normal distribution with unknown mean and variance.\n",
    " * The posterior distribution is of the same type as the prior – normal-inverse-gamma.\n",
    " * The posterior parameters below give you the exact posterior normal-inverse-gamma distribution.\n",
    " * No approximation or numerical integration is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array([1, 2, 3, 4])  # In class you will get a larger data set.\n",
    "                               # This is just to get you started.\n",
    "mu_post = ...\n",
    "nu_post = ...\n",
    "alpha_post = ...\n",
    "beta_post = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4 (optional)\n",
    "\n",
    "You are told that the prior information we used above is incorrect. Actually, the prior 95%\n",
    "confidence interval on the mean should be [-15, 15] and on the variance [0.5, 2]. So, the prior\n",
    "over the mean is less certain (broader) than we had before, but the prior over the variance is\n",
    "more certain (narrower).\n",
    "\n",
    "Determine prior hyperparameters for the normal-inverse-gamma distribution that match the\n",
    "prior information above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}