Winter School (Questions).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import subprocess\n",
    "\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor as GPR\n",
    "from sklearn.gaussian_process import kernels\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import from scipy import stats\n",
    "import emcee\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1 - Latin Hypercube Design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_lhs(npoints, ndim, seed):\n",
    "    \"\"\"\n",
    "    Generate a maximin Latin-hypercube sample (LHS) with the given number of\n",
    "    points, dimensions, and random seed.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    proc = subprocess.run(\n",
    "        ['R', '--slave'],\n",
    "        input=\"\"\"\n",
    "        library('lhs')\n",
    "        set.seed({})\n",
    "        write.table(maximinLHS({}, {}), col.names=FALSE, row.names=FALSE)\n",
    "        \"\"\".format(seed, npoints, ndim).encode(),\n",
    "        stdout=subprocess.PIPE,\n",
    "        check=True\n",
    "    )\n",
    "\n",
    "    lhs = np.array(\n",
    "        [l.split() for l in proc.stdout.splitlines()],\n",
    "        dtype=float\n",
    "    )\n",
    "\n",
    "\n",
    "    return lhs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create and plot design matrix below\n",
    "##### 20 points, 2 dimensions, set the seed to 80"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Fill in the parameters for generate_lhs(), and for plt.scatter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##Fill in the parameters to generate_lhs\n",
    "design = generate_lhs(npoints = ,\n",
    "                      ndim = ,\n",
    "                      seed = )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "###Fill in parameters to plt.scatter\n",
    "plt.scatter(x = ,\n",
    "           y = ) \n",
    "plt.title('Latin Hypercube Design')\n",
    "plt.xlabel('Input 1')\n",
    "plt.ylabel('Input 2')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 - Toy GP Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part a - Mean and Variance Estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the two cells below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def truth(x):\n",
    "    return(3*x+np.cos(5*x))    \n",
    "design = np.linspace(start =-1,stop=1,num=5)\n",
    "model_data = truth(design)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.scatter(x=design,y=model_data)\n",
    "plt.title('Computer Model Output at Design Points')\n",
    "plt.xlabel('Design')\n",
    "plt.ylabel('Model Output')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell below, training the GP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ptp = 2\n",
    "kernel = (\n",
    "    1. * kernels.RBF(\n",
    "        length_scale=ptp,\n",
    "        length_scale_bounds=np.outer(ptp, (.1, 10))\n",
    "    ) \n",
    ")\n",
    "gp = GPR(kernel=kernel,\n",
    "    #alpha=0,\n",
    "    n_restarts_optimizer=0,\n",
    "    copy_X_train=False).fit(design.reshape(-1,1), model_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Create the vector of points X on which we will predict. X should be a (n x 1) numpy array (you choose n)\n",
    "\n",
    "    Hint: We only want to predict within the bounds of our design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = \n",
    "X = X.reshape(-1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#This returns the predictive mean and covariance at all the points in X\n",
    "#mean is a (n,) numpy array, and cov is a (n,n) numpy array\n",
    "mean, cov = gp.predict(return_cov=True,X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Fill in the lines for plt.plot, top_var, and bot_var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Set up the figure by first plotting the output at the design points\n",
    "plt.scatter(x = design,y = model_data,color = 'black',label = 'Design Output')\n",
    "plt.title('Computer Model Output at Design Points')\n",
    "plt.xlabel('Input')\n",
    "plt.ylabel('Model Output')\n",
    "\n",
    "\n",
    "#Add the mean 95% uncertainty interval of the GP predictions at all the in-between points\n",
    "##Fill in parameters to plt.plot\n",
    "##Find expressions for top_var and bot_var\n",
    "###### Hint: For a normal distribution, how many standard deviations away from the mean is the 97.5% quantile?\n",
    "###### Hint: Use the diagonal of predictive covariance matrix`\n",
    "plt.plot( , ,color= 'blue',label = 'GP Mean')\n",
    "top_var =\n",
    "bot_var = \n",
    "plt.fill_between(X[:,0], bot_var, top_var, where=top_var >= bot_var, facecolor='lightgray', interpolate=True)\n",
    "\n",
    "plt.plot(X,truth(X),color='black',label = 'Truth')\n",
    "plt.legend(loc='best', fontsize=12)\n",
    "\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Part b - Random Draws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prior exercises displayed the mean and variance of our function at all the points in X. But what about actual samples of the function itself? Here we visualize what a random draw of the function would look like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercse**: Fill in the lines for top_var, bot_var, and rand_draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Hint: For rand_draw, examine np.random.multivariate_normal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Get the upper 95% quantile, and lower 95% quantile of the GP predictions at all the in-between points\n",
    "##Fill in same values as previous exercise\n",
    "top_var = \n",
    "bot_var = \n",
    "\n",
    "plt.scatter(design,model_data,color = 'black',label = 'Design Output')\n",
    "plt.title('Computer Model Output at Design Points')\n",
    "plt.xlabel('Input')\n",
    "plt.ylabel('Model Output')\n",
    "plt.fill_between(X[:,0], bot_var, top_var, where=top_var >= bot_var, facecolor='lightgray', interpolate=True)\n",
    "plt.plot(X,truth(X),color='black',label = 'Truth')\n",
    "\n",
    "ndraws = 10\n",
    "colors = cm.rainbow(np.linspace(0, 1, ndraws))\n",
    "\n",
    "#Get [ndraws] random draws from the predictive distribution of the GP at all of the in-between points \n",
    "#Use the predictive mean and covariance\n",
    "rand_draw = \n",
    "for i in range(ndraws):\n",
    "    plt.plot(X,rand_draw[i,:],color = colors[i],linestyle = \":\")\n",
    "\n",
    "plt.legend(loc='best', fontsize=14)\n",
    "\n",
    "pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3 - Principal Component Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset below contains 20 developmental indices for 132 countries around the world. The country labels (rows) can be found in 'countries.txt,' and the index descriptions can be found in 'indices.txt.' Here we look to explore some of the practical uses of PCA for the purposes of Computer Emulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "developmental_indices = np.loadtxt('dev_indices.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "scaler = StandardScaler(copy=False)\n",
    "pca = PCA(copy=False, whiten=True, svd_solver='full')\n",
    "Z = pca.fit_transform(scaler.fit_transform(developmental_indices))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Plot the cummulative fraction of variance explained. How many PCs would you recommend using?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Hint: Examine attributes of the object pca, as well as np.cumsum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F_r = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(range(len(F_r)),F_r,'-o')\n",
    "plt.title('Fraction of Variance Explained')\n",
    "plt.xlabel('Number of Components')\n",
    "plt.ylabel('F_r')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Find the correlation between the first two principal component vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "corr = \n",
    "corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Plot the second principal component against the first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.scatter(,)\n",
    "plt.title('Principle Components')\n",
    "plt.xlabel('Component 1')\n",
    "plt.ylabel('Component 2')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4 - Bayes Rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we're going to explore the mechanics of Bayesian analysis, and how prior choices can affect the posterior outcome. Let's revisit the coin example from the slides - we flipped 10 coins, and 7 came up heads. Our goal to make inference on the true probability of heads $\\theta$, using Bayesian analysis. To do this, we need a likelihood and prior. \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coin_tosses = np.loadtxt('coin_tosses.txt')\n",
    "num_heads = sum(coin_tosses)\n",
    "N = len(coin_tosses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Likelihood**:\n",
    "A common model for independent coin flips (or anything with a 0/1 outcome) is Binomial. So our likelihood is Binomial with size N = 10 and probability $\\theta$:\n",
    "\\begin{align}\n",
    "    y\\sim \\text{Binom}(N,\\theta)\n",
    "\\end{align}\n",
    "\n",
    "The mean of a Binomial random variable $Y$ with size $N$ (flips) and probability $\\theta$ (probability of heads) is\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathsf{E}(Y) = N\\theta.\n",
    "\\end{align}\n",
    "\n",
    "The variance of $Y$ is\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathsf{V}(Y) = N\\theta(1-\\theta)\n",
    "\\end{align}\n",
    "\n",
    "More information on Binomial distribution can found here: https://en.wikipedia.org/wiki/Binomial_distribution. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Prior**:\n",
    "The most widely used prior for Binomial data and unknown probability is a Beta distribution. A Beta random variable can take value in (0,1).  \n",
    "\n",
    "\\begin{align}\n",
    "    \\theta\\sim \\text{Beta}(a, b)\n",
    "\\end{align}\n",
    "\n",
    "The mean of a Beta random variable $\\theta$ with first parameter $a$ and second paramter $b$ is\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathsf{E}(\\theta) = \\frac{a}{a+b}\n",
    "\\end{align}\n",
    "\n",
    "The variance of $\\theta$ is\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathsf{V}(\\theta) = \\frac{a + b}{(a+b)^2(a+b+1)}\n",
    "\\end{align}\n",
    "\n",
    "More information on the Beta distribution can found here: https://en.wikipedia.org/wiki/Beta_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Goal: Explore different priors for $\\theta$ to see how the posterior of $\\theta$ responds. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 1: MCMC via python package emcee"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in missing pieces to calc_lnprior and calc_lnlike. These functions calculate the log of the prior and likelihood, respectively. For example, for the prior, you need to calculate the log of the pdf at theta for given parameters $a$ and $b$. Use the Wikipedia pages to find the pdfs.\n",
    "\n",
    "    Hint 1: For the Beta(a,b) prior, you will need to calculate the Beta function of parameters a and b (see the Wikipedia page). The cell below imports betln, which calculates the natural log of the Beta function.\n",
    "\n",
    "    Hint 2: For the Binomial likelihood, you'll need to calculate \n",
    "$$\\left(\\begin{matrix}N\\\\y\\end{matrix}\\right)$$\n",
    "    \n",
    "    This is the number of ways to get y heads in N flips. The function nCr(N,y) in the cell below calculates this.\n",
    "\n",
    "\n",
    "    Hint 3: If done correctly:\n",
    "\n",
    "    prior_ln_pdf(theta=0.4, a = 2,b = 3) = 0.5470\n",
    "\n",
    "    likelihood_ln_pdf(theta=0.5,y=num_heads,N=10) = -2.1440\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.special import betaln\n",
    "import math\n",
    "\n",
    "def nCr(n,r):\n",
    "    f = math.factorial\n",
    "    return f(n) // f(r) // f(n-r)\n",
    "\n",
    "print(nCr(4,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Calculate the log pdf of the prior of theta\n",
    "def prior_ln_pdf(theta,a,b):\n",
    "    #Ensure theta is within (0,1)\n",
    "    if theta <=0 or theta >=1:\n",
    "        return -np.inf\n",
    "    \n",
    "    ######\n",
    "    ##Add function return\n",
    "    ######\n",
    "    return \n",
    "\n",
    "#Calculate the log pdf of the likelihood of y successes (heads) in N trials (flips)\n",
    "def likelihood_ln_pdf(y, theta, N):\n",
    "    ######\n",
    "    ##Add function return\n",
    "    ######\n",
    "    return \n",
    "\n",
    "def posterior_ln_pdf(theta,a,b,y, N):\n",
    "    ln_pr = prior_ln_pdf(theta=theta,a=a,b=b)\n",
    "    \n",
    "    if not np.isfinite(ln_pr):\n",
    "        return -np.inf\n",
    "        \n",
    "    ln_like = likelihood_ln_pdf(y=y, theta=theta,  N = N)\n",
    "    return ln_pr + ln_like"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below sets up the sampler. We provide the number of walkers nwalkers (chosen), the number of parameters ndim (1 in our case, since we just have $\\theta$), and the posterior function which takes as argument the parameters we're inferrin on. \"args\" is a tuple of all of arguments to the probability function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = 1\n",
    "b = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nsamples = 100\n",
    "nburnin= 50\n",
    "ndim, nwalkers = 1, 200\n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, posterior_ln_pdf, args=(a, b, num_heads,N))\n",
    "p0 = np.random.rand(nwalkers,ndim)\n",
    "out_post = sampler.run_mcmc(p0,nsamples)\n",
    "samples = sampler.chain[:, nburnin:, :].reshape((-1, ndim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$nsteps$ is the number of steps each walker will run, while $nburnin$ is the number of \"burn-in\" or \"warmup\" steps for each walker.  Our total number of samples will be $nwalkers\\times (nsteps-nburnin)$ (recall we want 10,000 total samples). After running each walker for $nsteps$ and discarding the first $nburnin$, we reshape the output chain so each row is a draw from the posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nsamples = 100\n",
    "nburnin= 50\n",
    "out_post = sampler.run_mcmc(p0,nsamples)\n",
    "samples = sampler.chain[:, nburnin:, :].reshape((-1, ndim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the samples to find posterior means and variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Compare the prior density to the posterior samples\n",
    "def plot_comparison_mcmc(samples,prior_a, prior_b):\n",
    "    theta_grid = np.linspace(0.01,0.99,100)\n",
    "    \n",
    "    #Calculate the posterior density at the grid of theta values\n",
    "    prior_pdf_vals = stats.beta(a,b).pdf(theta_grid)\n",
    "    \n",
    "    \n",
    "    plt.plot(theta_grid,prior_pdf_vals,label = 'Prior Density' )\n",
    "    plt.hist(samples,normed = True,label = 'Posterior Samples')\n",
    "    plt.legend(loc='best', fontsize=12)\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: In the above code, change a and b to plot a comparison of the prior density to the histogram of posterior samples for each set of priors. Also report the prior mean versus the mean of the posterior samples\n",
    "\n",
    "a = 1, b = 1\n",
    "\n",
    "a = 1, b = 10\n",
    "\n",
    "a = 10, b = 1\n",
    "\n",
    "a = 30, b = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Method 2: Direct computation of the posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can find the analytics posterior in this case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, recall the likelihood and prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\ty\\mid\\theta&\\sim \\mathsf{Binom}(N,\\theta)\\\\\n",
    "\t\\theta&\\sim \\mathsf{Beta}(a,b)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we apply Bayes Rule and do a little algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\tp(\\theta\\mid y)&\\propto p(y\\mid\\theta)p(\\theta)\\\\\n",
    "\t &= \\left[\\begin{pmatrix}N \\\\ y\\end{pmatrix}\\theta^y(1-\\theta)^{N-y}\\right]\\frac{\\theta^{a-1}(1-\\theta)^{b-1}}{B(a,b)}\\\\\n",
    "\t &\\propto \\theta^{a+y-1}(1-\\theta)^{b+N-y-1}\n",
    "\\end{align}\n",
    "Note that when doing this algebra, we can disregard the constants $\\begin{pmatrix}N \\\\ y\\end{pmatrix}$ and $B(a,b)$ because they fall away in the proportionality.\n",
    "\n",
    "At the end, we recognize the kernel of a $\\mathsf{Beta}(a^*,b^*)$ distribution, where\n",
    "\\begin{align}\n",
    "\ta^* &= a + y\\\\\n",
    "\tb^* &= b + N - y\n",
    "\\end{align}\n",
    "This is posterior distribution of $\\theta$.\n",
    "\n",
    "Note that this math was very simple only because of the choice of prior and likelihood; a Beta distribution is a so-called \"conjugate prior\" for a Binomial likelihood. For more comlicated models the posterior will not be available analytically, and we must resort to sampling methods. We only used sampling methods in this model for pedagogical purposes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Find the parameters for the posterior distribution of $\\theta$ in our model. Compare the prior density to posterior density for the above choices of $a$ and $b$. Do the posterior densities align with the histogram of your posterior samples?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plot_comparison_analytic(prior_a,prior_b, posterior_a, posterior_b):\n",
    "\n",
    "    theta_grid = np.linspace(0.01,0.99,100)\n",
    "    \n",
    "    #Calculate the prior pdf at the grid of theta values\n",
    "    prior_pdf_vals = stats.beta(prior_a,prior_b).pdf(theta_grid)\n",
    "    plt.plot(theta_grid, prior_pdf_vals)\n",
    "    plt.xlabel('Theta')\n",
    "    plt.ylabel('Pdf')\n",
    "    \n",
    "    #Calculate the posterior pdf at the grid of theta values\n",
    "    posterior_pdf_vals = stats.beta(posterior_a,posterior_b).pdf(theta_grid)\n",
    "    plt.plot(theta_grid,posterior_pdf_vals, label = 'Posterior')\n",
    "    plt.legend(loc='best', fontsize=12)\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = \n",
    "b = \n",
    "N = len(coin_tosses)\n",
    "y = sum(coin_tosses)\n",
    "a_star = \n",
    "b_star = \n",
    "\n",
    "plot_comparison_analytic(a,b, a_star, b_star)"
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}