CS146 Session 4.1 Pre-Class Work

CS146 Session 4.1 Pre-Class Work.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CS146 Session 4.1 Pre-Class Work\n",
    "\n",
    "In the previous session, you were told which parameters to use for the prior normalinverse-gamma distribution. This time you have to choose the prior parameters yourself.\n",
    "\n",
    "Given the information below, find reasonable values for the prior parameters of the normal-inverse-gamma distribution — $μ₀, ν₀, α₀, β₀$. You will be asked to provide your values for the prior hyperparameters in class, and to explain how you came up with them.\n",
    "\n",
    "Information provided:\n",
    "- The data are normally distributed. The error margins given below represent 1 standard deviation from the mean of the parameter.\n",
    "- Constraint: the mean of the data is approximately 2.3 ± 0.5.\n",
    "- Constraint: the variance of the data is approximately 2.75 ± 1.\n",
    "- Find $μ₀, ν₀, α₀, β₀$ hyperparameters for the normal-inverse-gamma prior that match this information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as sts\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from scipy.optimize import minimize\n",
    "from statistics import variance "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norminvgamma_pdf(x, sigma2, mu, nu, alpha, beta):\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",
    "    # 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": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With the hyperparameters from 3.1:\n",
      "Started with mu at 0.000, nu at 0.054; alpha at 2.100;  beta at 0.400\n",
      "Ended with mu at 2.300, nu at 5.500; alpha at 9.562;  beta at 23.547\n",
      "With randomly chosen hyperparameters:\n",
      "Started with mu at 10.000, nu at 10.000; alpha at 20.000;  beta at 15.000\n",
      "Ended with mu at 2.300, nu at 5.500; alpha at 9.562;  beta at 23.547\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    mu, nu, alpha, beta = x[0], x[1], x[2], x[3] # Setting up the prior hyperparameters\n",
    "    x = mu # Equation for x taken from Wikipedia\n",
    "    sigma_2 = beta/(alpha - 1) # Equation for sigma2 taken from Wikipedia\n",
    "    var_x = beta/((alpha - 1)*nu) # Equation for var_x taken from Wikipedia\n",
    "    var_sigma_2 = (beta**2) / (((alpha - 1)**2)*(alpha - 2)) # Equation for var_sigma2 taken from Wikipedia\n",
    "    return (x - 2.3)**2 + (sigma_2 - 2.75)**2 + (var_x - 0.5)**2 + (var_sigma_2 - 1)**2\n",
    "\n",
    "print('With the hyperparameters from 3.1:')\n",
    "x_initial = ([0, 0.054, 2.1, 0.4]) # Using the normal-inverse-gamma prior hyperparameters from 3.1, with the value of alpha bigger than 2 to prevent any error values\n",
    "result = minimize(f, x_initial)\n",
    "x_final = result.x\n",
    "print('Started with mu at %.3f, nu at %.3f; alpha at %.3f;  beta at %.3f' % (x_initial[0], x_initial[1], x_initial[2], x_initial[3]))\n",
    "print('Ended with mu at %.3f, nu at %.3f; alpha at %.3f;  beta at %.3f' % (x_final[0], x_final[1], x_final[2], x_final[3]))\n",
    "\n",
    "print('With randomly chosen hyperparameters:')\n",
    "x_initial2 = ([10, 10, 20, 15]) # Using randomly chosen normal-inverse-gamma prior hyperparameters to make sure there's no local minima\n",
    "result2 = minimize(f, x_initial2)\n",
    "x_final2 = result.x\n",
    "print('Started with mu at %.3f, nu at %.3f; alpha at %.3f;  beta at %.3f' % (x_initial2[0], x_initial2[1], x_initial2[2], x_initial2[3]))\n",
    "print('Ended with mu at %.3f, nu at %.3f; alpha at %.3f;  beta at %.3f' % (x_final2[0], x_final2[1], x_final2[2], x_final2[3]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}