{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 1\n",
    "\n",
    "# Exercise 1\n",
    "$n$ data points $y=(y_1,\\ldots,y_n)$ are drawn from a Gaussian distribution with known dispersion $\\sigma$ and unknown mean $\\mu$. Assuming a flat prior on $\\mu$, what is the expression for the posterior probability distribution on the mean parameter of the Gaussian distribution from which the data is generated, $P(\\mu|y)$? Pick values for the true mean and dispersion, simulate a large number of inferences, then calculate the fraction of times the true value of the mean lies within the 68% credible region of the inference.\n",
    "\n",
    "# Exercise 2\n",
    "Repeat the exercise above, but now assuming a Gaussian prior on the mean of the distribution. Observe how the answer to the question above is sensitive to a) the mean and dispersion of the Gaussian prior, b) the number of data points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "### Exercise 1\n",
    "\n",
    "Assuming a flat prior on $\\mu$, the posterior probability distribution for the mean of a Gaussian with known dispersion is a Gaussian with mean $\\mu_n$ and dispersion $\\sigma_n$,\n",
    "$P(\\mu|y) \\propto \\dfrac{1}{\\sqrt{2\\pi}\\sigma_n}\\exp{\\left\\{-\\dfrac{(\\mu - \\mu_n)^2}{2\\sigma_n^2}\\right\\}}$,\n",
    "with\n",
    "\n",
    "$\\mu_n \\equiv \\bar{y} = \\dfrac{1}{n}\\sum_i^n y_i$\n",
    "\n",
    "and\n",
    "\n",
    "$\\sigma_n^2 \\equiv \\dfrac{\\sigma^2}{n}$.\n",
    "\n",
    "Let's fix $\\sigma = 1$, $n=10$.\n",
    "Then\n",
    "\n",
    "1. Draw a random value for the true mean of the distribution, $\\mu_{true}$.\n",
    "2. Generate the data: draw $n$ data points from a Gaussian with mean $\\mu_{true}$ and dispersion $\\sigma$.\n",
    "3. Calculate $\\mu_n$ and $\\sigma_n$. Verify if $\\mu_n - \\sigma_n < \\mu_{true} < \\mu_n + \\sigma_n$.\n",
    "4. Repeat for a large number of simulations. Count the fraction of times condition 3 is satisfied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "n = 10\n",
    "sigma = 1.\n",
    "\n",
    "nsim = 10000 # number of simulations\n",
    "\n",
    "n_enclosed = 0\n",
    "\n",
    "for i in range(nsim):\n",
    "    mu_true = 2.*np.random.rand(1) -1. # draw a random value of the true mean of the distribution\n",
    "    y = np.random.normal(mu_true, sigma, n) # generate data\n",
    "    \n",
    "    ybar = y.mean()\n",
    "    mu_n = ybar\n",
    "    sigma_n = sigma/n**0.5\n",
    "    \n",
    "    if (mu_true > mu_n - sigma_n) and (mu_true < mu_n + sigma_n):\n",
    "        n_enclosed += 1\n",
    "        \n",
    "print '%3.2f'%(n_enclosed/float(nsim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "If we assume a Gaussian prior on $\\mu$, centered on $\\mu_0$ and with dispersion $\\tau_0$,\n",
    "$P(\\mu) = \\dfrac{1}{\\sqrt{2\\pi}\\tau_0}\\exp{\\left\\{-\\dfrac{(\\mu - \\mu_0)^2}{2\\tau_0^2}\\right\\}}$\n",
    "\n",
    "the posterior probability distribution is still a Gaussian, but the mean and dispersion change as follow:\n",
    "\n",
    "$\\sigma_n^2 \\equiv \\dfrac{1}{\\dfrac{1}{\\sigma^2} + \\dfrac{1}{\\tau_0^2}}$\n",
    "\n",
    "$\\mu_n \\equiv \\left(\\dfrac{n\\bar{y}}{\\sigma^2} + \\dfrac{\\mu_0}{\\tau_0^2}\\right)\\sigma_n^2$,\n",
    "\n",
    "which, in the limit $\\tau_0 \\rightarrow \\infty$, reduce to the expression valid for a flat prior.\n",
    "\n",
    "Let's fix $\\mu_0 = 0$ and $\\tau_0 = 3$. Since the prior is broad ($\\tau_0 > \\sigma/\\sqrt{n}$), this should not change the answer to the question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_enclosed = 0\n",
    "\n",
    "mu_0 = 0.\n",
    "tau_0 = 3.\n",
    "\n",
    "for i in range(nsim):\n",
    "    mu_true = 2.*np.random.rand(1) -1. # draw a random value of the true mean of the distribution\n",
    "    y = np.random.normal(mu_true, sigma, n) # generate data\n",
    "    \n",
    "    sigma_n = (n/sigma**2 + 1./tau_0**2)**(-0.5)\n",
    "    \n",
    "    ybar = y.mean()\n",
    "    mu_n = (n*ybar/sigma**2 + mu_0/tau_0**2)*sigma_n**2\n",
    "    \n",
    "    if (mu_true > mu_n - sigma_n) and (mu_true < mu_n + sigma_n):\n",
    "        n_enclosed += 1\n",
    "        \n",
    "print '%3.2f'%(n_enclosed/float(nsim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order for the posterior probability distribution to be significantly affected by the prior, we should:\n",
    "1. Pick a low number of measurements\n",
    "2. Assume a narrow prior\n",
    "\n",
    "For instance, let's set $n=4$ and $\\tau_0 = 0.3$ (not a good prior choice). Let's plot the prior and posterior distribution for one simulation, and compare the posterior with the one obtained assuming a flat prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pylab\n",
    "\n",
    "mu_0 = 0.\n",
    "n = 4\n",
    "tau_0 = 0.3\n",
    "\n",
    "mu_true = 2.*np.random.rand(1) - 1.\n",
    "y = np.random.normal(mu_true, sigma, n)\n",
    "\n",
    "mu_arr = np.linspace(-2., 2.)\n",
    "\n",
    "ybar = y.mean()\n",
    "sigma_n = (n/sigma**2 + 1./tau_0**2)**(-0.5)\n",
    "mu_n = (n*ybar/sigma**2 + mu_0/tau_0**2)*sigma_n**2\n",
    "\n",
    "mu_n_flat = ybar\n",
    "sigma_n_flat = sigma/n**0.5\n",
    "\n",
    "pylab.plot(mu_arr, 1./(2.*np.pi)**0.5/tau_0 * np.exp(-0.5*(mu_arr - mu_0)**2/tau_0**2), label='Prior')\n",
    "pylab.plot(mu_arr, 1./(2.*np.pi)**0.5/sigma_n * np.exp(-0.5*(mu_arr - mu_n)**2/sigma_n**2), label='Posterior (Gaussian prior)')\n",
    "pylab.plot(mu_arr, 1./(2.*np.pi)**0.5/sigma_n_flat * np.exp(-0.5*(mu_arr - mu_n_flat)**2/sigma_n_flat**2), label='Posterior (Flat prior)')\n",
    "pylab.axvline(mu_true, linestyle='--', color='k', label='$\\mu_{true}$')\n",
    "pylab.legend(loc='upper left')\n",
    "pylab.xlabel('$\\mu$')\n",
    "pylab.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_enclosed = 0\n",
    "\n",
    "for i in range(nsim):\n",
    "    mu_true = 2.*np.random.rand(1) -1. # draw a random value of the true mean of the distribution\n",
    "    y = np.random.normal(mu_true, sigma, n) # generate data\n",
    "    \n",
    "    sigma_n = (n/sigma**2 + 1./tau_0**2)**(-0.5)\n",
    "    \n",
    "    ybar = y.mean()\n",
    "    mu_n = (n*ybar/sigma**2 + mu_0/tau_0**2)*sigma_n**2\n",
    "    \n",
    "    if (mu_true > mu_n - sigma_n) and (mu_true < mu_n + sigma_n):\n",
    "        n_enclosed += 1\n",
    "        \n",
    "print '%3.2f'%(n_enclosed/float(nsim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With few data points and the poor prior choice made above, the fraction of times the inferred 68% region contains the true value of the mean of the distribution is much lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}