{
 "metadata": {
  "name": "",
  "signature": "sha256:6a55e9ceaba7fd56fcc161d21a4d3879deffe27fe83a1cd84bd00e91ed549a21"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from __future__ import division\n",
      "rcParams['figure.figsize'] = (10.0, 8.0)\n",
      "rcParams['font.size'] = 14\n",
      "\n",
      "from scipy import stats"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Posterior Sampling Example"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Consider the case of a sample $\\{x_i\\}$ of size $n$ drawn from a normal distribution of unknown mean $\\mu$ and standard deviation $\\sigma$ with the improper prior\n",
      "$$\n",
      "f(\\mu,\\sigma|I) \\propto \\frac{1}{\\sigma} \\qquad -\\infty<\\mu<\\infty,\\ 0<\\sigma<\\infty\n",
      "$$\n",
      "You showed on problem set 7 that the joint posterior pdf can be written\n",
      "$$\n",
      "f(\\mu,\\sigma|\\{x_i\\},I) \\propto \\frac{1}{\\sigma^{n+1}} \\exp\n",
      "\\left(\n",
      "-\\frac{1}{2\\sigma^2}\n",
      "\\left[\n",
      "(n-1)s^2+n(\\mu-\\bar{x})^2\n",
      "\\right]\n",
      "\\right)\n",
      "$$\n",
      "If we change variables from $\\sigma$ to $\\ln\\sigma$, this can be written as\n",
      "$$\n",
      "f(\\mu,\\ln\\sigma|\\{x_i\\},I) \\propto \\frac{1}{\\sigma^{n}} \\exp\n",
      "\\left(\n",
      "-\\frac{1}{2\\sigma^2}\n",
      "\\left[\n",
      "(n-1)s^2+n(\\mu-\\bar{x})^2\n",
      "\\right]\n",
      "\\right)\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'd like to draw a sample $\\{\\mu_s,\\ln\\sigma_s\\}$ from this distribution.  It turns out that $\\bar{x}$ and $s$ just act as location and scale parameters for this distribution, so the only thing that affects the shape is the sample size $n$.  We take $n=5$ for this illustration and suppose for simplicity that $\\bar{x}=0$ and $s=1$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n=5; xbar=0; s=1;"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will evaluate the pdf on a grid of $\\mu$ and $\\ln\\sigma$ values and use that to draw the sample."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Ngrid = 101\n",
      "mumin = -6\n",
      "mumax = 6\n",
      "lnsigmamin = -1\n",
      "lnsigmamax = 2\n",
      "muvals = linspace(mumin,mumax,Ngrid)\n",
      "lnsigmavals = linspace(lnsigmamin,lnsigmamax,Ngrid)\n",
      "mugrid,lnsigmagrid = meshgrid(muvals,lnsigmavals)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `meshgrid` function produces arrays which repeat the vectors of mu or lnsigma values as the rows or columns of the output:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(muvals)\n",
      "print(mugrid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(lnsigmavals)\n",
      "print(lnsigmagrid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sigmagrid = exp(lnsigmagrid)\n",
      "posteriorgrid0 = sigmagrid**(-n) * exp(-0.5*sigmagrid**(-2)*\n",
      "                                      ((n-1)*s**2+n*(mugrid-xbar)**2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(posteriorgrid0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We run the risk of overflow or underflow errors because of all the exponentials; it's better to work with the logarithms to start with, subtract off a constant value (since the proportionality constant is irrelevant) and then exponentiate:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lnpostgrid = -n*lnsigmagrid - 0.5*sigmagrid**(-2)*((n-1)*s**2+n*(mugrid-xbar)**2)\n",
      "lnpostgrid -= numpy.max(lnpostgrid)\n",
      "posteriorgrid = exp(lnpostgrid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the posterior is unnormalized, and in fact scaled so the pdf has a maximum value of 1."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print numpy.min(posteriorgrid/posteriorgrid0)\n",
      "print numpy.max(posteriorgrid/posteriorgrid0)\n",
      "print numpy.max(posteriorgrid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can make a contour plot of the gridded posterior (but note that we don't know what fraction of the volume under the posterior is outside the last curve)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "contour(mugrid,lnsigmagrid,posteriorgrid)\n",
      "xlabel(r'$\\mu$');\n",
      "ylabel(r'$\\ln\\sigma$');\n",
      "title(r'Contours of $f(\\mu,\\ln\\sigma|\\{x_i\\},I)$');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that since we have a pseudo-random number generator, we supply it with a seed so that the results are reproducible (but not predictable)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "random.seed(20171130)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To do a grid-approximated sample, we set the probability of drawing a grid point proportional to the posterior pdf."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "probs = posteriorgrid.flatten()\n",
      "probs /= sum(probs)\n",
      "Nsample = 10**4\n",
      "idxsample=random.choice(arange(Ngrid*Ngrid),size=Nsample,p=probs)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We've generated an array of indices corresponding to a \"flattened\" one-dimensional version of the array.  (I.e., we can turn each $101\\times101$ array into a $10201$-element list."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "idxsample"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we use this list of ten thousand indices to pick out ten thousand pairs of $\\mu$ and $\\ln\\sigma$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "musample=mugrid.flatten()[idxsample]\n",
      "lnsigmasample=lnsigmagrid.flatten()[idxsample]\n",
      "print(musample)\n",
      "print(lnsigmasample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we scatter-plot this sample to see what it looks like:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(musample,lnsigmasample,'k.');\n",
      "xlabel(r'$\\mu$');\n",
      "ylabel(r'$\\ln\\sigma$');\n",
      "xlim(mumin,mumax);\n",
      "ylim(lnsigmamin,lnsigmamax);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This looks a bit odd because $\\mu$ and $\\ln\\sigma$ were only allowd to take on specific discrete values.  We can more accurately simulate a continuous distribution by adding a random \"jitter\" that drops the points randomly in a box of $\\pm$ one-half grid spacing:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dmu = muvals[1]-muvals[0]\n",
      "dlnsigma = lnsigmavals[1]-lnsigmavals[0]\n",
      "musamplejit = musample + random.uniform(-0.5*dmu,0.5*dmu,Nsample)\n",
      "lnsigmasamplejit = lnsigmasample + random.uniform(-0.5*dlnsigma,0.5*dlnsigma,Nsample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(musamplejit,lnsigmasamplejit,'k.',ms=3);\n",
      "xlabel(r'$\\mu$');\n",
      "ylabel(r'$\\ln\\sigma$');\n",
      "xlim(mumin,mumax);\n",
      "ylim(lnsigmamin,lnsigmamax);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The density pf points in the sample roughly tracks the posterior pdf from which we've sampled."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Change of variables"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can convert this sample $\\{(\\mu_s,\\ln\\sigma_s)\\}$ into an sample $\\{(\\mu_s,\\sigma_s)\\}$ by just exponentiating the $\\ln\\sigma$ values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sigmasamplejit = exp(lnsigmasamplejit)\n",
      "sigmamin = exp(lnsigmamin)\n",
      "sigmamax = exp(lnsigmamax)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(musamplejit,sigmasamplejit,'k.',ms=3);\n",
      "xlabel(r'$\\mu$');\n",
      "ylabel(r'$\\sigma$');\n",
      "xlim(mumin,mumax);\n",
      "ylim(sigmamin,sigmamax);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that we don't need to calculate the Jacobian associated with the transformation, which we'd need to do to convert the explicit form of the pdf.  When we use sampling, we just need to transform each point, and the change in density takes care of itself."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Marginalization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To estimate the marginal posterior $f(\\mu|\\{x_i\\},I)=\\int_{-\\infty}^{\\infty} f(\\mu,\\ln\\sigma|\\{x_i\\},I)\\,d\\ln\\sigma$, we can just discard the $\\ln\\sigma$ values from the pairs and histogram the $\\{\\mu_s\\}$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hist(musamplejit,bins=sqrt(Nsample),normed=True);\n",
      "xlabel(r'$\\mu$');\n",
      "ylabel(r'$f(\\mu|\\{x_i\\},I)$')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We saw on the homework that this should be a Student $t$ distribution."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can do likewise for $\\ln\\sigma$ and $\\sigma$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hist(lnsigmasamplejit,bins=sqrt(Nsample),normed=True);\n",
      "xlabel(r'$\\ln\\sigma$');\n",
      "ylabel(r'$f(\\ln\\sigma|\\{x_i\\},I)$')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hist(sigmasamplejit**(-2),bins=sqrt(Nsample),normed=True);\n",
      "xlabel(r'$\\ln\\sigma$');\n",
      "ylabel(r'$f(\\ln\\sigma|\\{x_i\\},I)$')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Expectation values"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can estimate expectation values generated from the posterior by converting them to averages over the sample:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Posterior mean of mu is %g and standard deviation is %g'%(mean(musamplejit),std(musamplejit)))\n",
      "print('Posterior mean of ln(sigma) is %g and standard deviation is %g'%(mean(lnsigmasamplejit),std(lnsigmasamplejit)))\n",
      "print('Posterior mean of sigma is %g and standard deviation is %g'%(mean(sigmasamplejit),std(sigmasamplejit)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Percentiles"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Likewise, the percentiles of the posterior are approximated by percentiles of the sample:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print('Posterior median of mu is %g and 90pct PI is %g to %g'\n",
      "      %(median(musamplejit),percentile(musamplejit,5),percentile(musamplejit,95)))\n",
      "print('Posterior median of ln(sigma) is %g and 90pct PI is %g to %g'\n",
      "      %(median(lnsigmasamplejit),percentile(lnsigmasamplejit,5),percentile(lnsigmasamplejit,95)))\n",
      "print('Posterior median of sigma is %g and 90pct PI is %g to %g'\n",
      "      %(median(sigmasamplejit),percentile(sigmasamplejit,5),percentile(sigmasamplejit,95)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}