Statistical Rethinking - Chapter 4 (Section 4.3) using R.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: rstan\n",
      "Loading required package: ggplot2\n",
      "Loading required package: StanHeaders\n",
      "rstan (Version 2.12.1, packaged: 2016-09-11 13:07:50 UTC, GitRev: 85f7a56811da)\n",
      "For execution on a local, multicore CPU with excess RAM we recommend calling\n",
      "rstan_options(auto_write = TRUE)\n",
      "options(mc.cores = parallel::detectCores())\n",
      "Loading required package: parallel\n",
      "rethinking (Version 1.59)\n"
     ]
    }
   ],
   "source": [
    "# R code 4.24\n",
    "\n",
    "library(rethinking)\n",
    "data(Howell1)\n",
    "\n",
    "d <- Howell1\n",
    "d2 <- d[ d$age >= 18, ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flist <- alist (\n",
    "    height ~ dnorm(mu, sigma),\n",
    "    mu ~ dnorm(178, 20),\n",
    "    sigma ~ dunif(0, 50)\n",
    ")\n",
    "\n",
    "m4.1 <- map(flist, data=d2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "        Mean StdDev   5.5%  94.5%\n",
       "mu    154.61   0.41 153.95 155.27\n",
       "sigma   7.73   0.29   7.26   8.20"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "precis(m4.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These number provide Gaussian approximations (normal/quadratic) for each parameter's marginal distribution. This means the plausibility of each value of $\\mu$, after averaging over the plausibilities of each value of $\\sigma$ is given by a Gaussian distribution with mean 154.6 and standard deviation 0.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m4.2 <- map(\n",
    "    alist(\n",
    "        height ~ dnorm(mu, sigma),\n",
    "        mu ~ dnorm(178, 0.1),\n",
    "        sigma ~ dunif(0, 50)\n",
    "    ),\n",
    "    data = d2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "        Mean StdDev   5.5%  94.5%\n",
       "mu    177.86   0.10 177.70 178.02\n",
       "sigma  24.52   0.93  23.03  26.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "precis(m4.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the estimate for $\\mu$ has hardly moved off the prior. The prior was very concentrated around 178, so this is not surprising. But also notice the estimate for $\\sigma$ has changed quite a lot, even though we didn't change its prior at all. Once the golem is certain that the mean is 178 -- as the prior insists -- then the golem has to estimate $\\sigma$ conditional on that fact. This results in a different posterior for $\\sigma$, even though all we changed was the prior for $\\mu$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 4.3.6 Sampling from a ```map``` fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like a mean and a variance are sufficient to describe a one-dimensional Gaussian distribution, a list of means and a matrix of variances and covariances are sufficient to describe a multi-dimensional Gaussian distribution. To see the matrix of variances and covariances for model m4.1: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>mu</th><th scope=col>sigma</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>mu</th><td>0.1696800422</td><td>0.0002247144</td></tr>\n",
       "\t<tr><th scope=row>sigma</th><td>0.0002247144</td><td>0.0848312814</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & mu & sigma\\\\\n",
       "\\hline\n",
       "\tmu & 0.1696800422 & 0.0002247144\\\\\n",
       "\tsigma & 0.0002247144 & 0.0848312814\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 0.169680042174359\n",
       "2. 0.000224714361676325\n",
       "3. 0.000224714361676325\n",
       "4. 0.0848312813634882\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "      mu           sigma       \n",
       "mu    0.1696800422 0.0002247144\n",
       "sigma 0.0002247144 0.0848312814"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vcov(m4.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A variance-covariance matrix can be factored into two elements:\n",
    "- a vector of variances for the parameters\n",
    "- a correlation matrix that tells us how changes in any parameter lead to correlated changes in the others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<dl class=dl-horizontal>\n",
       "\t<dt>mu</dt>\n",
       "\t\t<dd>0.169680042174359</dd>\n",
       "\t<dt>sigma</dt>\n",
       "\t\t<dd>0.0848312813634882</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[mu] 0.169680042174359\n",
       "\\item[sigma] 0.0848312813634882\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "mu\n",
       ":   0.169680042174359sigma\n",
       ":   0.0848312813634882\n",
       "\n"
      ],
      "text/plain": [
       "        mu      sigma \n",
       "0.16968004 0.08483128 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>mu</th><th scope=col>sigma</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>mu</th><td>1.000000000</td><td>0.001872999</td></tr>\n",
       "\t<tr><th scope=row>sigma</th><td>0.001872999</td><td>1.000000000</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & mu & sigma\\\\\n",
       "\\hline\n",
       "\tmu & 1.000000000 & 0.001872999\\\\\n",
       "\tsigma & 0.001872999 & 1.000000000\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 1\n",
       "2. 0.00187299851992746\n",
       "3. 0.00187299851992746\n",
       "4. 1\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "      mu          sigma      \n",
       "mu    1.000000000 0.001872999\n",
       "sigma 0.001872999 1.000000000"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "diag(vcov(m4.1))\n",
    "cov2cor(vcov(m4.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ var = std^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two element vector (```diag(vcov(m4.1))```) in the output is the list of variances. If you take the square root, you get the standard deviation of the parameters, as shown in the precis output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<dl class=dl-horizontal>\n",
       "\t<dt>mu</dt>\n",
       "\t\t<dd>0.411922373966696</dd>\n",
       "\t<dt>sigma</dt>\n",
       "\t\t<dd>0.291258100940537</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[mu] 0.411922373966696\n",
       "\\item[sigma] 0.291258100940537\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "mu\n",
       ":   0.411922373966696sigma\n",
       ":   0.291258100940537\n",
       "\n"
      ],
      "text/plain": [
       "       mu     sigma \n",
       "0.4119224 0.2912581 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sqrt(diag(vcov(m4.1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>mu</th><th scope=col>sigma</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>154.4174</td><td>7.662421</td></tr>\n",
       "\t<tr><td>154.7229</td><td>7.843495</td></tr>\n",
       "\t<tr><td>154.6446</td><td>7.472021</td></tr>\n",
       "\t<tr><td>154.9366</td><td>7.919850</td></tr>\n",
       "\t<tr><td>153.8720</td><td>7.710660</td></tr>\n",
       "\t<tr><td>155.2253</td><td>7.981645</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       " mu & sigma\\\\\n",
       "\\hline\n",
       "\t 154.4174 & 7.662421\\\\\n",
       "\t 154.7229 & 7.843495\\\\\n",
       "\t 154.6446 & 7.472021\\\\\n",
       "\t 154.9366 & 7.919850\\\\\n",
       "\t 153.8720 & 7.710660\\\\\n",
       "\t 155.2253 & 7.981645\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "  mu       sigma   \n",
       "1 154.4174 7.662421\n",
       "2 154.7229 7.843495\n",
       "3 154.6446 7.472021\n",
       "4 154.9366 7.919850\n",
       "5 153.8720 7.710660\n",
       "6 155.2253 7.981645"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# R code 4.32\n",
    "post <- extract.samples(m4.1, n=1e4)\n",
    "head(post)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "        Mean StdDev  |0.89  0.89|\n",
       "mu    154.60   0.41 153.94 155.26\n",
       "sigma   7.73   0.29   7.26   8.19"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# R code 4.33\n",
    "precis(post)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Chapter 10: debugging PYMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data(UCBadmit)\n",
    "d <- UCBadmit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   dept applicant.gender admit reject applications\n",
      "1     A             male   512    313          825\n",
      "2     A           female    89     19          108\n",
      "3     B             male   353    207          560\n",
      "4     B           female    17      8           25\n",
      "5     C             male   120    205          325\n",
      "6     C           female   202    391          593\n",
      "7     D             male   138    279          417\n",
      "8     D           female   131    244          375\n",
      "9     E             male    53    138          191\n",
      "10    E           female    94    299          393\n",
      "11    F             male    22    351          373\n",
      "12    F           female    24    317          341\n"
     ]
    }
   ],
   "source": [
    "show(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d$male <- ifelse( d$applicant.gender == \"male\" , 1, 0 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m10.6 <- map (\n",
    "    alist(\n",
    "        admit ~ dbinom( applications, p),\n",
    "        logit(p) <- a + bm * male,\n",
    "        a ~ dnorm(0, 10),\n",
    "        bm ~ dnorm(0, 10)\n",
    "    ), \n",
    "    data = d \n",
    ")\n",
    "\n",
    "m10.7 <- map (\n",
    "    alist(\n",
    "        admit ~ dbinom( applications, p),\n",
    "        logit(p) <- a,\n",
    "        a ~ dnorm(0, 10)\n",
    "    ), \n",
    "    data = d \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "        WAIC pWAIC dWAIC weight    SE   dSE\n",
       "m10.6 5954.9   2.0   0.0      1 35.05    NA\n",
       "m10.7 6046.2   0.9  91.4      0 29.89 19.19"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "compare(m10.6, m10.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "    Mean StdDev  5.5% 94.5%\n",
       "a  -0.83   0.05 -0.91 -0.75\n",
       "bm  0.61   0.06  0.51  0.71"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "precis(m10.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}