{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Eight Schools - An Edward/Stan Comparison\n", "\n", "In this tutorial we fit the 8 schools model to data made famous by __[Stan](https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started)__, __[Bayesian Data Analysis](http://www.stat.columbia.edu/~gelman/book/)__ and first appearing in __[Rubin (1981) - (paywalled)](http://journals.sagepub.com/doi/abs/10.3102/10769986006004377)__. The tutorial will be particularly suited to Stan users who are familiar with hierarchical models that would like to start using Edward.\n", "\n", "\n", "In this example we introduce a slightly modified 8 Schools example and we compare Stan using the No-U-Turn sampler (NUTS) and Edward using Automatic Differentiation Variational Inference (ADVI/KLQP) and Hamiltonian Monte Carlo (HMC). Both the Stan and Edward code are included in this tutorial and you may choose to install __[pystan](http://pystan.readthedocs.io/en/latest/getting_started.html)__ or __[rstan](https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started)__ to run the Stan portions of code. Comparison of Stan and Edward implementations can be a good way to check the algorithm and model is correct. As we use our own Stan implementation of 8 schools our results will differ slightly with the improper priors used in the Stan __[getting started](https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started)__ page and the Normal/Cauchy priors used in __[Diagnosing Biased Inference with Divergences](http://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html)__. \n", "\n", "We use the non-centered parameterization as advocated in __[Diagnosing Biased Inference with Divergences](http://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html)__ - which is highly recommended reading, just note all three models use slightly different priors when comparing results.\n", "\n", "In the centered (and simplified) form eight schools is a hierarchical model:\n", "\n", "$\\mu \\sim \\mathcal{N}(0,10)$\n", "\n", "$\\log \\tau \\sim \\mathcal{N}(5,1)$\n", "\n", "$\\theta \\sim \\mathcal{N}(\\mu, (e^{ \\log \\tau})^2)$\n", "\n", "$y_j \\sim \\mathcal{N}(\\theta, \\sigma^2_j)$\n", "\n", "where $j \\in [1..8]$ and $\\{ y_j, \\sigma_j \\}$ are observed.\n", "\n", "It may be easier to think of $\\tau$ as the root variance, which we give a log normal prior to.\n", "\n", "\n", "\n", "The non-centered model is equivalent and more suitable for inference and takes the following form:\n", "\n", "$\\mu \\sim \\mathcal{N}(0,10)$\n", "\n", "$\\log \\tau \\sim \\mathcal{N}(5,1)$\n", "\n", "$\\theta' \\sim \\mathcal{N}(0,1)$\n", "\n", "$y_j \\sim \\mathcal{N}(\\mu + (e^{ \\log \\tau}) \\theta', \\sigma_j^2)$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/d.rohde/tensorflow3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6\n", " return f(*args, **kwds)\n" ] } ], "source": [ "from __future__ import absolute_import\n", "from __future__ import division\n", "from __future__ import print_function\n", "import numpy as np\n", "import edward as ed\n", "import tensorflow as tf\n", "from edward.models import Normal, Empirical\n", "import matplotlib.pyplot as plt\n", "import pystan" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Data" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "J = 8\n", "data_y = np.array([28, 8, -3, 7, -1, 1, 18, 12])\n", "data_sigma = np.array([15, 10, 16, 11, 9, 11, 10, 18])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Edward Model Definition\n", "\n", "Let's now define the same model using Edward:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "mu = Normal(0., 10.)\n", "logtau = Normal(5., 1.)\n", "theta_prime = Normal(tf.zeros(J), tf.ones(J))\n", "sigma = tf.placeholder(tf.float32, J)\n", "y = Normal(mu + tf.exp(logtau) * theta_prime, sigma * tf.ones([J]))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "data = {y: data_y, sigma: data_sigma}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Edward Automatic Differentiation Variational Inference\n", "\n", "We now use ed.KLqp to perform ADVI. To do this we define the variational approximation family (as Gaussian) and then optimize the lower bound using Edwards ed.KLqp function." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "20000/20000 [100%] ██████████████████████████████ Elapsed: 56s | Loss: 37.511\n", "==== ed.KLqp inference ====\n", "E[mu] = 6.207960\n", "E[logtau] = 2.298177\n", "E[theta_prime]=\n", "[ 0.65755504 0.08147718 -0.25415799 0.0358108 -0.36784816 -0.22100581\n", " 0.55378842 0.13417724]\n", "==== end ed.KLqp inference ====\n" ] } ], "source": [ "with tf.variable_scope('q_logtau'):\n", " q_logtau = Normal(tf.get_variable('loc', []),\n", " tf.nn.softplus(tf.get_variable('scale', [])))\n", "\n", "with tf.variable_scope('q_mu'):\n", " q_mu = Normal(tf.get_variable('loc', []),\n", " tf.nn.softplus(tf.get_variable('scale', [])))\n", "\n", "with tf.variable_scope('q_theta_prime'):\n", " q_theta_prime = Normal(tf.get_variable('loc', [J]),\n", " tf.nn.softplus(tf.get_variable('scale', [J])))\n", "\n", "\n", "inference = ed.KLqp({logtau: q_logtau, mu: q_mu, theta_prime: q_theta_prime},\n", " data=data)\n", "inference.run(n_samples=15, n_iter=20000)\n", "\n", "print(\"==== ed.KLqp inference ====\")\n", "print(\"E[mu] = %f\" % (q_mu.mean().eval()))\n", "print(\"E[logtau] = %f\" % (q_logtau.mean().eval()))\n", "print(\"E[theta_prime]=\")\n", "print((q_theta_prime.mean().eval()))\n", "print(\"==== end ed.KLqp inference ====\")" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# compute 95% interval for theta'\n", "klqp_theta_prime_med = np.array(q_theta_prime.mean().eval())\n", "klqp_theta_prime_hi = np.array(q_theta_prime.mean().eval() + 1.96 *\n", " q_theta_prime.stddev().eval()\n", " )\n", "klqp_theta_prime_low = np.array(q_theta_prime.mean().eval() - 1.96 *\n", " q_theta_prime.stddev().eval()\n", " )\n", "\n", "# transform theta' to theta with S Monte Carlo samples\n", "S = 400000\n", "klqp_theta = np.tile(q_mu.sample(S).eval(), [8, 1]\n", " ).T + np.tile(np.exp(q_logtau.sample(S).eval()), [8, 1]\n", " ).T * q_theta_prime.sample(S).eval()\n", "# compute 95% interval for theta\n", "klqp_theta_low = np.array([np.percentile(klqp_theta[:, ii], 2.5)\n", " for ii in range(8)])\n", "klqp_theta_med = np.array([np.percentile(klqp_theta[:, ii], 50)\n", " for ii in range(8)])\n", "klqp_theta_hi = np.array([np.percentile(klqp_theta[:, ii], 97.5)\n", " for ii in range(8)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Edward HMC inference\n", "\n", "We now fit the same Edward model using HMC which is similar to Stan's NUTS algorithm but has more algorithm parameters to tune (although here we just use Edward's defaults).\n", "\n", "This is similar to ADVI except that the approximating family is Edward's Empirical type and the sampling is performed using Edwards ed.HMC function. Note: unlike Stan's self tuning NUTS algorithm HMC does need to be tuned for a particular problem. Here we find the defaults work acceptably, but this cannot be assumed in general." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100000/100000 [100%] ██████████████████████████████ Elapsed: 114s | Acceptance Rate: 0.914\n", "==== ed.HMC inference ====\n", "E[mu] = 5.225341\n", "E[logtau] = 2.454356\n", "E[theta_prime]=\n", "[ 0.66309202 0.11263941 -0.23707606 0.06023347 -0.31419596 -0.16834836\n", " 0.57439238 0.14914863]\n", "==== end ed.HMC inference ====\n" ] } ], "source": [ "S = 100000\n", "burn = S // 2\n", "hq_logtau = Empirical(tf.Variable(tf.zeros([S])))\n", "hq_mu = Empirical(tf.Variable(tf.zeros([S])))\n", "hq_theta_prime = Empirical(tf.Variable(tf.zeros([S, J])))\n", "\n", "inference = ed.HMC({logtau: hq_logtau, mu: hq_mu,\n", " theta_prime: hq_theta_prime}, data=data)\n", "inference.run()\n", "\n", "print(\"==== ed.HMC inference ====\")\n", "print(\"E[mu] = %f\" % (hq_mu.params.eval()[burn:].mean()))\n", "print(\"E[logtau] = %f\" % (hq_logtau.params.eval()[burn:].mean()))\n", "print(\"E[theta_prime]=\")\n", "print(hq_theta_prime.params.eval()[burn:, ].mean(0))\n", "print(\"==== end ed.HMC inference ====\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# compute 95% interval for theta'\n", "hmc_theta_prime_low = np.array([np.percentile(\n", " hq_theta_prime.params.eval()[burn:, ii], 2.5) for ii in range(8)])\n", "hmc_theta_prime_hi = np.array([np.percentile(\n", " hq_theta_prime.params.eval()[burn:, ii], 97.5) for ii in range(8)])\n", "hmc_theta_prime_med = np.array([np.percentile(\n", " hq_theta_prime.params.eval()[burn:, ii], 50) for ii in range(8)])\n", "\n", "# transform theta' to theta\n", "hmc_theta = np.tile(hq_mu.params.eval()[burn:], [8, 1]\n", " ).T + np.tile(np.exp(hq_logtau.params.eval()[burn:]), [8, 1]\n", " ).T * hq_theta_prime.params.eval()[burn:]\n", "\n", "# compute 95% interval for theta\n", "hmc_theta_low = np.array([np.percentile(hmc_theta[:, ii], 2.5)\n", " for ii in range(8)])\n", "hmc_theta_med = np.array([np.percentile(hmc_theta[:, ii], 50)\n", " for ii in range(8)])\n", "hmc_theta_hi = np.array([np.percentile(hmc_theta[:, ii], 97.5)\n", " for ii in range(8)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Stan Model Definition\n", "\n", "The Stan implementation of this non-centered model is:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "eight_schools_stan = \"\"\"\n", "data {\n", " int J;\n", " real y[J];\n", " real sigma[J];\n", "}\n", "\n", "parameters {\n", " real mu;\n", " real logtau;\n", " real theta_prime[J];\n", "}\n", "\n", "model {\n", " mu ~ normal(0, 10);\n", " logtau ~ normal(5, 1);\n", " theta_prime ~ normal(0, 1);\n", " for (j in 1:J) {\n", " y[j] ~ normal(mu + exp(logtau) * theta_prime[j], sigma[j]);\n", " }\n", "}\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Stan code quite closely follows the model definition. Note that putting the Stan code in a string like this is convenient for this tutorial, but usually you should instead put the code in a separate file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Stan NUTS inference \n", "\n", "This allows us to run the Stan code (make sure pystan is installed for this step)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_d553225c2fc13d05db15363950aa04c0 NOW.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "==== Stan/NUTS inference ====\n", "E[mu] = 5.836806\n", "E[logtau] = 2.450053\n", "E[theta_prime]=\n", "[ 0.64940756 0.09001582 -0.23279844 0.04471902 -0.33542507 -0.2041105\n", " 0.53249937 0.14456798]\n", "==== end Stan/NUTS inference ====\n" ] } ], "source": [ "standata = dict(J=J, y=data_y, sigma=data_sigma)\n", "\n", "fit = pystan.stan(model_code=eight_schools_stan,\n", " data=standata, iter=20000)\n", "\n", "s = fit.extract()\n", "print(\"==== Stan/NUTS inference ====\")\n", "print(\"E[mu] = %f\" % (s['mu'].mean()))\n", "print(\"E[logtau] = %f\" % (s['logtau'].mean()))\n", "print(\"E[theta_prime]=\")\n", "print(s['theta_prime'].mean(0))\n", "print(\"==== end Stan/NUTS inference ====\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# compute 95% interval for theta'\n", "stan_theta_prime_low = np.array([np.percentile(s['theta_prime'][:, ii], 2.5)\n", " for ii in range(8)])\n", "stan_theta_prime_hi = np.array([np.percentile(s['theta_prime'][:, ii], 97.5)\n", " for ii in range(8)])\n", "stan_theta_prime_med = np.array([np.percentile(s['theta_prime'][:, ii], 50)\n", " for ii in range(8)])\n", "# transform theta' to theta\n", "stan_theta = np.tile(s['mu'],\n", " [8, 1]\n", " ) + np.tile(np.exp(s['logtau']), [8, 1]\n", " ) * s['theta_prime'].T\n", "# compute 95% interval for theta\n", "stan_theta_low = np.array([np.percentile(stan_theta[ii, ],\n", " 2.5) for ii in range(8)])\n", "stan_theta_hi = np.array([np.percentile(stan_theta[ii, ],\n", " 97.5) for ii in range(8)])\n", "stan_theta_med = np.array([np.percentile(stan_theta[ii, ],\n", " 50) for ii in range(8)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Comparison of the approximate posteriors for $\\theta'$\n", "\n", "\n", "We can graphically show the the three methods have very similar $95\\%$ intervals for $\\theta'$:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAF3CAYAAADzW6zjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XuUVdWZ7/3vQ4FA0AYVRKRAoUWIEYSy1FcFmjRGIRiCdgKaYafhfY9kaAwkJ/EMPUkrJyMO7ZaOTVpzYSRBc2nBlyZEBC/EaEBQAwWIooJG7VglKmggmoBAMc8fe1NSUFVA3fZexfczxh7sNdfaaz+1vNSPudacM1JKSJIkKRvaFboASZIkHT7DmyRJUoYY3iRJkjLE8CZJkpQhhjdJkqQMMbxJkiRliOFNkiQpQwxvkiRJGWJ4kyRJyhDDmyRJUoa0L3QBLal79+7ptNNOK3QZkiRJh1RRUbE1pdTjUMe16fB22mmnsXr16kKXIUmSdEgR8d+Hc5y3TSVJkjLE8CZJkpQhhjdJkqQMadPPvEmSdDTbvXs3lZWV7Ny5s9ClaD+dOnWitLSUDh06NOrzhjdJktqoyspKjjvuOE477TQiotDlCEgp8e6771JZWUm/fv0adQ5vm0qS1Ebt3LmTE0880eBWRCKCE088sUm9oYY3SZLasCMNbpN+9BSTfvRUC1UjOPJ/JgcyvEmSpBZTUlLC0KFDa1633377Qcc88cQTXHbZZS1ax4wZM5g5c2bN9tNPP80111xT53dPnjyZ+fPnAzBq1Cj69u1LSqlm/4QJEzj22GNrtjdt2sSnP/1pBgwYQFlZGRMnTuTtt99usZ/FZ94kSRIAC9dWsfaP29hVvZeLbv8tN1w6kAnDejfpnJ07d2bdunXNVOHh2bNnD+3bNxxxHnroIcaMGXNY5+vWrRsrVqxg+PDhbNu2jc2bN9fs27lzJ+PGjeO73/0un/nMZ4BcGN2yZQs9e/Zs/A/RAHveJEkSC9dWcdOC59hVvReAqm07uGnBcyxcW9Ui3/fwww8zaNAgysrKWLBgQU374MGD2bZtGyklTjzxRH72s58B8MUvfpGlS5fy+uuvM2LECMrKyigrK2PlypVALjCNGDGC8ePHc+aZZwJw6623csYZZzB8+HA2btxY6/sfe+wxLr744sOq9corr2Tu3LkALFiwgCuuuKJm33/+539ywQUX1AQ3yPXWnXXWWY24KofH8CZJkrjjkY3s2F1dq23H7mrueGRjPZ84PDt27Kh123TevHns3LmTa665hkWLFlFRUcFbb71Vc/xFF13EihUr2LBhA/3792f58uUAPPXUU1x44YWcdNJJLF26lDVr1jBv3jymTZtW89k1a9Ywa9YsNm3aREVFBXPnzmXdunUsWbKEVatW1Ry3detWOnToQNeuXQFYvnx5rRofeOCBWj/D6NGjWbZsGdXV1cydO5dJkybV7Hv++ec555xzmnSNjpS3TSVJEm9u23FE7Yerrtum69ato1+/fgwYMACAq6++mtmzZwMwYsQIli1bxqmnnsq1117L7Nmzqaqq4vjjj6dLly5s376d66+/nnXr1lFSUsKmTZtqznveeefVTL+xfPlyLr/8cj72sY8BMH78+JrjHn30US655JKa7REjRvDggw/WbE+ePLlWvSUlJQwfPpy5c+eyY8cOTjvttCZdk6ay562lzBmXe0mSlAGndOt8RO0tZeTIkSxfvpzly5czatQoevTowfz58xkxYgQAd955Jz179uTZZ59l9erV7Nq1q+azXbp0OazvOJLn3fa58sormTZtGhMnTqzV/olPfIKKioojOldTGd6awOHUkqS24oZLB9K5Q0mtts4dSrjh0oHN/l2DBg3i9ddf5w9/+AMA9913X82+Pn36sHXrVl5++WX69+/P8OHDmTlzJiNHjgRg+/bt9OrVi3bt2vHzn/+c6urqOr9j5MiRLFy4kB07dvD++++zaNEiIDdJ7vr16xk6dOgR1TxixAhuuukmrrrqqlrtX/jCF1i5ciWLFy+uaVu2bBnPP//8EZ3/SBjeCmDKw1OY8vCUQpchSVKNCcN6c9sVgzmmJBcNenfrzG1XDG7yaNMDn3m78cYb6dSpE7Nnz2bcuHGUlZVx0kkn1frM+eefzxlnnAHkQlNVVRXDhw8H4LrrruPee+/l7LPP5qWXXqq3t62srIxJkyZx9tlnM3bsWM4991wAKioqGDZs2BHPtRYRfOMb36B79+612jt37syDDz7If/zHfzBgwADOPPNMvv/979OjR48jOv8R1bL/vCVtTXl5eVq9enWLnX9fr9u8L11w8M59t0ynLD5o177gNmfMnBarTZKkF198kY9//ONH9JkGf7e1Ad/5znc4/fTTufLKKwtaR13/bCKiIqVUfqjPOmBBkiTVaKuhbZ9vfetbhS6hybxtKkmSlCGGN0mSpAwxvEmSJGWI4U2SJClDDG+SJOkjTjJf9AxvkiSpxZSUlNSa5+32228/6JgnnniCyy67rEXrmDFjBjNnzqzZfvrpp7nmmmsO+u5vfetbjBkzhg8//JBRo0ZR15RjTz75JOeddx6DBg1i4MCBfP/736/1Pb1792bo0KGcddZZB62T2hycKkSSJOWsvx8qV0H1h3DnWTD6Zhgy8dCfa0Bda5u2tD179tC+fcMRp64lsr7zne+wYsUKlixZQseOHev83FtvvcUXvvAFFi5cSFlZGVu3buXSSy+lV69eXH755QB87Wtf4xvf+AYvvvgiI0aM4J133qFdu+brL7PnTZIk5YLbomm54Aaw/Y3c9vr7W+TrHn74YQYNGkRZWRkLFiyoaR88eDDbtm0jpcSJJ57Iz372MwC++MUvsnTpUl5//XVGjBhBWVkZZWVlrFy5Esj13o0YMYLx48dz5plnAnDrrbdyxhlnMHz4cDZu3Fjr+x977DEuvvjimu1/+7d/46GHHmLRokV07lz/eq533303kydPpqysDIDu3bvzr//6r9xxxx0HHfvxj3+c9u3bs3Xr1kZepboVPLxFRJ+IeDwiXoiIDRExvY5jIiK+FxGvRMT6iCgrRK2SJLVZj30bdu+o3bZ7R669CQ5cHmvevHns3LmTa665hkWLFlFRUcFbb71Vc/xFF13EihUr2LBhA/3792f58uUAPPXUU1x44YWcdNJJLF26lDVr1jBv3jymTZtW89k1a9Ywa9YsNm3aREVFBXPnzmXdunUsWbKEVatW1Ry3detWOnToQNeuXQFYsWIFP/zhD3nooYc49thjG/x5NmzYwDnnnFOrrby8nBdeeOGgY5955hnatWvX7EtlFcNt0z3A11NKayLiOKAiIpamlPa/CmOBAfnX+cAP8n9KkqTmsL3yyNoPU123TdetW0e/fv0YMGAAAFdffTWzZ88GcmuZLlu2jFNPPZVrr72W2bNnU1VVxfHHH0+XLl3Yvn07119/PevWraOkpIRNmzbVnPe8886jX79+ACxfvpzLL7+cj33sYwCMHz++5rhHH32USy65pGb79NNP509/+hNLly7lH/7hH5r08wLceeed/OIXv+C4445j3rx5R7yO6qEUvOctpbQ5pbQm//594EXgwFVwPwv8LOU8DXSLiF6tXKokSW1X19Ija28hI0eOZPny5SxfvpxRo0bRo0cP5s+fz4gRI4BcMOrZsyfPPvssq1evZteuXTWfrW+R+gMd+Lxbz549WbJkCV/96ld5/PHHG/zsmWeeSUVFRa22iooKyss/WpL0a1/7GuvWrWP58uU1dTengoe3/UXEacAw4JkDdvUG3thvu5KDA54kSWqs0TdDhwOe9erQOdfezAYNGsTrr7/OH/7wBwDuu+++mn19+vRh69atvPzyy/Tv35/hw4czc+ZMRo4cCcD27dvp1asX7dq14+c//znV1dV1fsfIkSNZuHAhO3bs4P3332fRokUApJRYv349Q4cOrXX8GWecwYIFC7j66qsbHGDx5S9/mXvuuafmmHfffZdvfvOb/PM//3PjL8gRKprwFhHHAv8FfDWl9OcmnGdqRKyOiNVbtmxpvgIlSWrLhkyEz3wPSvKjLLv2yW03cbTpgc+83XjjjXTq1InZs2czbtw4ysrKOOmkk2p95vzzz+eMM84AcrdRq6qqGD58OADXXXcd9957L2effTYvvfRSvb1tZWVlTJo0ibPPPpuxY8dy7rnnArlesmHDhtV5K/Pcc89lzpw5jB8/viZYjhs3jtLSUkpLS/n85z9Pr169+MUvfsHUqVMZOHAgp5xyCtOmTePv/u7vmnSdjkSklFrty+otIqID8CDwSErpu3Xs/xHwRErpvvz2RmBUSmlzQ+ctLy9Pdc3P0lwm/egpAOZ96YKDd+6b4HDK4oN2TXl4Su6QMXNarDZJkl588UU+/vGPH9mHGvj91RZ85zvf4fTTT+fKK69slvN9//vf5wc/+AHLli3j+OOPP+zP1fXPJiIqUkrl9XykRsEHLEQu+v4EeLGu4Jb3AHB9RMwlN1Bh+6GCmyRJaoQ2Gtr2+da3vtWs57vuuuu47rrrmvWch1Lw8AZcBPwj8FxE7LvJ/L+BvgAppR8CS4BPA68AfwWmFKBOSZKKXoN3hRrw2vbXAOjXtV+z16TmVfDwllJ6EmhwDG3K3dv9cutUJEmSVLyKZsCCJEmSDs3wJkmSlCGGN0mSVGPKw1NqZkVQcTK8SZKkFnXrrbfyiU98giFDhjB06FCeeeYZ/v3f/52//vWvTTrv7t27axaIjwi+/vWv1+ybOXMmM2bMAGDy5MnMnz+/1mePPfZYnnvuuZr550444QT69evH0KFDufjii9m7dy/Tpk3jrLPOYvDgwZx77rm89tprTaq3uRR8wIIkSSoOT7zxBOu3rGfX3l1cMv8SppdNZ1z/cU0651NPPcWDDz7ImjVr6NixI1u3bmXXrl1MmjSJq6++umbt0cZ48sknueiiiwDo2LEjCxYs4KabbqJ79+6H9fnBgwfXrJQwefJkLrvsMj73uc8BuVUf3nzzTdavX0+7du2orKw87OW3Wpo9b5IkHS22vpx71eGJN57g7nV3s2tvbq3QzX/ZzIyVM1j8atPmfdu8eTPdu3enY8fcyg3du3dn/vz5vPnmm3zyk5/kk5/8JADXXnst5eXlfOITn+CWW26p+fxpp53GLbfcQllZGYMHD+all16q2ffwww8zduxYANq3b8/UqVO58847m1Tv/nXvW4YLoLS09Igm4W1JhjdJksTPX/g5H1Z/WKttZ/VOZq2Z1aTzXnLJJbzxxhucccYZXHfddfzud79j2rRpnHLKKTz++OM1C8HfeuutrF69mvXr1/O73/2O9evX15yje/furFmzhmuvvZaZM2fWtD/++OOMGjWqZvvLX/4yv/zlL9m+fXuTagaYOHEiixYtYujQoXz9619n7dq1TT5nczG8SZIktu7YWmf7W395q0nnPfbYY6moqGD27Nn06NGDSZMmcc899xx03P33309ZWRnDhg1jw4YNvPDCCzX7rrjiCgDOOeccXn/9dQCqqqo44YQTat12/Zu/+Ru++MUv8r3vfa/Wuetax7Sutv2VlpayceNGbrvtNtq1a8fo0aN57LHHDvfHblE+8yZJkujeuTtbdmw5qP3kLic3+dwlJSWMGjWKUaNGMXjwYO69995a+1977TVmzpzJqlWrOP7445k8eTI7d+6s2b/vlmtJSQl79uwBcrdML7300oO+66tf/SplZWVMmfLRiNkTTzyRP/3pTzXb77333mE9F9exY0fGjh3L2LFj6dmzJwsXLmT06NFH9sO3AHveJEkS/3jmP9KxpGOttk4lnZheNr1J5924cSMvv/zRc3br1q3j1FNP5bjjjuP9998H4M9//jNdunSha9euvP322zz00EOHPO/+z7vt74QTTmDixIn85Cc/qWkbNWoU8+bNY9eu3PN899xzT82zdkDuOcCdf651njVr1vDmm28CsHfvXtavX8+pp556+D94C7LnTZIkMarPKADuWnsXu/buoleXXs0y2vSDDz7gK1/5Ctu2baN9+/acfvrpzJ49m/vuu48xY8bUPPs2bNgwBg0aRJ8+fWpGkNanurqaV155hUGDBtW5/+tf/zp33XVXzfZll11GRUUF55xzDiUlJfzt3/4tP/zhDxv8jnfeeYdrrrmGDz/MPQd43nnncf311x/hT98yDG+SJAnIBbhllcsAmDNmTrOc85xzzmHlypUHtX/lK1/hK1/5Ss12Xc/BATXPuAGUl5fzxBNP8OSTT3L++efXOu6DDz6oed+zZ8+D5pC75ZZbao1iPdA9d/0LdB9Qsz1mzBjGjBnDa9tzc7v169qv3s+2NsObJEmq0VyhrSUNHz6c4cOHF7qMgvGZN0mS1Ob9YcsH/GHLB4c+MAMMb5IkSRlieJMkqQ1LKRW6BB2gqf9MDG+SJLVRnTp14t133zXAFZGUEu+++y6dOnVq9DkcsCBJUhtVWlpKZWUlW7bkJ9/94J3cn1v2HHTsvhUWdnbeedC+tmDL+7kpP3Zt7Xjwzla+Lp06daK0tLTRnze8SZLURnXo0IF+/fab4mLON3J/Tjl4sfkpD+dWJMjCaNPGmPGjpwCY96WhB+/M2HXxtqkkSVKGGN4kSZIyxPAmSZKUIYY3SZKkDDG8SZIkZYjhTZIkKUMMb5IkSRlieJMkScqQoghvEfHTiHgnIp6vZ/+oiNgeEevyr5tbu0ZJkqRiUCwrLNwD3AX8rIFjlqeULmudciRJkopTUfS8pZSWAe8Vug5JkqRiVxTh7TBdEBHPRsRDEfGJQhcjSZJUCMVy2/RQ1gCnppQ+iIhPAwuBAXUdGBFTgakAffv2bb0KVWNSzeK/FxzR54px8V+1jsb+OyNJR6NM9LyllP6cUvog/34J0CEiutdz7OyUUnlKqbxHjx6tWqek1jXl4Sk1oV+SjhaZCG8RcXJERP79eeTqfrewVUmSJLW+orhtGhH3AaOA7hFRCdwCdABIKf0Q+BxwbUTsAXYAV6aUUoHKlSRJKpiiCG8ppasOsf8uclOJSJIkHdUycdtUkiRJOYY3SZKkDDG8SSpuc8blXpIkwPCm1uYv4kZxSgxJ0j6GN0mSpAwxvEmSJGWI4U0qFt5SliQdBsOb1Iom/eipmnU8JUlqDMObJElShhjeJEmSMsTwJkmSlCGGN0mSpAwxvEmSJGWI4U2SJClDDG+SpMxx2h0dzQxvkqSjhusEqy0wvEmSJGWI4U2SJClDDG+SJEkZYniTJEnKEMObJElShhjeJElty5xxuZfURhneJEmSMsTwJkmSlCGGN0mSpAwxvEmSJGWI4U2SJClDiiK8RcRPI+KdiHi+nv0REd+LiFciYn1ElLV2jZIkScWgKMIbcA8wpoH9Y4EB+ddU4AetUJMkSVLRKYrwllJaBrzXwCGfBX6Wcp4GukVEr9apTpIkqXgURXg7DL2BN/bbrsy3SZLqMOXhKUx5eEqhy5DUArIS3g5bREyNiNURsXrLli2FLkeSJKlZZSW8VQF99tsuzbcdJKU0O6VUnlIq79GjR6sUJ0ktZdKPnmLSj54qdBmSikhWwtsDwBfzo07/H2B7SmlzoYuSJElqbe0LXQBARNwHjAK6R0QlcAvQASCl9ENgCfBp4BXgr4APckiSpKNSUYS3lNJVh9ifgC+3UjmSJElFKyu3TSVJkoThTZIkKVMMb5IkSRlieJMkScoQw5skSVKGGN4kSZIyxPAmSVk1Z1zuJemoYniTJEnKEMObJElShhjeJEmSMsTwJkmSlCGGN0mSpAwxvEmSJGWI4U2SJClDDG+SJEkZYniTJEnKEMNbIy1cW8XaP27jmdfe46Lbf8vCtVWFLkmSJB0FDG+NsHBtFTcteI5d1XsBqNq2g5sWPPdRgFt/P1Sugv9+Eu48K7edt/jVxazfsp7Vb6/mkvmXsPjVxYX4ESRJUka1L3QBWXTHIxvZsbu6VtuO3dXc8chGJpSsgEXToPrD3I7tb+S2gcXHdmHGyhns2rsLgM1/2cyMlTMAGNff9QklSdKh2fPWCG9u21F/+2Pfht0H7N+da5+1ZhY7q3fW2rWzeiez1sxqqVIlSVIbY3hrhFO6da6/fXtl3R/aXslbf3mrzl31tUuSJB3I8NYIN1w6kM4dSmq1de5Qwg2XDoSupXV/qGspJ3c5uc5d9bVLkiQdyPDWCBOG9ea2KwZzTEnu8vXu1pnbrhjMhGG9YfTN0OGAnrkOnWH0zUwvm06nkk61dnUq6cT0sumtVbokSco4Byw00oRhvbnv938EYN6XLvhox5CJuT9/fX1u0ELXPrlAN2Qi+4Yk3LziZnbt3UWvLr2YXjbdwQqSJOmwGd5awpCJUHFv7v2U2lOBjOs/jvmb5gMwZ8yc1q5MktSG7ZuDdFf1Xi66/bfccOnA3F0htSneNpUkqQ1oyhykyhbDmyRJbUBDc5Cy/v665yA1wGWS4U2SpDagsXOQgqv/ZE1RhLeIGBMRGyPilYi4sY79kyNiS0Ssy7/+RyHqlCSpWDV2DtLFry6uc/UfA1zxKnh4i4gS4G5gLHAmcFVEnFnHofNSSkPzrx+3apGSJBW5xs5B6uo/2VPw8AacB7ySUno1pbQLmAt8tsA1SZKUKY2dg9TVf7KnGMJbb+CN/bYr820H+oeIWB8R8yOiT+uUJklSdkwY1pthfbtxfr8TWHHj3380TciQifCZ70FJx9x21z657SETXf0ng4ohvB2ORcBpKaUhwFLg3voOjIipEbE6IlZv2bKl1QqUJKmoDZkIpefCqcPha8/XTCrv6j/ZUwzhrQrYvyetNN9WI6X0bkopP76ZHwPn1HeylNLslFJ5Sqm8R48ezV6sJEltybj+45hx4QyOaXcMAL269GLGhTNc/aeIFcMKC6uAARHRj1xouxL4wv4HRESvlNLm/OZ44MXWLVGSpLbL1X+ypeDhLaW0JyKuBx4BSoCfppQ2RMS3gdUppQeAaRExHtgDvAdMLljBkiRJBVTw8AaQUloCLDmg7eb93t8E3NTadUmSJBWbYnjmTZIkSYfJ8CZJkpQhhjdJkqQMMbxJkiRliOFNkiQpQwxvkiRJGWJ4kyRJyhDDmyRJUoYY3iRJkjLE8CZJkpQhhjdJkqQMMbxJkiRliOFNkiQpQwxvkiRJGWJ4kyRJypAGw1tETIiIuyNiYGsVJEmSpPo1GN5SSguB7wKfjIiprVOSJElS81m4toq1f9zGM6+9x0W3/5aFa6sKXVKTHM5t02rgb4GxEfHjiLg+Ik5t4bokSZKabOHaKm5a8By7qvcCULVtBzcteO6wAtziVxezfst6Vr+9mkvmX8LiVxe3dLmH5XDC26+Bl4C7gE8BZwPL8rdTO7ZkcZIkSU1xxyMb2bG7ulbbjt3V3PHIxtzG+vuhchX895Nw51m5bXLBbcbKGezauwuAzX/ZzIyVM4oiwB1OeCtJKf0kpfQY8F5K6RpyPXGvA7NbsjhlT1vrmpYkZdub23bU377+flg0Dao/zDVufyO3vf5+Zq2Zxc7qnbU+s7N6J7PWzGrpkg/pcMLbbyLi+vz7BJBS2pNSugO4oMUqU+Y0pWtakqSWcEq3zvW3P/Zt2H1AuNu9Ax77Nm/95a06P1dfe2s6nPD2P4GuEbEaOCUipkbE1RFxN/Buy5anLDlk13QDivW5AklStt1w6UA6dyip1da5Qwk3XDoQtlfW/aHtlZzc5eQ6d9XX3poOGd5SSntTSrcCI4GpwMnAOcDzwNiWLU9Z0mDXNGTyuQJJUrZNGNab264YzDElucjTu1tnbrtiMBOG9YaupXV/qGsp08um06mkU63mTiWdmF42vaVLPqT2h3tgSumvwAP5l3SQU7p1pqqOAHdKt871P1cAzNr043qfKxjXf1yL163C2vec5K7qvVx0+2+54dKBuf+pSlIzmTCsN/f9/o8AzPvSfk98jb4597to/1unHTrD6Jtrfv/cvOJmdu3dRa8uvZheNr0ofi+5woKaTYNd0xl9rkAtqy0O4ZeUIUMmwme+ByX5yTO69sltD5kIwLj+4xjSYwjlPct59HOPFkVwA8ObmlGDXdMZfa5ALastDuGXlDFDJkLpuXDqcPja8zXBrZgZ3tSsJgzrzbC+3Ti/3wmsuPHvP7r9ldHnCtSy2uIQfklqaUUR3iJiTERsjIhXIuLGOvZ3jIh5+f3PRMRprV+lmmT0zbnnCPa333MFMy6cwTHtjgGgV5dezLhwRtF0T6vltMUh/JLU0goe3iKiBLib3MjVM4GrIuLMAw77/4A/pZROB+4E/qV1q1STZfS5ArWstjiEX5JaWsHDG3Ae8EpK6dWU0i5gLvDZA475LHBv/v18YHRERCvWqOaQwecK1LLa4hB+SWppxRDeegNv7LddmW+r85iU0h5gO3Biq1QnqUXV+5ykt9olqU6HPc9bVkTEVHKTCdO3b98CVyOp0fb1zP76+tygha59coFuv1vt8zfNB2DOmDmFqlKSWl0xhLcqoM9+26X5trqOqYyI9kBX6lmaK6U0G5gNUF5enpq9WkmtZ8hEqMg/MTHFaUAkCYrjtukqYEBE9IuIY4ArOXgVhweAf8q//xzw25SSwUySJB11Ct7zllLaExHXA48AJcBPU0obIuLbwOqU0gPAT4CfR8QrwHvkAp4kSdJRp+DhDSCltARYckDbzfu93wl8vrXrkiRJKjbFcNtUkiRJh8nwJkmSlCGGN0mSpAwxvEmSMmXh2irW/nEbz7z2Hhfd/lsWrj1wdimpbTO8SZIyY+HaKm5a8By7qvcCULVtBzcteM4Ap6OK4U2SlBl3PLKRHbura7Xt2F3NHY9sLFBFUuszvEmSMuPNbTuOqH1/i19dzPot61n99moumX8Ji1911Q5lk+FNkpQZp3TrfETt+yx+dTEzVs5g195dAGz+y2ZmrJxhgFMmGd4kSZlxw6UD6dyhpFZb5w4l3HDpwNzG+vuhchX895Nw51m5bWDWmlnsrN5Z63M7q3cya82sVqlbak6GN6mVNHaEnLd6pI9MGNab264YzDEluV9fvbt15rYrBjNhWO9cUFs0Dao/zB28/Y3c9vr7eesvb9V5vvrapWJmeJNawSFHyNX5XAQLAAAO9UlEQVTTW+CtHulgE4b1Zljfbpzf7wRW3Pj3ueAG8Ni3YfcBz77t3gGPfZuTu5xc57nqa5eKmeFNagUNjpBroLfAWz3SEdheWW/79LLpdCrpVKu5U0knppdNb4XCpOZVFAvTS21dgyPkGugteOuEuv9+5a0eqQ5dS3N/+amjfVz/cQDcvOJmdu3dRa8uvZheNr2mXcoSe96kVtDgCLkGegu81XN0cyWBIzT6ZuhwwH9rHTrn2oFx/ccxpMcQynuW8+jnHjW4KbMMb1IraHCEXNfSuj/UtdRbPUcxVxJohCET4TPfg5KOue2ufXLbQyYWti6pmRnepFbQ4Ai5BnoLxvUfx4wLZ3BMu2MA6NWlFzMunGGPwVGgKSsJHNUjlIdMhNJz4dTh8LXnDW5qk3zmTWolE4b15r7f/xGAeV+64KMd+365/Pr63KCFrn1ygS7fPq7/OOZvmg/AnDFzWrVmFU5jVxKob4QyYOiX2gh73qRiYG+BDtDYlQQcoSy1fYY3SSpCh1xJoB5ORiu1fYY3SSpCDT4nCfVO7OwIZantM7xJUpGqdyWBBiZ2doSy1PYZ3iQpaxqY2NkRylLb52hTScqaBiZ2BkcoS22dPW+SlDUNTOwsqe0zvElS1hxiGShJbZvhTZKyxmWgpKOaz7xJUhYNmQgV9+beTzmKlr+SVNiet4g4ISKWRsTL+T+Pr+e46ohYl3890Np1SpIkFYtC3za9EXgspTQAeCy/XZcdKaWh+df41itPkiSpuBQ6vH0WyPf7cy8woYC1SJIkFb1Ch7eeKaXN+fdvAT3rOa5TRKyOiKcjwoAnSZKOWi0+YCEifgPUtajeN/ffSCmliEj1nObUlFJVRPQHfhsRz6WU/lDP900FpgL07du3CZVLkiQVnxYPbymli+vbFxFvR0SvlNLmiOgFvFPPOaryf74aEU8Aw4A6w1tKaTYwG6C8vLy+MChJkpRJhb5t+gDwT/n3/wT8+sADIuL4iOiYf98duAh4odUqlCRJKiKFDm+3A5+KiJeBi/PbRER5RPw4f8zHgdUR8SzwOHB7SsnwJkmSjkoFnaQ3pfQuMLqO9tXA/8i/XwkMbuXSJEmSilKhe94kSZJ0BAxvkiRJGWJ4kyRJyhDDmyRJUoYY3iRJkjLE8CZJkpQhhjdJkqQMMbxJkiRliOFNkiQpQwxvkiRJGWJ4kyRJyhDDmyRJUoYY3iRJkjKkfaELkCRJrWTK4kJXoGZgz5skSVKGGN4kSZIyxPAmSZKUIYY3SZKkDHHAgiRJbci8L11Q6BLUwux5kyRJyhDDmyRJUoYY3iRJkjLE8CZJkpQhDliQJEnMGTOn0CXoMNnzJkmSlCH2vLUU14+TJEktwJ43SZKkDDG8SZIkZUhBw1tEfD4iNkTE3ogob+C4MRGxMSJeiYgbW7NGSZKkYlLoZ96eB64AflTfARFRAtwNfAqoBFZFxAMppRdap8T6uQSJJElqbQUNbymlFwEioqHDzgNeSSm9mj92LvBZoODhTZIkqbUVuuftcPQG3thvuxI4v0C1NAvn0pEkSY3V4uEtIn4DnFzHrm+mlH7dAt83FZgK0Ldv3+Y+vSRJUkG1eHhLKV3cxFNUAX322y7Nt9X3fbOB2QDl5eWpid8tSZJUVLIwVcgqYEBE9IuIY4ArgQcKXJMkSVJBFHqqkMsjohK4AFgcEY/k20+JiCUAKaU9wPXAI8CLwP0ppQ2FqlmSJKmQCj3a9FfAr+pofxP49H7bS4AlrViaJElSUcrCbVNJkiTlGd4kSZIyxPAmSZKUIYY3SZKkDDG8SZIkZYjhTZIkKUMMb5IkSRlieJMkScoQw5skSVKGGN4kSZIyxPAmSZKUIQVd21RHoSmLC12BJEmZZs+bJElShhjeJEmSMsTwJkmSlCGGN0mSpAwxvEmSJGWIo03V7OZ96YJClyBJUptleFPRmDNmTqFLkCSp6BneJBU35waUpFp85k2SJClDDG+SJEkZYniTJEnKEMObJElShjhgQVLBOb2MJB0+w5ukzHJ6GUlHI2+bSpIkZUhBw1tEfD4iNkTE3ogob+C41yPiuYhYFxGrW7NGSZKkYlLo26bPA1cAPzqMYz+ZUtrawvVIkiQVtYKGt5TSiwARUcgyJEmSMiMrz7wl4NGIqIiIqYUuRpIkqVBavOctIn4DnFzHrm+mlH59mKcZnlKqioiTgKUR8VJKaVk93zcVmArQt2/fRtUsSZLalsZOSVSMo9pbPLyllC5uhnNU5f98JyJ+BZwH1BneUkqzgdkA5eXlqanfLUmSVEwKPWDhkCKiC9AupfR+/v0lwLcLXJYkSWorpiwudAVHpNBThVweEZXABcDiiHgk335KRCzJH9YTeDIingV+DyxOKT1cmIolSZIKq9CjTX8F/KqO9jeBT+ffvwqc3cqlSZIkFaWsjDaVJEkSGXjmTZKkI5Kx55ekI2V4k6SsMqRIRyVvm0qSJGWIPW+SpMxp7ISrUltgz5skSVKGGN4kSZIyxNumkqSjRjGuUykdKXveJEmSMsSeN0kqYj6YL+lA9rxJkiRliOFNkiQpQwxvkiRJGeIzb5LUBjmqUmq77HmTJEnKEMObJElShhjeJEmSMsTwJkmSlCGGN0mSpAxxtKnUiho7W74jByVJ+9jzJkmSlCH2vEnFYsriQlcgScoAe94kSZIyxPAmSZKUIYY3SZKkDDG8SZIkZYjhTZIkKUMKGt4i4o6IeCki1kfEryKiWz3HjYmIjRHxSkTc2Np1SpIkFYtC97wtBc5KKQ0BNgE3HXhARJQAdwNjgTOBqyLizFatUpIkqUgUNLyllB5NKe3Jbz4NlNZx2HnAKymlV1NKu4C5wGdbq0ZJkqRiUuiet/39v8BDdbT3Bt7Yb7sy3yZJknTUafEVFiLiN8DJdez6Zkrp1/ljvgnsAX7ZDN83FZgK0Ldv36aeTpIkqai0eHhLKV3c0P6ImAxcBoxOKaU6DqkC+uy3XZpvq+/7ZgOzAcrLy+s6nyRJUmYVerTpGOB/AeNTSn+t57BVwICI6BcRxwBXAg+0Vo2SJEnFpNDPvN0FHAcsjYh1EfFDgIg4JSKWAOQHNFwPPAK8CNyfUtpQqIIlSZIKqcVvmzYkpXR6Pe1vAp/eb3sJsKS16pIkSSpWUfdjZm1DRGwB/ruFv6Y7sLWFvyOLvC7189rUzetSP69N3bwu9fPa1K3Yr8upKaUehzqoTYe31hARq1NK5YWuo9h4Xerntamb16V+Xpu6eV3q57WpW1u5LoV+5k2SJElHwPAmSZKUIYa3pptd6AKKlNelfl6bunld6ue1qZvXpX5em7q1ieviM2+SJEkZYs+bJElShhjeGikixkTExoh4JSJuLHQ9xSIifhoR70TE84WupZhERJ+IeDwiXoiIDRExvdA1FYuI6BQRv4+IZ/PX5v8UuqZiEhElEbE2Ih4sdC3FJCJej4jn8hO8ry50PcUiIrpFxPyIeCkiXoyICwpdUzGIiIH5f1f2vf4cEV8tdF2N5W3TRoiIEmAT8CmgktwSXlellF4oaGFFICJGAh8AP0spnVXoeopFRPQCeqWU1kTEcUAFMMF/ZyAiAuiSUvogIjoATwLTU0pPF7i0ohAR/xMoB/4mpXRZoespFhHxOlCeUirmObtaXUTcCyxPKf04v6Tkx1JK2wpdVzHJ/w6vAs5PKbX0XLAtwp63xjkPeCWl9GpKaRcwF/hsgWsqCimlZcB7ha6j2KSUNqeU1uTfv09uqbfeha2qOKScD/KbHfIv/1YJREQpMA74caFrUfGLiK7ASOAnACmlXQa3Oo0G/pDV4AaGt8bqDbyx33Yl/iLWYYqI04BhwDOFraR45G8NrgPeAZamlLw2Of8O/C9gb6ELKUIJeDQiKiJiaqGLKRL9gC3AnPyt9h9HRJdCF1WErgTuK3QRTWF4k1pRRBwL/Bfw1ZTSnwtdT7FIKVWnlIYCpcB5EXHU33KPiMuAd1JKFYWupUgNTymVAWOBL+cf2TjatQfKgB+klIYBfwF8Jns/+VvJ44H/v9C1NIXhrXGqgD77bZfm26R65Z/n+i/glymlBYWupxjlb/E8DowpdC1F4CJgfP7ZrrnA30fELwpbUvFIKVXl/3wH+BW5x1mOdpVA5X491/PJhTl9ZCywJqX0dqELaQrDW+OsAgZERL98ir8SeKDANamI5R/K/wnwYkrpu4Wup5hERI+I6JZ/35ncQKCXCltV4aWUbkoplaaUTiP3/5jfppSuLnBZRSEiuuQH/pC/LXgJcNSPcE8pvQW8ERED802jgaN+UNQBriLjt0wh18WqI5RS2hMR1wOPACXAT1NKGwpcVlGIiPuAUUD3iKgEbkkp/aSwVRWFi4B/BJ7LP9sF8L9TSksKWFOx6AXcmx8B1g64P6XktBhqSE/gV7m/E9Ee+M+U0sOFLalofAX4Zb5j4VVgSoHrKRr5oP8p4EuFrqWpnCpEkiQpQ7xtKkmSlCGGN0mSpAwxvEmSJGWI4U2SJClDDG+SJEkZYniTdFSKiG9GxIaIWB8R6yLi/HqOmxwRdzXTd74eEd2b41ySjl7O8ybpqBMRFwCXAWUppQ/zgeqYApclSYfFnjdJR6NewNaU0ocAKaWtKaU3I+LciFgZEc9GxO/3zeIPnBIRD0fEyxHxr/tOEhFXRcRzEfF8RPzLodolqTk4Sa+ko05EHAs8CXwM+A0wD3iK3LJck1JKqyLib4C/AlcDNwPDgA+BjcBwoBp4GjgH+BPwKPA94Pd1taeUFubXKS1PKW1tnZ9UUltkz5uko05K6QNy4WoqsIVcePsSsDmltCp/zJ9TSnvyH3kspbQ9pbST3FqRpwLnAk+klLbkj/slMLKBdklqFj7zJumolFKqBp4AnoiI54AvN3D4h/u9r8b/d0oqIHveJB11ImJgRAzYr2ko8CLQKyLOzR9zXEQ0FNJ+D/xdRHSPiBLgKuB3DbRLUrPwb4+SjkbHAv8REd2APcAr5G6hzsm3dwZ2ABfXd4KU0uaIuBF4HAhgcUrp1wD1tUtSc3DAgiRJUoZ421SSJClDDG+SJEkZYniTJEnKEMObJElShhjeJEmSMsTwJkmSlCGGN0mSpAwxvEmSJGXI/wVfOjoJS2RoaAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%pylab inline\n", "pylab.rcParams['figure.figsize'] = (10, 6)\n", "\n", "fig, ax = plt.subplots(nrows=1, sharex=True)\n", "\n", "ax.errorbar(np.array(range(8)), hmc_theta_prime_med,\n", " yerr=[hmc_theta_prime_med - hmc_theta_prime_low,\n", " hmc_theta_prime_hi - hmc_theta_prime_med], fmt='o')\n", "ax.errorbar(np.array(range(8)) + 0.1, klqp_theta_prime_med,\n", " yerr=[klqp_theta_prime_med - klqp_theta_prime_low,\n", " klqp_theta_prime_hi - klqp_theta_prime_med], fmt='o')\n", "ax.errorbar(np.array(range(8)) + 0.2, stan_theta_prime_med,\n", " yerr=[stan_theta_prime_med - stan_theta_prime_low,\n", " stan_theta_prime_hi - stan_theta_prime_med], fmt='o')\n", "ax.legend(('Edward/HMC', 'Edward/KLQP', 'Stan/NUTS'))\n", "plt.ylabel('$\\\\theta\\'$')\n", "plt.xlabel('School')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Comparison of the approximate posteriors for $\\theta$ and $y$\n", "\n", "We can graphically show the the three methods have very similar $95\\%$ intervals for $\\theta$:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAF3CAYAAAD6sAyZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4V+Wd///nTUCDqKhssilYAUVBiIF+VbAoKtiAom1B+3MsdEYdV/BrneqolPbSqTM6A3RGW/nVIl1GsJSiEWulLhUFFwIxgopYl8omoBI3YiDc3z8+IRBIkCXJOeTzfFwXFzn3WT7vHEFeue9z3yfEGJEkSVI6NUm6AEmSJNXOsCZJkpRihjVJkqQUM6xJkiSlmGFNkiQpxQxrkiRJKWZYkyRJSjHDmiRJUooZ1iRJklLMsCZJkpRiTZMuoC61bt06dunSJekyJEmSvlJRUdH6GGObrzquUYW1Ll26sHDhwqTLkCRJ+kohhPd257jUDIOGEHJCCItDCI9WbncNIbwYQngrhDAjhHBA0jVKkiQ1tNSENWAs8Pp22/8OTIwxHgt8DPxjIlVJkiQlKBVhLYTQCSgAflm5HYAzgZmVh0wDRiRTnSRJUnLS8szaJOBfgEMqt1sBG2KMmyu3VwAdkyhMkqT9yaZNm1ixYgVlZWVJl6JKubm5dOrUiWbNmu3V+YmHtRDCMGBtjLEohDBoL86/HLgc4Kijjqrj6iRJ2r+sWLGCQw45hC5dupAZqFKSYox8+OGHrFixgq5du+7VNdIwDHoacF4I4V1gOpnhz8nAYSGErWGyE7CyppNjjFNijPkxxvw2bb5y9qskSY1aWVkZrVq1MqilRAiBVq1a7VNPZ+JhLcZ4c4yxU4yxC3AR8FSM8f8Dnga+XXnY94CHEypRkqT9ikEtXfb1v0fiw6C78ENgegjhdmAxcH/C9UiSpD0wcOBAPv30U5o2beo6qPsgVWEtxvgM8Ezl128D/ZOsR5Ik7b158+YlXUKjkPgwqCRJkmpnWJMkSXVu0KBBvPHGGwB8+OGHnHjiiQlXtP8yrEmSpDr31ltv0b17dwBKSkro1atXwhXtv1L1zJq2GfP4GACmDp2acCWSpP3VuHHjKC4urtNr9unTh0mTJu3ymPfee4+OHTvSpEmmT6ikpITevXvXaR3ZxJ41SZJUp1555ZVq4ayoqMiwtg/sWWsAo+5bAMCMK05JuBJJUjb5qh6w+lJcXFy1COzy5ct5+OGHuf322xOppTGwZ02SJNWpV155hS1btnDSSSfxk5/8hJ49ezJt2rSky9pv2bMmSZLqVElJCYsWLeKQQw5JupRGwZ41SZJUZz799FNCCAa1OmRYkyRJdeaQQw7hzTffTLqMRsWwJkmSlGKGNUmSpBQzrEmSJKWYYU2SJCnFDGuSJEkpZliTJElKMcOaJElSihnWJElSncrJyaFPnz5Vv+68886djnnmmWcYNmxYvdYxYcIE7r777qrtF154gcsuu6zGzx49ejQzZ84EYNCgQRx11FHEGKv2jxgxgoMPPrhq+8033+Sb3/wm3bp1Iy8vj5EjR/LBBx/Uy/fh66YkSVKdat68OcXFxQ36mZs3b6Zp013Hmj/96U8MHTp0t6532GGH8fzzzzNgwAA2bNjA6tWrq/aVlZVRUFDAf/3XfzF8+HAgEz7XrVtHu3bt9v6bqIU9a5IkZbHZi1dy2p1P0fWmOZx251PMXryy3j7r8ccf57jjjiMvL49Zs2ZVtffq1YsNGzYQY6RVq1b8+te/BuDSSy9l7ty5vPvuuwwcOJC8vDzy8vKYP38+kAlIAwcO5LzzzqNnz54A3HHHHXTv3p0BAwawbNmyap//5JNPctZZZ+1WrRdddBHTp08HYNasWVx44YVV+/73f/+XU045pSqoQaY37sQTT9yLu/LVDGuSJGWp2YtXcvOsV1m5YSMRWLlhIzfPenWfA9vGjRurDYPOmDGDsrIyLrvsMgoLCykqKmLNmjVVx5922mk8//zzLF26lGOOOYZ58+YBsGDBAk499VTatm3L3LlzWbRoETNmzOC6666rOnfRokVMnjyZN998k6KiIqZPn05xcTGPPfYYL7/8ctVx69evp1mzZrRs2RKAefPmVavxkUceqfY9DB48mGeffZaKigqmT5/OqFGjqvYtWbKEk08+eZ/u0Z5wGDRpUwsyv4+Zk2wdkqSsc9efl7FxU0W1to2bKrjrz8sY0bfjXl+3pmHQ4uJiunbtSrdu3QC45JJLmDJlCgADBw7k2Wef5eijj+bKK69kypQprFy5ksMPP5wWLVpQWlrKNddcQ3FxMTk5OdXePdq/f3+6du0KZALYBRdcwEEHHQTAeeedV3XcE088wTnnnFO1PXDgQB599NGq7dGjR1erNycnhwEDBjB9+nQ2btxIly5d9vp+7Ct71iRJylKrNmzco/b6cvrppzNv3jzmzZvHoEGDaNOmDTNnzmTgwIEATJw4kXbt2vHKK6+wcOFCysvLq85t0aLFbn3GnjyvttVFF13Eddddx8iRI6u1n3DCCRQVFe3RtfaFYU2SpCzV4bDme9S+L4477jjeffdd/va3vwHw4IMPVu3r3Lkz69evZ/ny5RxzzDEMGDCAu+++m9NPPx2A0tJS2rdvT5MmTfjNb35DRUVFjZ9x+umnM3v2bDZu3Minn35KYWEhADFGSkpK6NOnzx7VPHDgQG6++WYuvvjiau3f/e53mT9/PnPmbBsVe/bZZ1myZMkeXX93GdYkScpSNw7pQfNmOdXamjfL4cYhPfbpujs+s3bTTTeRm5vLlClTKCgoIC8vj7Zt21Y75+tf/zrdu3cHMiFp5cqVDBgwAICrrrqKadOmcdJJJ/HGG2/U2puWl5fHqFGjOOmkkzj33HPp168fAEVFRfTt25cQwh59HyEEfvCDH9C6detq7c2bN+fRRx/lv//7v+nWrRs9e/bk3nvvpU2bNnt0/d2uY/s1RPZ3+fn5ceHChUmXsZNR9y0AYMYVp+y8s5Zn1sY8Piaze+jUeq1NktS4vP766xx//PG7ffzsxSu568/LWLVhIx0Oa86NQ3rs0/NqaXT77bdz7LHHctFFFyVWQ03/XUIIRTHG/K861wkGkiRlsRF9Oza6cLajW2+9NekS9onDoJIkSSlmWJMkSUoxw5okSVKKGdYkSZJSzLAmSZKUYoY1SZJUp3Jycqqts3bnnXfudMwzzzzDsGHD6rWOCRMmcPfdd1dtv/DCC1x22WU7ffatt97K0KFD+fLLLxk0aBA1LQP23HPP0b9/f4477jh69OjBvffeW+1zOnbsSJ8+fTjxxBN3es/ovnLpDkmSVKdqejdofdu8eTNNm+461tT0yqnbb7+d559/nscee4wDDzywxvPWrFnDd7/7XWbPnk1eXh7r169nyJAhtG/fngsuuACA66+/nh/84Ae8/vrrDBw4kLVr19KkSd30idmzJklSNit5CCaeCBMOy/xe8lC9fdTjjz/OcccdR15eHrNmzapq79WrFxs2bCDGSKtWrfj1r38NwKWXXsrcuXN59913GThwIHl5eeTl5TF//nwg0zs3cOBAzjvvPHr27AnAHXfcQffu3RkwYADLli2r9vlPPvkkZ511VtX2f/7nf/KnP/2JwsJCmjev/RVb99xzD6NHjyYvLw+A1q1b8x//8R/cddddOx17/PHH07RpU9avX7+Xd2ln9qxJkpStSh6CwutgU+WL20vfz2wD9B5Z+3lfYevrpra6+eabOf/887nssst46qmnOPbYYxk1alTV/tNOO43nn3+eo48+mmOOOYZ58+Zx6aWXsmDBAn7+858TQmDu3Lnk5uayfPlyLr744qqhykWLFrFkyRK6du1KUVER06dPp7i4mM2bN5OXl8fJJ58MwPr162nWrBktW7YE4Pnnn2fZsmUUFRVx8MEH7/L7Wbp0Kd/73veqteXn5/Paa6/tdOyLL75IkyZN6vTVU4Y1SZKy1ZM/2RbUttq0MdO+D2GtpmHQ4uJiunbtSrdu3QC45JJLmDJlCpB5F+izzz7L0UcfzZVXXsmUKVNYuXIlhx9+OC1atKC0tJRrrrmG4uJicnJyePPNN6uu279/f7p27QrAvHnzuOCCCzjooIMAOO+886qOe+KJJzjnnHOqto899lg+/vhj5s6dy7e+9a29/l63mjhxIr/97W855JBDmDFjxh6/h3RXHAaVJClbla7Ys/Z6cvrppzNv3jzmzZvHoEGDaNOmDTNnzmTgwIFAJgi1a9eOV155hYULF1JeXl51bm0vdd/Rjs+rtWvXjscee4xx48bx9NNP7/Lcnj17UlRUVK2tqKiI/Pxtr/W8/vrrKS4uZt68eVV11xXDmiRJ2aplpz1r3wfHHXcc7777Ln/7298AePDBB6v2de7cmfXr17N8+XKOOeYYBgwYwN13383pp58OQGlpKe3bt6dJkyb85je/oaKiosbPOP3005k9ezYbN27k008/pbCwEIAYIyUlJdWGZgG6d+/OrFmzuOSSS3Y5IeLqq6/mgQceqDrmww8/5JZbbuG2227b+xuyBwxrkiRlq8HjodkOD9Y3a55p3wdbn1nb+uumm24iNzeXKVOmUFBQQF5eHm3btq12zte//nW6d+8OZIZFV65cyYABAwC46qqrmDZtGieddBJvvPFGrb1peXl5jBo1ipNOOolzzz2Xfv36AZlesL59+9Y4NNmvXz+mTp3KeeedVxUkCwoK6NSpE506deI73/kO7du357e//S2XX345PXr0oEOHDlx33XV84xvf2Kf7tLtCjLFBPqgh5Ofnx5rWRknaqPsWADDjilN23jm1IPP7mDnVmsc8Piaze+jUeq1NktS4vP766xx//PG7f0LJQ5ln1EpXZHrUBo/fp+fV0uj222/n2GOP5aKLLqqT69177738/Oc/59lnn+Xwww/frXNq+u8SQiiKMebXckoVJxhIkpTNeo9sdOFsR7feemudXu+qq67iqquuqtNr7kriw6AhhNwQwkshhFdCCEtDCD+ubO8aQngxhPBWCGFGCOGApGuVJElqaImHNeBL4MwY40lAH2BoCOH/AP8OTIwxHgt8DPxjgjVKkpQuUwu2PUqjRi3xsBYzPqvcbFb5KwJnAjMr26cBIxIoT5Ik7afeKX2Hd0rfSbqMfZZ4WAMIIeSEEIqBtcBc4G/Ahhjj5spDVgAdk6pPkiQpKakIazHGihhjH6AT0B84bnfPDSFcHkJYGEJYuG7dunqrUZKk/UVj6VFSRirC2lYxxg3A08ApwGEhhK2zVTsBK2s5Z0qMMT/GmF+X7+GSJClpo+5bULX80/7kjjvu4IQTTqB379706dOHF198kUmTJvHFF1/s03U3bdpU9TL1EAI33HBD1b67776bCRMmADB69GhmzpxZ7dyDDz6YV199tWrttyOOOIKuXbvSp08fzjrrLLZs2cJ1113HiSeeSK9evejXrx/vvJOOwJv40h0hhDbAphjjhhBCc+BsMpMLnga+DUwHvgc8nFyVkiRpdyxYsIBHH32URYsWceCBB7J+/XrKy8sZNWoUl1xySdV7O/fU39Z9xoLnnuW0004D4MADD2TWrFncfPPNtG7dereu0atXr6q3EIwePZphw4bx7W9/G8i8UWHVqlWUlJTQpEkTVqxYsduvsqpvaehZaw88HUIoAV4G5sYYHwV+CPzfEMJbQCvg/gRrlCSpUZrz9hzOmXkOvaf15pyZ5zDn7TlffdIurF69mtatW3PggQcC0Lp1a2bOnMmqVas444wzOOOMMwC48soryc/P54QTTuBHP/pR1fldunThRz/6EXl5efTq1Ys33nijat+8p+Zy7rnnAtC0aVMuv/xyJk6cuE/1bl/31ldaAXTq1Gm3F7ytb4mHtRhjSYyxb4yxd4zxxBjjTyrb344x9o8xHhtj/E6M8cuka5UkqTGZ8/YcJsyfwOrPVxOJrP58NRPmT9inwHbOOefw/vvv0717d6666ir++te/ct1119GhQweefvrpqpem33HHHSxcuJCSkhL++te/UlJSUnWN1q1bs2jRIq688kruvvvuqvYXns+86H2rq6++mt/97neUlpbudb1bjRw5ksLCQvr06cMNN9zA4sWL9/madSXxsCZJkpIxedFkyirKqrWVVZQxedHkvb7mwQcfTFFREVOmTKFNmzaMGjWKBx54YKfjHnroIfLy8ujbty9Lly7ltddeq9p34YUXAnDyySfz7rvvArBm9SpaHnZ4tWHUQw89lEsvvZSf/exn1a5d0ztAa2rbXqdOnVi2bBk//elPadKkCYMHD+bJJ5/c3W+7XiX+zJokSUrGms/X7FH77srJyWHQoEEMGjSIXr16MW3atGr733nnHe6++25efvllDj/8cEaPHk1Z2bbQuHUINScnh82bM6t4PfvUXAaeMXinzxo3bhx5eXmMGTOmqq1Vq1Z8/PHHVdsfffTRbj3XduCBB3Luuedy7rnn0q5dO2bPns3gwTt/ZkOzZy1BpYWFfPFKMZ+//BLLzxxMaWFh0iVJkrLIkS2O3KP23bFs2TKWL19etV1cXMzRRx/NIYccwqeffgrAJ598QosWLWjZsiUffPABf/rTn77yus8+9Re+MfjsndqPOOIIRg4/m/v///uq2gYNGsSMGTMoLy8H4IEHHqh6Vq42ixYtYtWqVQBs2bKFkpISjj766K/+hhuAPWsJKS0sZPVt4+l8auYP0uZVq1h923gAWg4fnmRpkqQsMTZvLBPmT6g2FJqbk8vYvLF7fc3PPvuMa6+9lg0bNtC0aVOOPfZYpkyZwoMPPsjQoUOrnl3r27cvxx13HJ07d66a4VmbiooK3nvnbb7WrUeN+2+48vv8z/2/rdoeNmwYRUVFnDfoPHKa5HB89+P5xS9+scvPWLt2LZdddhlffpl5RL5///5cc801e/jd1w/DWkLWTpxELKv+nEAsK2PtxEmGNUlSgyg4JvNu0cmLJrPm8zUc2eJIxuaNrWrfGyeffDLz58/fqf3aa6/l2muvrdqu6Tk2oOoZNYD8/HyeeeYZnnvuOfrk5Vc77rPPPqv6ul3b1nzx9xJo3a2q7Uc/+hGXjrsUgK4tu+70OTt+/tChQxk6dGit31eSDGsJ2bx69R61S5JUHwqOKdincNYQBgwYQPsefZIuIzE+s5aQpu3b71G7JEnKToa1hLS9fhwhN7daW8jNpe314xKqSJIkpZHDoAnZ+lxaePKfieXlNO3QgbbXj/N5NUmSVI1hLUEthw+H9ZnZKd3u3bfXe0iSpMbJYVBJkqQUM6xJkqQ6tWLFCs4//3y6devG1772NcaOHVttgdq0rF+2vYMPPjjpEmplWJMkSXUmxsiFF17IiBEjWL58OW+++SafffYZt9xyS7195tZXUjVWhjVJkrJYaWEhy88czOvH96yTVx8+9dRT5ObmVr2rMycnh4kTJ/KrX/2KL774AoD333+fQYMG0a1bN3784x8D8Pnnn1NQUMBJJ53EiSeeyIwZMwAoKiriG9/4BuefNZDRI0ewunI90kGDBjFu3Djy8/O5Y+LPObrvILZs2VJ1rc6dO7Np0ybee+c9hg4dysknn8zAgQN54403gMz7SU855RR69erFrbfeWuP3Mn78eCZNmlS1fcsttzB58t6/5H5vOcFAkqQstfXVh1vfqFMXrz5cunQpJ598crW2Qw89lKOOOoq33noLgJdeeoklS5Zw0EEH0a9fPwoKCnjvvffo0KEDc+ZkJtyVlpayadMmrr32Wh5++GE+oTlzZv+BW265hV/96lcAlJeXs3DhQli/nEUlr/HXv/6VM844g0cffZQhQ4bQrFkz/nXsvzLtl9Po1q0bL774IldddRVPPfUUY8eO5corr+TSSy/lnnvuqfF7+f73v8+FF17IuHHj2LJlC9OnT+ell17aq/uyLwxrkiRlqaRefXj22WfTqlUrAC688EKee+45vvnNb3LDDTfwwx/+kGHDhjFw4ECWLFnCkiVLOPvssynfvIWKLRUc3alj1XVGjRq17esR32TGjBmcccYZTJ8+nauuuorPP/ucRS8t4jvf+U7VcVvf/fn888/zhz/8AYB/+Id/4Ic//OFOdXbp0oVWrVqxePFiPvjgA/r27VtVd0MyrEmSlKXq49WHPXv2ZObMmdXaPvnkE/7+979z7LHHsmjRIkII1faHEOjevTuLFi3iscce49Zbb2Xw4MFccMEFnHDCCSxYsIC/rcu8C/RrbbZNBGjRokXV1+cNOZN//enP+OijjygqKuLMM89k6aqlHNryUIqLi2usdcc6avJP//RPPPDAA6xZs4bvf//7u30f6pLPrEmSlKXq49WHgwcP5osvvuDXv/41ABUVFdxwww2MHj2agw46CIC5c+fy0UcfsXHjRmbPns1pp53GqlWrOOigg7jkkku48cYbWbRoET169GDdunUsWLAAgE2bNrF06dIaP/fgg1vQr18/xo4dy7Bhw8jJyeGQQw+h89Gd+f3vfw9kJj+88sorAJx22mlMnz4dgN/97ne1fj8XXHABjz/+OC+//DJDhgzZ6/uyLwxrkiRlqfp49WEIgT/+8Y/8/ve/p1u3bnTv3p3c3Fz+7d/+reqY/v37861vfYvevXvzrW99i/z8fF599VX69+9Pnz59+PGPf8ytt97KAQccwMyZMzNDo4NOYfiZpzJ//vxaP3vUqFH89re/rTY8OnHKRO6//35OOukkTjjhBB5++GEAJk+ezD333EOvXr1YuXJlrdc84IADOOOMMxg5ciQ5OTl7fV/2hcOgkiRlqa3Ppa2dOInNq1fTtH37Onn1YefOnSmsZVbp6NGjGT169E7tQ4YMqbHnqk+fPjz77LM7DYM+88wzOx377W9/mxhj9Vq6dObxxx/f6diuXbtW9dgB3H777TXWu2XLFl544YWq3rkkGNYkScpiLYcP973UtXjttdcYNmwYF1xwAd26dUusDsOaJElSDXr27Mnbb7+ddBk+syZJkpRmhjVJkhqZHZ/bUrL29b+HYU2SpEYkNzeXstIyA1tKxBj58MMPyd1h1u2e8Jk1Sek1tSDz+5g5ydYh7Uc6derEo4sepfWHrfmy+ZdJl1Nn1n2a+V7K1x+4887P1lYeVP2F7us3rgegrHnZjmc0qNzcXDp16rTX5xvWlA7+oyxJdaJZs2Y88vEjAEwdOjXhaurOhPsyy2zMuKLPzjun/iDz+w7/hox5PPMy+f39PjgMKkmSlGKGNUmSpBQzrEmSJKWYYU2SJCnFDGtKvTGPj6l6SDQrTC3YNuEiC4y6bwGj7lvw1QfuIOv+XEjKWoY1Nai9/YdZkqRsZViT9iP2JklS9jGsSZIkpZhhTZIkKcUMa1ICfHZPkrS7DGuSJEkpZliTJElKMcOaJElSihnWJEmSUsywJkmSlGKJh7UQQucQwtMhhNdCCEtDCGMr248IIcwNISyv/P3wpGuVJCUsy17HJkEKwhqwGbghxtgT+D/A1SGEnsBNwJMxxm7Ak5XbkiRJWSXxsBZjXB1jXFT59afA60BH4HxgWuVh04ARyVQoSZKUnMTD2vZCCF2AvsCLQLsY4+rKXWuAdgmVJUmSlJjUhLUQwsHAH4BxMcZPtt8XY4xArOW8y0MIC0MIC9etW9cAlUqSJDWcVIS1EEIzMkHtdzHGWZXNH4QQ2lfubw+srencGOOUGGN+jDG/TZs2DVOwJElSA0k8rIUQAnA/8HqM8b+22/UI8L3Kr78HPNzQtUmS9h9jHh/DmMfHJF2GVOeaJl0AcBrwD8CrIYTiyrZ/Be4EHgoh/CPwHjAyofokSQ1s1H0LAJhxxSkJVyIlL/GwFmN8Dgi17B7ckLVIkiSlTeLDoJIkSaqdYU2SJCnFDGuSJO1nSgsL+eKVYj5/+SWWnzmY0sLCpEtSPTKsSZK0HyktLGT1beOJ5eUAbF61itW3jTewNWKGNUmS9iNrJ04ilpVVa4tlZaydOCmhilTfDGuSJO1HNq9evUft2v8Z1iRJ2o80bd9+j9q1/zOsSdL+YGpB5peyXtvrxxFyc6u1hdxc2l4/LqGKVN8SXxRXkiTtvpbDhwMQnvxnYnk5TTt0oO3146ra1fgY1iRJ2s+0HD4c1v8CgG73zkm4GtU3h0ElSZJSzLAmSZKUYoY1SZKkFDOsSZIkpZhhTZIkKcUMa5IkSSlmWJOk/diYx8cw5vExSZchqR4Z1iQpJUbdt4BR9y1IugxJKWNYkyRJSjHDmiRJUooZ1iRJklLMsCZJkpRihjVJkqQUM6yl0Jy351CyroSFHyzknJnnMOftOUmXJEmSEmJYS5k5b89hwvwJlG8pB2D156uZMH+CgU2SpCxlWKtnsxevZPHfN/DiOx9x2p1PMXvxyl0eP3nRZMoqyqq1lVWUMXnR5PosU5IkpZRhrR7NXrySm2e9SnnFFgBWbtjIzbNe3WVgW/P5mj1qlyRJjZthrR7d9edlbNxUUa1t46YK7vrzslrPObLFkXvULkmSGjfDWj1atWHjHrUDjM0bS25ObrW23JxcxuaNrdPaJEnS/sGwVo86HNZ8j9oBCo4pYMKpEzigyQEAtG/RngmnTqDgmIJ6qVGSJKWbYa0e3TikB82b5VRra94shxuH9MhslDwEK16G956DiSdmtskEtt5tepPfLp8nvv2EQU2SpCxmWKtHI/p25KcX9uKAnMxt7nhYc356YS9G9O2YCWaF10HFl5mDS9/PbFcGNklSdtvT1QS2cq3Oxqdp0gU0diP6duTBl/4OwIwrTtm248mfwKYdnl3btDHT3ntkA1YoSUqb2lYTgMy/K7Wpba1OwFGa/Zg9a0kpXbFn7ZKkrLE3qwlAdq7VWVpYyBevFPP5yy+x/MzBlBYWJl1SnTOsJaVlpz1rlyRljb1ZTQCyb63O0sJCVt82nlie6UncvGoVq28b3+gCm2EtKYPHQ7MdZoU2a55plyRlta9cTaCWCWrZtlbn2omTiGXVexJjWRlrJ05KqKL6YVhLSu+RMPxnkHNgZrtl58y2z6tJUtbb5WoCu5iglm1rdW5evXqP2vdXTjBIUu+RUDQt8/UYZ+tIkjK2TiL4l5kllFdsoeNhzblxSI9M+8TaJ6hp9IS2AAAV40lEQVQVXL8EgPHPj6d8SzntW7RnbN7YRju5oGn79mxetarG9sbEsCZJUgrVuprAV0xQKzimgJlvzgRg6tCp9Vpj0tpeP47Vt1V/fCjk5tL2+nEJVVQ/DGuSJO1PWnbKDH3W1J5lWg4fDkB48p+J5eU07dCBttePq2pvLAxrkiTtTwaPzzyjtv1QaBZPUGs5fDis/wUA3e5tnI8UOcFAkqT9iRPUsk4qetZCCL8ChgFrY4wnVrYdAcwAugDvAiNjjB8nVaMkSanhBLWskpaetQeAoTu03QQ8GWPsBjxZuS1JkpRVUhHWYozPAh/t0Hw+UPljA9OAEQ1alCRJUgqkIqzVol2MceuqdmuAdkkWI0mSlIQ0h7UqMcYIxJr2hRAuDyEsDCEsXLduXQNXJkmSVL/SHNY+CCG0B6j8fW1NB8UYp8QY82OM+W3atGnQAiVJkurbV84GDSF0Aa4GvkbmubJioDDG+F69VgaPAN8D7qz8/eF6/jxJkqTU2Z2etYeBN4B7gLOBk4BnQwj3hBAOrIsiQggPAguAHiGEFSGEfyQT0s4OISwHzqrcliRJyiq7s85aTozxfoAQwkcxxstCCE2B64EpZHq99kmM8eJadg3e12tLkiTtz3anZ+0vIYRrKr+OADHGzTHGu4BTaj9NkiRJ+2p3etb+L3BzCGEh0CGEcDnwBZmg9mF9FidJkpTtvrJnLca4JcZ4B3A6cDlwJHAysAQ4t37LUzYoLSzki1eK+fzll1h+5mBKCwuTLkmSpNTY7XeDxhi/IDND85H6K0fZprSwkNW3jafzqeUAbF61itW3jQeg5fDhSZYmSVIqpHmdNWWBtRMnEcvKqrXFsjLWTpyUUEWSpMZgzttzKFlXwsIPFnLOzHOY8/b++8J7w5oStXn16j1qlyRln9mLV7L47xt48Z2POO3Op5i9eOW2nSUPwYqX4b3nYOKJUPIQc96ew4T5Eyjfkhm1Wf35aibMn7DfBrbdHgaV6kPT9u3ZvGpVje2SJM1evJKbZ71KecUWAFZu2MjNs14FYETO81B4HVR8mTm49H0ovI7JXY+lrKL6qE1ZRRmTF02m4JiCBq2/LtizpkS1vX4c4YBm1drCAc1oe/04oHF1Y2vPOPFEEsBdf17Gxk0V1do2bqrgrj8vgyd/Aps2Vj9h00bWlG+o8VprPl9TX2XWK8OaGkxN3dgtj95I+34fE5pEAJoetJn2/T6m5dEbG103tnbf1oknsbz6xBMDm5R9Vm3YWHt76Yoa9x25uaLm9hZH1lldDcmwpgZRWzf2F38aT8vOn3BQ6020aFtOt/PW0rLzJ/DkT5i8aHKt3dhq3L5q4ok9rlL26HBY89rbW3aqcd/YL3PIzcmt1pabk8vYvLF1Xl9DMKypQdTWjZ27sZYu6dIVtXZX76/d2Np9u5p4Yo+rlF1uHNKD5s1yqrU1b5bDjUN6wODx0GyHMNesOQUDxzPh1Akc0OQAANq3aM+EUyfsl8+rgWFNDaTWbuwtrWo+oWWnWrur99dubO2sthleTVsdWuPxTVsdao+rlGVG9O3ITy/sxQE5mcjS8bDm/PTCXozo2xF6j4ThP4OcAzMHt+yc2e49koJjCujdpjf57fJ54ttP7LdBDQxraiC1dWP/8oBLavypiMHjGZs3tlF1Y6u62obGZy9eSdvenxBytlQ7PuRsoW3vT+xxlbLQiL4d6XvUYXy96xE8f9OZmaC2Ve+R0KkfHD0Arl+S2W5kDGtqELV1Y/cpuHyXPxU1pm5sVberGV4t266ifb/SHSaelNKy7Sp7XCVlHddZU4PY+lPQv8wsobxiCx0Pa86NQ3pUto+EommZA8dUf+6o4JgCZr45E4CpQ6c2ZMmqZ7uc4dWuEy27vF/V1u28tZkvWnZmbN5YJsyfUG0o1B5XSY2ZPWtqMLvsxlbW2eUMr1oeGmbweHtcJWUdw5qkROxyhtcuHhoGGtWDw5L0VRwGlZSIXQ+NkwlmtQyPS1I2CTHGpGuoM/n5+XHhwoX1+hnjxo2juLh4j855bdUnAPTsUMNyBGsy7zfjyF7Vmt/46A0AjjviuD0vMsVqvRe13AdonPdib/5MgPdie96LjMZ4H8D/V2zl349tGvrvR58+fZg0adIen7cnQghFMcb8rzrOYVBJkqQUcxh0D+1Nyh513wIAZlxxys47p1Y+a7PDMM+Yx8dkdjeyGZC13ota7gM0znuxN38mwHuxPe9F5WYjvA/g/yu28u/HNtn898OeNUmSpBQzrEmS9g8lD8GKl+G952DiiZltKQsY1iRJ6VfyEBReBxVfZrZL389sG9iUBQxrUprYcyDV7MmfwKYd3nqxaWOmXWrkDGtSWthzIAEwe/FKFv99Ay++8xGn3fkUsxevhNIVNR9c2T7n7TmUrCth4QcLOWfmOcx527X51HgY1qS0sOdAYvbildw861XKK7YAsHLDRm6e9SpfND+y5hNadmLO23OYMH8C5VvKAVj9+WomzJ9gYFOjYViT0sKeA4m7/ryMjZsqqrVt3FTBf2waVev7YicvmkxZRVm1XWUVZUxeNLm+y5UahGFNSouWnWptt+dA2WLVho01tk/7rH+t74td8/maGs+prV3a3xjWpAZW4/M4AIPH23OgrNfhsOa1t/ceCZ36wdED4PolmW3gyBY1D5HW1i7tbwxrUgOq7Xmc2YtXZv7hsedANcmiWcI3DulB82Y51dqaN8vhxiE9aj1nbN5YcnNyq7Xl5uQyNm9svdQoNTTDmtSAanse564/L8ts2HOgHWXZLOERfTvy0wt7cUBO5p+njoc156cX9mJE3461nlNwTAETTp3AAU0OAKB9i/ZMOHUCBccUNEjNUn0zrEkNqLbncWpr38qegyy2i1nCjXXSyYi+Hel71GF8vesRPH/TmbsMalsVHFNA7za9yW+XzxPffsKgpkbFsCY1oF0+j7ML9hxksVpmCc/Z/JGTTqQsYViTGtDePI+zlT0HjVutE09qmSU8udURTjqRsoRhTWpAe/M8jhq/XU48qWWW8JqcUOO1nHQiNT6GNamB7c3zOGrcdjnxpJZZwke2aF/jtZx0IjU+hjVJSthXTjypYZawk06k7GFYk6SE7c3Ek2ycdFJaWMgXrxTz+csvsfzMwZQWFiZdktQgDGuSlLC9nXiSTZNOSgsLWX3beGJ5Zvbr5lWrWH3beAObsoJhTZIS5sSTr7Z24iRiWfXZr7GsjLUTJyVUkdRwmiZdgCQpE9gefOnvAMy44pSEq0mfzatX71G71JikvmcthDA0hLAshPBWCOGmpOuRJDW8pu1rnv1aW7vUmKQ6rIUQcoB7gHOBnsDFIYSeyVYlSWpoba8fR8itPvs15ObS9vpxCVUkNZy0D4P2B96KMb4NEEKYDpwPvJZoVZKkBtVy+HAAwpP/TCwvp2mHDrS9flxVu9SYpT2sdQTe3257BfD1hGqRJCWo5fDhsP4XAHS713egKnukehh0d4QQLg8hLAwhLFy3bl3S5UiSJNWptIe1lUDn7bY7VbZViTFOiTHmxxjz27Rp06DFSZIk1be0h7WXgW4hhK4hhAOAi4BHEq5JkhqUK/dL2S3VYS3GuBm4Bvgz8DrwUIxxabJVSVLDceV+SakOawAxxsdijN1jjF+LMd6RdD2S1JBcuV9S6sOaJGUzV+6XZFiTpBRz5X5JhjVJSjFX7peU9kVxJSmruXK/JMOaJKWcK/dL2c1hUEmSpBQzrEmSJKWYw6BSipQWFtLslWJieTmrfjPYZ5OkLDfjilOSLkEpYFiTUmLrSvWdT62+Uj1gYJO0R6YOnZp0CapDDoNKKeFK9ZKkmhjWpJRwpXpJUk0Ma1JKuFK9JKkmhjUpJVypXpJUEycYSCnhSvWS9sgYF0jOFoY1KUVcqV6StCOHQSVJklLMsCZJkpRihjVJkqQU85m1BrA3rwtx9WlJkgT2rEmSJKWaYU2SJCnFHAZNmuvkSJKkXbBnTZIkKcUMa5IkSSlmWJMkSUoxw5okSVKKGdYkSZJSzLAmSZKUYoY1SZKkFDOsSZIkpZhhTZIkKcUMa5IkSSlmWJMkSUoxw5okSVKKGdYkSZJSzLAmSZKUYoY1SZKkFDOsSZIkpVjTpAtQdplxxSl7fM7UoVProRJJkvYP9qxJkiSlmGFNkiQpxQxrkiRJKZZoWAshfCeEsDSEsCWEkL/DvptDCG+FEJaFEIYkVaMkSVKSkp5gsAS4ELhv+8YQQk/gIuAEoAPwlxBC9xhjRcOXqAYxZk7SFUiSlEqJ9qzFGF+PMS6rYdf5wPQY45cxxneAt4D+DVudJElS8pLuWatNR+CF7bZXVLZJEuCSLpKyR72HtRDCX4Aja9h1S4zx4Tq4/uXA5QBHHXXUvl5OkiQpVeo9rMUYz9qL01YCnbfb7lTZVtP1pwBTAPLz8+NefJYkSVJqpXXpjkeAi0IIB4YQugLdgJcSrkmSJKnBJfrMWgjhAuC/gTbAnBBCcYxxSIxxaQjhIeA1YDNwtTNBpcZpl68gc5awJCUb1mKMfwT+WMu+O4A7GrYiSZKkdEnrMKgkSZIwrEmSJKWaYU2SJCnFDGuSJEkpZliTJElKMcOaJElSihnWJEmSUsywJkmSlGKJLoorSZK0zxr5207sWZMkSUoxw5okSVKKGdYkSZJSzGfWJGk/NnXo1KRLkFTPDGuSJKlRaiw/zBjWJCklZlxxStIl7Ncayz/M0o4Ma5IkKfWy+YcZw5okKXVq/Ye5ka+nJdXE2aCSJEkpZliTJElKMYdBJWl/4PCflLXsWZMkSUoxw5okSVKKGdYkSZJSzLAmSZKUYoY1SZKkFDOsSZIkpZhhTZIkKcUMa5IkSSlmWJMkSUoxw5okSVKKGdYkSZJSzLAmSZKUYr7IXUobX9gtSdqOPWuSJEkpZliTJElKMcOaJElSihnWJEmSUsywJkmSlGKGNUmSpBQzrEmSJKWYYU2SJCnFDGuSJEkplmhYCyHcFUJ4I4RQEkL4YwjhsO323RxCeCuEsCyEMCTJOiVJkpKS9Oum5gI3xxg3hxD+HbgZ+GEIoSdwEXAC0AH4Swihe4yxIsFapcRNHTo16RIkSQ0s0Z61GOMTMcbNlZsvAJ0qvz4fmB5j/DLG+A7wFtA/iRolSZKSlHTP2va+D8yo/LojmfC21YrKNqlRmHHFKUmXIEnaT9R7WAsh/AU4soZdt8QYH6485hZgM/C7vbj+5cDlAEcdddQ+VCpJkpQ+9R7WYoxn7Wp/CGE0MAwYHGOMlc0rgc7bHdapsq2m608BpgDk5+fHmo6RJEnaXyU9G3Qo8C/AeTHGL7bb9QhwUQjhwBBCV6Ab8FISNUqSJCUp6WfW/gc4EJgbQgB4Icb4zzHGpSGEh4DXyAyPXu1MUEmSlI0SDWsxxmN3se8O4I4GLEeSJCl1fIOBJElSihnWJEmSUsywJkmSlGKGNUmSpBQzrEmSJKWYYU2SJCnFDGuSJEkpZliTJElKMcOaJElSioVt707f/4UQ1gHvJV1HLVoD65MuIiW8Fxneh228F9t4L7bxXmR4H7ZpbPfi6Bhjm686qFGFtTQLISyMMeYnXUcaeC8yvA/beC+28V5s473I8D5sk633wmFQSZKkFDOsSZIkpZhhreFMSbqAFPFeZHgftvFebOO92MZ7keF92CYr74XPrEmSJKWYPWuSJEkpZlirZyGEoSGEZSGEt0IINyVdT5JCCL8KIawNISxJupYkhRA6hxCeDiG8FkJYGkIYm3RNSQkh5IYQXgohvFJ5L36cdE1JCiHkhBAWhxAeTbqWJIUQ3g0hvBpCKA4hLEy6niSFEA4LIcwMIbwRQng9hHBK0jUlIYTQo/LPw9Zfn4QQxiVdV0NxGLQehRBygDeBs4EVwMvAxTHG1xItLCEhhNOBz4BfxxhPTLqepIQQ2gPtY4yLQgiHAEXAiGz8cxFCCECLGONnIYRmwHPA2BjjCwmXlogQwv8F8oFDY4zDkq4nKSGEd4H8GGNjWk9rr4QQpgHzYoy/DCEcABwUY9yQdF1Jqvy3dSXw9RhjWtdWrVP2rNWv/sBbMca3Y4zlwHTg/IRrSkyM8Vngo6TrSFqMcXWMcVHl158CrwMdk60qGTHjs8rNZpW/svInyBBCJ6AA+GXStSgdQggtgdOB+wFijOXZHtQqDQb+li1BDQxr9a0j8P522yvI0n+UVbMQQhegL/BispUkp3LorxhYC8yNMWbrvZgE/AuwJelCUiACT4QQikIIlyddTIK6AuuAqZXD478MIbRIuqgUuAh4MOkiGpJhTUpICOFg4A/AuBjjJ0nXk5QYY0WMsQ/QCegfQsi6IfIQwjBgbYyxKOlaUmJAjDEPOBe4uvIRimzUFMgDfh5j7At8DmT7s88HAOcBv0+6loZkWKtfK4HO2213qmxTlqt8PusPwO9ijLOSricNKod3ngaGJl1LAk4Dzqt8Vms6cGYI4bfJlpScGOPKyt/XAn8k80hJNloBrNiut3kmmfCWzc4FFsUYP0i6kIZkWKtfLwPdQghdK38auAh4JOGalLDKh+rvB16PMf5X0vUkKYTQJoRwWOXXzclMxnkj2aoaXozx5hhjpxhjFzL/n3gqxnhJwmUlIoTQonLiDZVDfucAWTmDPMa4Bng/hNCjsmkwkHUTkXZwMVk2BAqZLlbVkxjj5hDCNcCfgRzgVzHGpQmXlZgQwoPAIKB1CGEF8KMY4/3JVpWI04B/AF6tfFYL4F9jjI8lWFNS2gPTKmd3NQEeijFm9bIVoh3wx8zPNDQF/jfG+HiyJSXqWuB3lT/wvw2MSbiexFSG97OBK5KupaG5dIckSVKKOQwqSZKUYoY1SZKkFDOsSZIkpZhhTZIkKcUMa5IkSSlmWJOUFUIIt4QQloYQSkIIxSGEr9dy3OgQwv/U0We+G0JoXRfXkpS9XGdNUqMXQjgFGAbkxRi/rAxQByRcliTtFnvWJGWD9sD6GOOXADHG9THGVSGEfiGE+SGEV0IIL21dOR/oEEJ4PISwPITwH1svEkK4OITwaghhSQjh37+qXZLqgoviSmr0QggHA88BBwF/AWYAC8i82mpUjPHlEMKhwBfAJcB4oC/wJbAMGABUAC8AJwMfA08APwNeqqk9xji78l2f+THG9Q3znUpqjOxZk9ToxRg/IxOmLgfWkQlrVwCrY4wvVx7zSYxxc+UpT8YYS2OMZWTexXg00A94Jsa4rvK43wGn76JdkuqEz6xJygoxxgrgGeCZEMKrwNW7OPzL7b6uwP9XSkqQPWuSGr0QQo8QQrftmvoArwPtQwj9Ko85JISwq1D2EvCNEELryhfPXwz8dRftklQn/GlRUjY4GPjvEMJhwGbgLTJDolMr25sDG4GzartAjHF1COEm4GkgAHNijA8D1NYuSXXBCQaSJEkp5jCoJElSihnWJEmSUsywJkmSlGKGNUmSpBQzrEmSJKWYYU2SJCnFDGuSJEkpZliTJElKsf8Hn9Nno0Kt6eoAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%pylab inline\n", "pylab.rcParams['figure.figsize'] = (10, 6)\n", "\n", "fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True)\n", "ax.scatter(np.array(range(8)), hmc_theta_med)\n", "ax.scatter(np.array(range(8)) + 0.1, klqp_theta_med)\n", "ax.scatter(np.array(range(8)) + 0.2, stan_theta_med)\n", "ax.scatter(np.array(range(8)) + 0.1, data_y)\n", "\n", "mu = hq_mu.params.eval()[burn:].mean()\n", "\n", "plt.plot([-0.2, 7.4], [mu, mu], 'k')\n", "\n", "ax.errorbar(np.array(range(8)), hmc_theta_med,\n", " yerr=[hmc_theta_med - hmc_theta_low,\n", " hmc_theta_hi - hmc_theta_med], fmt='none')\n", "\n", "ax.errorbar(np.array(range(8)) + 0.1, klqp_theta_med,\n", " yerr=[klqp_theta_med - klqp_theta_low,\n", " klqp_theta_hi - klqp_theta_med],\n", " fmt='none')\n", "\n", "ax.errorbar(np.array(range(8)) + 0.2, stan_theta_med,\n", " yerr=[stan_theta_med - stan_theta_low,\n", " stan_theta_hi - stan_theta_med],\n", " fmt='none')\n", "\n", "\n", "ax.legend(('$\\\\bar{\\\\mu}$', 'Edward/HMC', 'Edward/KLQP',\n", " 'Stan/NUTS', 'Observed y'))\n", "\n", "\n", "plt.xlabel('School')\n", "plt.ylabel('$\\\\theta$')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model makes an estimate of $\\theta_j$ for each school which combines information from the score for that school $y_j$ and the global average schore $\\mu$. The prior on $\\tau$ influences how much pooling occurs. If $\\tau$ is small then the pooling is strong and all schools have a similar $\\theta_j$ and the uncertanity on $\\theta_j$ decreases due to the sharing of information, if $\\tau$ is large all schools will have a different $\\theta_j$ with more uncertainty.\n", "\n", "For the prior used here we observe a compromise, the median of $\\theta_j$ is between $y_j$ and $\\bar{\\mu}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Appendix R code for running through rstan\n", "\n", "library(rstan)\n", "\n", "rstan_options(auto_write = TRUE)\n", "\n", "options(mc.cores = parallel::detectCores())\n", "\n", "schools_dat <- list(J = 8,\n", " y = c(28, 8, -3, 7, -1, 1, 18, 12),\n", " sigma = c(15, 10, 16, 11, 9, 11, 10, 18))\n", "\n", "\n", "fit <- stan(file = 'eight_schools.stan', data = schools_dat,iter = 100000, chains = 4)\n", "\n", "print(fit)" ] }, { "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.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }