notebooks/ch04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "Setup"
   },
   "outputs": [],
   "source": [
    "remotes::install_github('rmcelreath/rethinking', upgrade=F)\n",
    "library(rethinking)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.1"
   },
   "outputs": [],
   "source": [
    "pos <- replicate( 1000 , sum( runif(16,-1,1) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.2"
   },
   "outputs": [],
   "source": [
    "prod( 1 + runif(12,0,0.1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.3"
   },
   "outputs": [],
   "source": [
    "growth <- replicate( 10000 , prod( 1 + runif(12,0,0.1) ) )\n",
    "dens( growth , norm.comp=TRUE )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.4"
   },
   "outputs": [],
   "source": [
    "big <- replicate( 10000 , prod( 1 + runif(12,0,0.5) ) )\n",
    "small <- replicate( 10000 , prod( 1 + runif(12,0,0.01) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.5"
   },
   "outputs": [],
   "source": [
    "log.big <- replicate( 10000 , log(prod(1 + runif(12,0,0.5))) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.6"
   },
   "outputs": [],
   "source": [
    "w <- 6; n <- 9;\n",
    "p_grid <- seq(from=0,to=1,length.out=100)\n",
    "posterior <- dbinom(w,n,p_grid)*dunif(p_grid,0,1)\n",
    "posterior <- posterior/sum(posterior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.7"
   },
   "outputs": [],
   "source": [
    "library(rethinking)\n",
    "data(Howell1)\n",
    "d <- Howell1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.8"
   },
   "outputs": [],
   "source": [
    "str( d )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.9"
   },
   "outputs": [],
   "source": [
    "precis( d )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.10"
   },
   "outputs": [],
   "source": [
    "d$height"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.11"
   },
   "outputs": [],
   "source": [
    "d2 <- d[ d$age >= 18 , ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.12"
   },
   "outputs": [],
   "source": [
    "curve( dnorm( x , 178 , 20 ) , from=100 , to=250 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.13"
   },
   "outputs": [],
   "source": [
    "curve( dunif( x , 0 , 50 ) , from=-10 , to=60 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.14"
   },
   "outputs": [],
   "source": [
    "sample_mu <- rnorm( 1e4 , 178 , 20 )\n",
    "sample_sigma <- runif( 1e4 , 0 , 50 )\n",
    "prior_h <- rnorm( 1e4 , sample_mu , sample_sigma )\n",
    "dens( prior_h )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.15"
   },
   "outputs": [],
   "source": [
    "sample_mu <- rnorm( 1e4 , 178 , 100 )\n",
    "prior_h <- rnorm( 1e4 , sample_mu , sample_sigma )\n",
    "dens( prior_h )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.16"
   },
   "outputs": [],
   "source": [
    "mu.list <- seq( from=150, to=160 , length.out=100 )\n",
    "sigma.list <- seq( from=7 , to=9 , length.out=100 )\n",
    "post <- expand.grid( mu=mu.list , sigma=sigma.list )\n",
    "post$LL <- sapply( 1:nrow(post) , function(i) sum(\n",
    "    dnorm( d2$height , post$mu[i] , post$sigma[i] , log=TRUE ) ) )\n",
    "post$prod <- post$LL + dnorm( post$mu , 178 , 20 , TRUE ) +\n",
    "    dunif( post$sigma , 0 , 50 , TRUE )\n",
    "post$prob <- exp( post$prod - max(post$prod) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.17"
   },
   "outputs": [],
   "source": [
    "contour_xyz( post$mu , post$sigma , post$prob )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.18"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.18"
   },
   "outputs": [],
   "source": [
    "image_xyz( post$mu , post$sigma , post$prob )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.19"
   },
   "outputs": [],
   "source": [
    "sample.rows <- sample( 1:nrow(post) , size=1e4 , replace=TRUE ,\n",
    "    prob=post$prob )\n",
    "sample.mu <- post$mu[ sample.rows ]\n",
    "sample.sigma <- post$sigma[ sample.rows ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.20"
   },
   "outputs": [],
   "source": [
    "plot( sample.mu , sample.sigma , cex=0.5 , pch=16 , col=col.alpha(rangi2,0.1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.21"
   },
   "outputs": [],
   "source": [
    "dens( sample.mu )\n",
    "dens( sample.sigma )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.22"
   },
   "outputs": [],
   "source": [
    "PI( sample.mu )\n",
    "PI( sample.sigma )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.23"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.23"
   },
   "outputs": [],
   "source": [
    "d3 <- sample( d2$height , size=20 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.24"
   },
   "outputs": [],
   "source": [
    "mu.list <- seq( from=150, to=170 , length.out=200 )\n",
    "sigma.list <- seq( from=4 , to=20 , length.out=200 )\n",
    "post2 <- expand.grid( mu=mu.list , sigma=sigma.list )\n",
    "post2$LL <- sapply( 1:nrow(post2) , function(i)\n",
    "    sum( dnorm( d3 , mean=post2$mu[i] , sd=post2$sigma[i] ,\n",
    "    log=TRUE ) ) )\n",
    "post2$prod <- post2$LL + dnorm( post2$mu , 178 , 20 , TRUE ) +\n",
    "    dunif( post2$sigma , 0 , 50 , TRUE )\n",
    "post2$prob <- exp( post2$prod - max(post2$prod) )\n",
    "sample2.rows <- sample( 1:nrow(post2) , size=1e4 , replace=TRUE ,\n",
    "    prob=post2$prob )\n",
    "sample2.mu <- post2$mu[ sample2.rows ]\n",
    "sample2.sigma <- post2$sigma[ sample2.rows ]\n",
    "plot( sample2.mu , sample2.sigma , cex=0.5 ,\n",
    "    col=col.alpha(rangi2,0.1) ,\n",
    "    xlab=\"mu\" , ylab=\"sigma\" , pch=16 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.25"
   },
   "outputs": [],
   "source": [
    "dens( sample2.sigma , norm.comp=TRUE )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.26"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.26"
   },
   "outputs": [],
   "source": [
    "library(rethinking)\n",
    "data(Howell1)\n",
    "d <- Howell1\n",
    "d2 <- d[ d$age >= 18 , ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.27"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.27"
   },
   "outputs": [],
   "source": [
    "flist <- alist(\n",
    "    height ~ dnorm( mu , sigma ) ,\n",
    "    mu ~ dnorm( 178 , 20 ) ,\n",
    "    sigma ~ dunif( 0 , 50 )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.28"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.28"
   },
   "outputs": [],
   "source": [
    "m4.1 <- quap( flist , data=d2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.29"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.29"
   },
   "outputs": [],
   "source": [
    "precis( m4.1 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.30"
   },
   "outputs": [],
   "source": [
    "start <- list(\n",
    "    mu=mean(d2$height),\n",
    "    sigma=sd(d2$height)\n",
    ")\n",
    "m4.1 <- quap( flist , data=d2 , start=start )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.31"
   },
   "outputs": [],
   "source": [
    "m4.2 <- quap(\n",
    "    alist(\n",
    "        height ~ dnorm( mu , sigma ) ,\n",
    "        mu ~ dnorm( 178 , 0.1 ) ,\n",
    "        sigma ~ dunif( 0 , 50 )\n",
    "    ) , data=d2 )\n",
    "precis( m4.2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.32"
   },
   "outputs": [],
   "source": [
    "vcov( m4.1 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.33"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.33"
   },
   "outputs": [],
   "source": [
    "diag( vcov( m4.1 ) )\n",
    "cov2cor( vcov( m4.1 ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.34"
   },
   "outputs": [],
   "source": [
    "library(rethinking)\n",
    "post <- extract.samples( m4.1 , n=1e4 )\n",
    "head(post)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.35"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.35"
   },
   "outputs": [],
   "source": [
    "precis(post)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.36"
   },
   "outputs": [],
   "source": [
    "library(MASS)\n",
    "post <- mvrnorm( n=1e4 , mu=coef(m4.1) , Sigma=vcov(m4.1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.37"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.37"
   },
   "outputs": [],
   "source": [
    "library(rethinking)\n",
    "data(Howell1); d <- Howell1; d2 <- d[ d$age >= 18 , ]\n",
    "plot( d2$height ~ d2$weight )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.38"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.38"
   },
   "outputs": [],
   "source": [
    "set.seed(2971)\n",
    "N <- 100                   # 100 lines\n",
    "a <- rnorm( N , 178 , 20 )\n",
    "b <- rnorm( N , 0 , 10 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.39"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.39"
   },
   "outputs": [],
   "source": [
    "plot( NULL , xlim=range(d2$weight) , ylim=c(-100,400) ,\n",
    "    xlab=\"weight\" , ylab=\"height\" )\n",
    "abline( h=0 , lty=2 )\n",
    "abline( h=272 , lty=1 , lwd=0.5 )\n",
    "mtext( \"b ~ dnorm(0,10)\" )\n",
    "xbar <- mean(d2$weight)\n",
    "for ( i in 1:N ) curve( a[i] + b[i]*(x - xbar) ,\n",
    "    from=min(d2$weight) , to=max(d2$weight) , add=TRUE ,\n",
    "    col=col.alpha(\"black\",0.2) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.40"
   },
   "outputs": [],
   "source": [
    "b <- rlnorm( 1e4 , 0 , 1 )\n",
    "dens( b , xlim=c(0,5) , adj=0.1 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.41"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.41"
   },
   "outputs": [],
   "source": [
    "set.seed(2971)\n",
    "N <- 100                   # 100 lines\n",
    "a <- rnorm( N , 178 , 20 )\n",
    "b <- rlnorm( N , 0 , 1 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.42"
   },
   "outputs": [],
   "source": [
    "# load data again, since it's a long way back\n",
    "library(rethinking)\n",
    "data(Howell1); d <- Howell1; d2 <- d[ d$age >= 18 , ]\n",
    "\n",
    "# define the average weight, x-bar\n",
    "xbar <- mean(d2$weight)\n",
    "\n",
    "# fit model\n",
    "m4.3 <- quap(\n",
    "    alist(\n",
    "        height ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + b*( weight - xbar ) ,\n",
    "        a ~ dnorm( 178 , 20 ) ,\n",
    "        b ~ dlnorm( 0 , 1 ) ,\n",
    "        sigma ~ dunif( 0 , 50 )\n",
    "    ) , data=d2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.43"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.43"
   },
   "outputs": [],
   "source": [
    "m4.3b <- quap(\n",
    "    alist(\n",
    "        height ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + exp(log_b)*( weight - xbar ),\n",
    "        a ~ dnorm( 178 , 20 ) ,\n",
    "        log_b ~ dnorm( 0 , 1 ) ,\n",
    "        sigma ~ dunif( 0 , 50 )\n",
    "    ) , data=d2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.44"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.44"
   },
   "outputs": [],
   "source": [
    "precis( m4.3 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.45"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.45"
   },
   "outputs": [],
   "source": [
    "round( vcov( m4.3 ) , 3 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.46"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.46"
   },
   "outputs": [],
   "source": [
    "plot( height ~ weight , data=d2 , col=rangi2 )\n",
    "post <- extract.samples( m4.3 )\n",
    "a_map <- mean(post$a)\n",
    "b_map <- mean(post$b)\n",
    "curve( a_map + b_map*(x - xbar) , add=TRUE )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.47"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.47"
   },
   "outputs": [],
   "source": [
    "post <- extract.samples( m4.3 )\n",
    "post[1:5,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.48"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.48"
   },
   "outputs": [],
   "source": [
    "N <- 10\n",
    "dN <- d2[ 1:N , ]\n",
    "mN <- quap(\n",
    "    alist(\n",
    "        height ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + b*( weight - mean(weight) ) ,\n",
    "        a ~ dnorm( 178 , 20 ) ,\n",
    "        b ~ dlnorm( 0 , 1 ) ,\n",
    "        sigma ~ dunif( 0 , 50 )\n",
    "    ) , data=dN )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.49"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.49"
   },
   "outputs": [],
   "source": [
    "# extract 20 samples from the posterior\n",
    "post <- extract.samples( mN , n=20 )\n",
    "\n",
    "# display raw data and sample size\n",
    "plot( dN$weight , dN$height ,\n",
    "    xlim=range(d2$weight) , ylim=range(d2$height) ,\n",
    "    col=rangi2 , xlab=\"weight\" , ylab=\"height\" )\n",
    "mtext(concat(\"N = \",N))\n",
    "\n",
    "# plot the lines, with transparency\n",
    "for ( i in 1:20 )\n",
    "    curve( post$a[i] + post$b[i]*(x-mean(dN$weight)) ,\n",
    "        col=col.alpha(\"black\",0.3) , add=TRUE )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.50"
   },
   "outputs": [],
   "source": [
    "post <- extract.samples( m4.3 )\n",
    "mu_at_50 <- post$a + post$b * ( 50 - xbar )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.51"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.51"
   },
   "outputs": [],
   "source": [
    "dens( mu_at_50 , col=rangi2 , lwd=2 , xlab=\"mu|weight=50\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.52"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.52"
   },
   "outputs": [],
   "source": [
    "PI( mu_at_50 , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.53"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.53"
   },
   "outputs": [],
   "source": [
    "mu <- link( m4.3 )\n",
    "str(mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.54"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.54"
   },
   "outputs": [],
   "source": [
    "# define sequence of weights to compute predictions for\n",
    "# these values will be on the horizontal axis\n",
    "weight.seq <- seq( from=25 , to=70 , by=1 )\n",
    "\n",
    "# use link to compute mu\n",
    "# for each sample from posterior\n",
    "# and for each weight in weight.seq\n",
    "mu <- link( m4.3 , data=data.frame(weight=weight.seq) )\n",
    "str(mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.55"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.55"
   },
   "outputs": [],
   "source": [
    "# use type=\"n\" to hide raw data\n",
    "plot( height ~ weight , d2 , type=\"n\" )\n",
    "\n",
    "# loop over samples and plot each mu value\n",
    "for ( i in 1:100 )\n",
    "    points( weight.seq , mu[i,] , pch=16 , col=col.alpha(rangi2,0.1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.56"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.56"
   },
   "outputs": [],
   "source": [
    "# summarize the distribution of mu\n",
    "mu.mean <- apply( mu , 2 , mean )\n",
    "mu.PI <- apply( mu , 2 , PI , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.57"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.57"
   },
   "outputs": [],
   "source": [
    "# plot raw data\n",
    "# fading out points to make line and interval more visible\n",
    "plot( height ~ weight , data=d2 , col=col.alpha(rangi2,0.5) )\n",
    "\n",
    "# plot the MAP line, aka the mean mu for each weight\n",
    "lines( weight.seq , mu.mean )\n",
    "\n",
    "# plot a shaded region for 89% PI\n",
    "shade( mu.PI , weight.seq )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.58"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.58"
   },
   "outputs": [],
   "source": [
    "post <- extract.samples(m4.3)\n",
    "mu.link <- function(weight) post$a + post$b*( weight - xbar )\n",
    "weight.seq <- seq( from=25 , to=70 , by=1 )\n",
    "mu <- sapply( weight.seq , mu.link )\n",
    "mu.mean <- apply( mu , 2 , mean )\n",
    "mu.CI <- apply( mu , 2 , PI , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.59"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.59"
   },
   "outputs": [],
   "source": [
    "sim.height <- sim( m4.3 , data=list(weight=weight.seq) )\n",
    "str(sim.height)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.60"
   },
   "outputs": [],
   "source": [
    "height.PI <- apply( sim.height , 2 , PI , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.61"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.61"
   },
   "outputs": [],
   "source": [
    "# plot raw data\n",
    "plot( height ~ weight , d2 , col=col.alpha(rangi2,0.5) )\n",
    "\n",
    "# draw MAP line\n",
    "lines( weight.seq , mu.mean )\n",
    "\n",
    "# draw HPDI region for line\n",
    "shade( mu.HPDI , weight.seq )\n",
    "\n",
    "# draw PI region for simulated heights\n",
    "shade( height.PI , weight.seq )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.62"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.62"
   },
   "outputs": [],
   "source": [
    "sim.height <- sim( m4.3 , data=list(weight=weight.seq) , n=1e4 )\n",
    "height.PI <- apply( sim.height , 2 , PI , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.63"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.63"
   },
   "outputs": [],
   "source": [
    "post <- extract.samples(m4.3)\n",
    "weight.seq <- 25:70\n",
    "sim.height <- sapply( weight.seq , function(weight)\n",
    "    rnorm(\n",
    "        n=nrow(post) ,\n",
    "        mean=post$a + post$b*( weight - xbar ) ,\n",
    "        sd=post$sigma ) )\n",
    "height.PI <- apply( sim.height , 2 , PI , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.64"
   },
   "outputs": [],
   "source": [
    "library(rethinking)\n",
    "data(Howell1)\n",
    "d <- Howell1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.65"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.65"
   },
   "outputs": [],
   "source": [
    "d$weight_s <- ( d$weight - mean(d$weight) )/sd(d$weight)\n",
    "d$weight_s2 <- d$weight_s^2\n",
    "m4.5 <- quap(\n",
    "    alist(\n",
    "        height ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + b1*weight_s + b2*weight_s2 ,\n",
    "        a ~ dnorm( 178 , 20 ) ,\n",
    "        b1 ~ dlnorm( 0 , 1 ) ,\n",
    "        b2 ~ dnorm( 0 , 1 ) ,\n",
    "        sigma ~ dunif( 0 , 50 )\n",
    "    ) , data=d )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.66"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.66"
   },
   "outputs": [],
   "source": [
    "precis( m4.5 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.67"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.67"
   },
   "outputs": [],
   "source": [
    "weight.seq <- seq( from=-2.2 , to=2 , length.out=30 )\n",
    "pred_dat <- list( weight_s=weight.seq , weight_s2=weight.seq^2 )\n",
    "mu <- link( m4.5 , data=pred_dat )\n",
    "mu.mean <- apply( mu , 2 , mean )\n",
    "mu.PI <- apply( mu , 2 , PI , prob=0.89 )\n",
    "sim.height <- sim( m4.5 , data=pred_dat )\n",
    "height.PI <- apply( sim.height , 2 , PI , prob=0.89 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.68"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.68"
   },
   "outputs": [],
   "source": [
    "plot( height ~ weight_s , d , col=col.alpha(rangi2,0.5) )\n",
    "lines( weight.seq , mu.mean )\n",
    "shade( mu.PI , weight.seq )\n",
    "shade( height.PI , weight.seq )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.69"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.69"
   },
   "outputs": [],
   "source": [
    "d$weight_s3 <- d$weight_s^3\n",
    "m4.6 <- quap(\n",
    "    alist(\n",
    "        height ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + b1*weight_s + b2*weight_s2 + b3*weight_s3 ,\n",
    "        a ~ dnorm( 178 , 20 ) ,\n",
    "        b1 ~ dlnorm( 0 , 1 ) ,\n",
    "        b2 ~ dnorm( 0 , 10 ) ,\n",
    "        b3 ~ dnorm( 0 , 10 ) ,\n",
    "        sigma ~ dunif( 0 , 50 )\n",
    "    ) , data=d )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.70"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.70"
   },
   "outputs": [],
   "source": [
    "plot( height ~ weight_s , d , col=col.alpha(rangi2,0.5) , xaxt=\"n\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.71"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.71"
   },
   "outputs": [],
   "source": [
    "at <- c(-2,-1,0,1,2)\n",
    "labels <- at*sd(d$weight) + mean(d$weight)\n",
    "axis( side=1 , at=at , labels=round(labels,1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.72"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.72"
   },
   "outputs": [],
   "source": [
    "library(rethinking)\n",
    "data(cherry_blossoms)\n",
    "d <- cherry_blossoms\n",
    "precis(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.73"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.73"
   },
   "outputs": [],
   "source": [
    "d2 <- d[ complete.cases(d$doy) , ] # complete cases on doy\n",
    "num_knots <- 15\n",
    "knot_list <- quantile( d2$year , probs=seq(0,1,length.out=num_knots) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.74"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.74"
   },
   "outputs": [],
   "source": [
    "library(splines)\n",
    "B <- bs(d2$year,\n",
    "    knots=knot_list[-c(1,num_knots)] ,\n",
    "    degree=3 , intercept=TRUE )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.75"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.75"
   },
   "outputs": [],
   "source": [
    "plot( NULL , xlim=range(d2$year) , ylim=c(0,1) , xlab=\"year\" , ylab=\"basis\" )\n",
    "for ( i in 1:ncol(B) ) lines( d2$year , B[,i] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.76"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.76"
   },
   "outputs": [],
   "source": [
    "m4.7 <- quap(\n",
    "    alist(\n",
    "        D ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + B %*% w ,\n",
    "        a ~ dnorm(100,10),\n",
    "        w ~ dnorm(0,10),\n",
    "        sigma ~ dexp(1)\n",
    "    ), data=list( D=d2$doy , B=B ) ,\n",
    "    start=list( w=rep( 0 , ncol(B) ) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.77"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.77"
   },
   "outputs": [],
   "source": [
    "post <- extract.samples( m4.7 )\n",
    "w <- apply( post$w , 2 , mean )\n",
    "plot( NULL , xlim=range(d2$year) , ylim=c(-6,6) ,\n",
    "    xlab=\"year\" , ylab=\"basis * weight\" )\n",
    "for ( i in 1:ncol(B) ) lines( d2$year , w[i]*B[,i] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.78"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "name": "R code 4.78"
   },
   "outputs": [],
   "source": [
    "mu <- link( m4.7 )\n",
    "mu_PI <- apply(mu,2,PI,0.97)\n",
    "plot( d2$year , d2$doy , col=col.alpha(rangi2,0.3) , pch=16 )\n",
    "shade( mu_PI , d2$year , col=col.alpha(\"black\",0.5) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R code 4.79"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "R code 4.79"
   },
   "outputs": [],
   "source": [
    "m4.7alt <- quap(\n",
    "    alist(\n",
    "        D ~ dnorm( mu , sigma ) ,\n",
    "        mu <- a + sapply( 1:827 , function(i) sum( B[i,]*w ) ) ,\n",
    "        a ~ dnorm(100,1),\n",
    "        w ~ dnorm(0,10),\n",
    "        sigma ~ dexp(1)\n",
    "    ),\n",
    "    data=list( D=d2$doy , B=B ) ,\n",
    "    start=list( w=rep( 0 , ncol(B) ) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "name,-all",
   "main_language": "R",
   "notebook_metadata_filter": "-all"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}