## R code 4.1
library(rethinking)
## Loading required package: rstan
## Loading required package: Rcpp
## Loading required package: ggplot2
## rstan (Version 2.8.0, packaged: 2015-09-19 14:48:38 UTC, GitRev: 05c3d0058b6a)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
## Loading required package: parallel
## rethinking (Version 1.57)
pos <- replicate( 1000 , sum( runif(16,-1,1) ) )
## R code 4.2
prod( 1 + runif(12,0,0.1) )
## [1] 1.778027
## R code 4.3
growth <- replicate( 10000 , prod( 1 + runif(12,0,0.1) ) )
dens( growth , norm.comp=TRUE )
## R code 4.4
big <- replicate( 10000 , prod( 1 + runif(12,0,0.5) ) )
small <- replicate( 10000 , prod( 1 + runif(12,0,0.01) ) )
## R code 4.5
log.big <- replicate( 10000 , log(prod(1 + runif(12,0,0.5))) )
## R code 4.6
w <- 6; n <- 9;
p_grid <- seq(from=0,to=1,length.out=100)
posterior <- dbinom(w,n,p_grid)*dunif(p_grid,0,1)
posterior <- posterior/sum(posterior)
## R code 4.7
library(rethinking)
data(Howell1)
d <- Howell1
## R code 4.8
str( d )
## 'data.frame': 544 obs. of 4 variables:
## $ height: num 152 140 137 157 145 ...
## $ weight: num 47.8 36.5 31.9 53 41.3 ...
## $ age : num 63 63 65 41 51 35 32 27 19 54 ...
## $ male : int 1 0 0 1 0 1 0 1 0 1 ...
## R code 4.9
d$height
## [1] 151.7650 139.7000 136.5250 156.8450 145.4150 163.8300 149.2250
## [8] 168.9100 147.9550 165.1000 154.3050 151.1300 144.7800 149.9000
## [15] 150.4950 163.1950 157.4800 143.9418 121.9200 105.4100 86.3600
## [22] 161.2900 156.2100 129.5400 109.2200 146.4000 148.5900 147.3200
## [29] 137.1600 125.7300 114.3000 147.9550 161.9250 146.0500 146.0500
## [36] 152.7048 142.8750 142.8750 147.9550 160.6550 151.7650 162.8648
## [43] 171.4500 147.3200 147.9550 144.7800 121.9200 128.9050 97.7900
## [50] 154.3050 143.5100 146.7000 157.4800 127.0000 110.4900 97.7900
## [57] 165.7350 152.4000 141.6050 158.8000 155.5750 164.4650 151.7650
## [64] 161.2900 154.3050 145.4150 145.4150 152.4000 163.8300 144.1450
## [71] 129.5400 129.5400 153.6700 142.8750 146.0500 167.0050 158.4198
## [78] 91.4400 165.7350 149.8600 147.9550 137.7950 154.9400 160.9598
## [85] 161.9250 147.9550 113.6650 159.3850 148.5900 136.5250 158.1150
## [92] 144.7800 156.8450 179.0700 118.7450 170.1800 146.0500 147.3200
## [99] 113.0300 162.5600 133.9850 152.4000 160.0200 149.8600 142.8750
## [106] 167.0050 159.3850 154.9400 148.5900 111.1250 111.7600 162.5600
## [113] 152.4000 124.4600 111.7600 86.3600 170.1800 146.0500 159.3850
## [120] 151.1300 160.6550 169.5450 158.7500 74.2950 149.8600 153.0350
## [127] 96.5200 161.9250 162.5600 149.2250 116.8400 100.0760 163.1950
## [134] 161.9250 145.4150 163.1950 151.1300 150.4950 141.6050 170.8150
## [141] 91.4400 157.4800 152.4000 149.2250 129.5400 147.3200 145.4150
## [148] 121.9200 113.6650 157.4800 154.3050 120.6500 115.6000 167.0050
## [155] 142.8750 152.4000 96.5200 160.0000 159.3850 149.8600 160.6550
## [162] 160.6550 149.2250 125.0950 140.9700 154.9400 141.6050 160.0200
## [169] 150.1648 155.5750 103.5050 94.6150 156.2100 153.0350 167.0050
## [176] 149.8600 147.9550 159.3850 161.9250 155.5750 159.3850 146.6850
## [183] 172.7200 166.3700 141.6050 142.8750 133.3500 127.6350 119.3800
## [190] 151.7650 156.8450 148.5900 157.4800 149.8600 147.9550 102.2350
## [197] 153.0350 160.6550 149.2250 114.3000 100.9650 138.4300 91.4400
## [204] 162.5600 149.2250 158.7500 149.8600 158.1150 156.2100 148.5900
## [211] 143.5100 154.3050 131.4450 157.4800 157.4800 154.3050 107.9500
## [218] 168.2750 145.4150 147.9550 100.9650 113.0300 149.2250 154.9400
## [225] 162.5600 156.8450 123.1900 161.0106 144.7800 143.5100 149.2250
## [232] 110.4900 149.8600 165.7350 144.1450 157.4800 154.3050 163.8300
## [239] 156.2100 153.6700 134.6200 144.1450 114.3000 162.5600 146.0500
## [246] 120.6500 154.9400 144.7800 106.6800 146.6850 152.4000 163.8300
## [253] 165.7350 156.2100 152.4000 140.3350 158.1150 163.1950 151.1300
## [260] 171.1198 149.8600 163.8300 141.6050 93.9800 149.2250 105.4100
## [267] 146.0500 161.2900 162.5600 145.4150 145.4150 170.8150 127.0000
## [274] 159.3850 159.4000 153.6700 160.0200 150.4950 149.2250 127.0000
## [281] 142.8750 142.1130 147.3200 162.5600 164.4650 160.0200 153.6700
## [288] 167.0050 151.1300 147.9550 125.3998 111.1250 153.0350 139.0650
## [295] 152.4000 154.9400 147.9550 143.5100 117.9830 144.1450 92.7100
## [302] 147.9550 155.5750 150.4950 155.5750 154.3050 130.6068 101.6000
## [309] 157.4800 168.9100 150.4950 111.7600 160.0200 167.6400 144.1450
## [316] 145.4150 160.0200 147.3200 164.4650 153.0350 149.2250 160.0200
## [323] 149.2250 85.0900 84.4550 59.6138 92.7100 111.1250 90.8050
## [330] 153.6700 99.6950 62.4840 81.9150 96.5200 80.0100 150.4950
## [337] 151.7650 140.6398 88.2650 158.1150 149.2250 151.7650 154.9400
## [344] 123.8250 104.1400 161.2900 148.5900 97.1550 93.3450 160.6550
## [351] 157.4800 167.0050 157.4800 91.4400 60.4520 137.1600 152.4000
## [358] 152.4000 81.2800 109.2200 71.1200 89.2048 67.3100 85.0900
## [365] 69.8500 161.9250 152.4000 88.9000 90.1700 71.7550 83.8200
## [372] 159.3850 142.2400 142.2400 168.9100 123.1900 74.9300 74.2950
## [379] 90.8050 160.0200 67.9450 135.8900 158.1150 85.0900 93.3450
## [386] 152.4000 155.5750 154.3050 156.8450 120.0150 114.3000 83.8200
## [393] 156.2100 137.1600 114.3000 93.9800 168.2750 147.9550 139.7000
## [400] 157.4800 76.2000 66.0400 160.7000 114.3000 146.0500 161.2900
## [407] 69.8500 133.9850 67.9450 150.4950 163.1950 148.5900 148.5900
## [414] 161.9250 153.6700 68.5800 151.1300 163.8300 153.0350 151.7650
## [421] 132.0800 156.2100 140.3350 158.7500 142.8750 84.4550 151.9428
## [428] 161.2900 127.9906 160.9852 144.7800 132.0800 117.9830 160.0200
## [435] 154.9400 160.9852 165.9890 157.9880 154.9400 97.9932 64.1350
## [442] 160.6550 147.3200 146.7000 147.3200 172.9994 158.1150 147.3200
## [449] 124.9934 106.0450 165.9890 149.8600 76.2000 161.9250 140.0048
## [456] 66.6750 62.8650 163.8300 147.9550 160.0200 154.9400 152.4000
## [463] 62.2300 146.0500 151.9936 157.4800 55.8800 60.9600 151.7650
## [470] 144.7800 118.1100 78.1050 160.6550 151.1300 121.9200 92.7100
## [477] 153.6700 147.3200 139.7000 157.4800 91.4400 154.9400 143.5100
## [484] 83.1850 158.1150 147.3200 123.8250 88.9000 160.0200 137.1600
## [491] 165.1000 154.9400 111.1250 153.6700 145.4150 141.6050 144.7800
## [498] 163.8300 161.2900 154.9000 161.3000 170.1800 149.8600 123.8250
## [505] 85.0900 160.6550 154.9400 106.0450 126.3650 166.3700 148.2852
## [512] 124.4600 89.5350 101.6000 151.7650 148.5900 153.6700 53.9750
## [519] 146.6850 56.5150 100.9650 121.9200 81.5848 154.9400 156.2100
## [526] 132.7150 125.0950 101.6000 160.6550 146.0500 132.7150 87.6300
## [533] 156.2100 152.4000 162.5600 114.9350 67.9450 142.8750 76.8350
## [540] 145.4150 162.5600 156.2100 71.1200 158.7500
## R code 4.10
d2 <- d[ d$age >= 18 , ]
## R code 4.11
curve( dnorm( x , 178 , 20 ) , from=100 , to=250 )
## R code 4.12
curve( dunif( x , 0 , 50 ) , from=-10 , to=60 )
## R code 4.13
sample_mu <- rnorm( 1e4 , 178 , 20 )
sample_sigma <- runif( 1e4 , 0 , 50 )
prior_h <- rnorm( 1e4 , sample_mu , sample_sigma )
dens( prior_h )
## R code 4.14
mu.list <- seq( from=140, to=160 , length.out=200 )
sigma.list <- seq( from=4 , to=9 , length.out=200 )
post <- expand.grid( mu=mu.list , sigma=sigma.list )
post$LL <- sapply( 1:nrow(post) , function(i) sum( dnorm(
d2$height ,
mean=post$mu[i] ,
sd=post$sigma[i] ,
log=TRUE ) ) )
post$prod <- post$LL + dnorm( post$mu , 178 , 20 , TRUE ) +
dunif( post$sigma , 0 , 50 , TRUE )
post$prob <- exp( post$prod - max(post$prod) )
## R code 4.15
contour_xyz( post$mu , post$sigma , post$prob )
## R code 4.16
image_xyz( post$mu , post$sigma , post$prob )
## R code 4.17
sample.rows <- sample( 1:nrow(post) , size=1e4 , replace=TRUE ,
prob=post$prob )
sample.mu <- post$mu[ sample.rows ]
sample.sigma <- post$sigma[ sample.rows ]
## R code 4.18
plot( sample.mu , sample.sigma , cex=0.5 , pch=16 ,
col=col.alpha(rangi2,0.1) )
## R code 4.19
dens( sample.mu )
dens( sample.sigma )
## R code 4.20
HPDI( sample.mu )
## |0.89 0.89|
## 153.8693 155.1759
HPDI( sample.sigma )
## |0.89 0.89|
## 7.266332 8.195980
## R code 4.21
d3 <- sample( d2$height , size=20 )
## R code 4.22
mu.list <- seq( from=150, to=170 , length.out=200 )
sigma.list <- seq( from=4 , to=20 , length.out=200 )
post2 <- expand.grid( mu=mu.list , sigma=sigma.list )
post2$LL <- sapply( 1:nrow(post2) , function(i)
sum( dnorm( d3 , mean=post2$mu[i] , sd=post2$sigma[i] ,
log=TRUE ) ) )
post2$prod <- post2$LL + dnorm( post2$mu , 178 , 20 , TRUE ) +
dunif( post2$sigma , 0 , 50 , TRUE )
post2$prob <- exp( post2$prod - max(post2$prod) )
sample2.rows <- sample( 1:nrow(post2) , size=1e4 , replace=TRUE ,
prob=post2$prob )
sample2.mu <- post2$mu[ sample2.rows ]
sample2.sigma <- post2$sigma[ sample2.rows ]
plot( sample2.mu , sample2.sigma , cex=0.5 ,
col=col.alpha(rangi2,0.1) ,
xlab="mu" , ylab="sigma" , pch=16 )
## R code 4.23
dens( sample2.sigma , norm.comp=TRUE )
## R code 4.24
library(rethinking)
data(Howell1)
d <- Howell1
d2 <- d[ d$age >= 18 , ]
head(d2)
## height weight age male
## 1 151.765 47.82561 63 1
## 2 139.700 36.48581 63 0
## 3 136.525 31.86484 65 0
## 4 156.845 53.04191 41 1
## 5 145.415 41.27687 51 0
## 6 163.830 62.99259 35 1
## R code 4.25
flist <- alist(
height ~ dnorm( mu , sigma ) ,
mu ~ dnorm( 178 , 20 ) ,
sigma ~ dunif( 0 , 50 )
)
## R code 4.26
m4.1 <- map( flist , data=d2 )
## R code 4.27
precis( m4.1 )
## Mean StdDev 5.5% 94.5%
## mu 154.61 0.41 153.95 155.27
## sigma 7.73 0.29 7.27 8.20
## R code 4.28
start <- list(
mu=mean(d2$height),
sigma=sd(d2$height)
)
## R code 4.29
m4.2 <- map(
alist(
height ~ dnorm( mu , sigma ) ,
mu ~ dnorm( 178 , 0.1 ) ,
sigma ~ dunif( 0 , 50 )
) ,
data=d2 )
precis( m4.2 )
## Mean StdDev 5.5% 94.5%
## mu 177.86 0.10 177.70 178.02
## sigma 24.52 0.93 23.03 26.00
## R code 4.30
vcov( m4.1 )
## mu sigma
## mu 0.169739634 0.000218089
## sigma 0.000218089 0.084905854
## R code 4.31
diag( vcov( m4.1 ) )
## mu sigma
## 0.16973963 0.08490585
cov2cor( vcov( m4.1 ) )
## mu sigma
## mu 1.000000000 0.001816659
## sigma 0.001816659 1.000000000
## R code 4.32
library(rethinking)
post <- extract.samples( m4.1 , n=1e4 )
head(post)
## mu sigma
## 1 154.6518 7.998280
## 2 154.0942 7.083426
## 3 154.3286 7.237438
## 4 155.0338 8.015964
## 5 154.4840 7.852386
## 6 154.7178 7.766070
## R code 4.33
precis(post)
## Mean StdDev |0.89 0.89|
## mu 154.61 0.41 153.91 155.23
## sigma 7.73 0.29 7.24 8.17
## R code 4.34
library(MASS)
post <- mvrnorm( n=1e4 , mu=coef(m4.1) , Sigma=vcov(m4.1) )
## R code 4.35
m4.1_logsigma <- map(
alist(
height ~ dnorm( mu , exp(log_sigma) ) ,
mu ~ dnorm( 178 , 20 ) ,
log_sigma ~ dnorm( 2 , 10 )
) , data=d2 )
## R code 4.36
post <- extract.samples( m4.1_logsigma )
sigma <- exp( post$log_sigma )
## R code 4.37
plot( d2$height ~ d2$weight )
## R code 4.38
# load data again, since it's a long way back
library(rethinking)
data(Howell1)
d <- Howell1
d2 <- d[ d$age >= 18 , ]
head(d2)
## height weight age male
## 1 151.765 47.82561 63 1
## 2 139.700 36.48581 63 0
## 3 136.525 31.86484 65 0
## 4 156.845 53.04191 41 1
## 5 145.415 41.27687 51 0
## 6 163.830 62.99259 35 1
summary(d2)
## height weight age male
## Min. :136.5 Min. :31.07 Min. :18.00 Min. :0.0000
## 1st Qu.:148.6 1st Qu.:40.26 1st Qu.:28.00 1st Qu.:0.0000
## Median :154.3 Median :44.79 Median :39.00 Median :0.0000
## Mean :154.6 Mean :44.99 Mean :41.14 Mean :0.4688
## 3rd Qu.:160.7 3rd Qu.:49.29 3rd Qu.:51.00 3rd Qu.:1.0000
## Max. :179.1 Max. :62.99 Max. :88.00 Max. :1.0000
dim(d2)
## [1] 352 4
# fit model
m4.3 <- map(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + b*weight ,
a ~ dnorm( 156 , 100 ) ,
b ~ dnorm( 0 , 10 ) ,
sigma ~ dunif( 0 , 50 )
) ,
data=d2 )
## R code 4.39
m4.3 <- map(
alist(
height ~ dnorm( a + b*weight , sigma ) ,
a ~ dnorm( 178 , 100 ) ,
b ~ dnorm( 0 , 10 ) ,
sigma ~ dunif( 0 , 50 )
) ,
data=d2 )
## R code 4.40
precis( m4.3 )
## Mean StdDev 5.5% 94.5%
## a 113.90 1.91 110.86 116.95
## b 0.90 0.04 0.84 0.97
## sigma 5.07 0.19 4.77 5.38
## R code 4.41
precis( m4.3 , corr=TRUE )
## Mean StdDev 5.5% 94.5% a b sigma
## a 113.90 1.91 110.86 116.95 1.00 -0.99 0
## b 0.90 0.04 0.84 0.97 -0.99 1.00 0
## sigma 5.07 0.19 4.77 5.38 0.00 0.00 1
## R code 4.42
d2$weight.c <- d2$weight - mean(d2$weight)
## R code 4.43
m4.4 <- map(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + b*weight.c ,
a ~ dnorm( 178 , 100 ) ,
b ~ dnorm( 0 , 10 ) ,
sigma ~ dunif( 0 , 50 )
) ,
data=d2 )
## R code 4.44
precis( m4.4 , corr=TRUE )
## Mean StdDev 5.5% 94.5% a b sigma
## a 154.60 0.27 154.17 155.03 1 0 0
## b 0.91 0.04 0.84 0.97 0 1 0
## sigma 5.07 0.19 4.77 5.38 0 0 1
## R code 4.45
plot( height ~ weight , data=d2 )
abline( a=coef(m4.3)["a"] , b=coef(m4.3)["b"] )
## R code 4.46
post <- extract.samples( m4.3 )
## R code 4.47
post[1:5,]
## a b sigma
## 1 113.1740 0.9226801 5.022161
## 2 116.2994 0.8548417 5.019622
## 3 110.9260 0.9684564 5.544658
## 4 113.4436 0.9021136 4.893569
## 5 111.9005 0.9514062 4.793823
## R code 4.48
N <- 10
dN <- d2[ 1:N , ]
mN <- map(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + b*weight ,
a ~ dnorm( 178 , 100 ) ,
b ~ dnorm( 0 , 10 ) ,
sigma ~ dunif( 0 , 50 )
) , data=dN )
## R code 4.49
# extract 20 samples from the posterior
post <- extract.samples( mN , n=20 )
# display raw data and sample size
plot( dN$weight , dN$height ,
xlim=range(d2$weight) , ylim=range(d2$height) ,
col=rangi2 , xlab="weight" , ylab="height" )
mtext(concat("N = ",N))
# plot the lines, with transparency
for ( i in 1:20 )
abline( a=post$a[i] , b=post$b[i] , col=col.alpha("black",0.3) )
## R code 4.50
mu_at_50 <- post$a + post$b * 50
## R code 4.51
dens( mu_at_50 , col=rangi2 , lwd=2 , xlab="mu|weight=50" )
## R code 4.52
HPDI( mu_at_50 , prob=0.89 )
## |0.89 0.89|
## 153.6535 159.4628
## R code 4.53
mu <- link( m4.3 )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
str(mu)
## num [1:1000, 1:352] 44.7 43.3 44.2 45 45.5 ...
## R code 4.54
# define sequence of weights to compute predictions for
# these values will be on the horizontal axis
weight.seq <- seq( from=25 , to=70 , by=1 )
# use link to compute mu
# for each sample from posterior
# and for each weight in weight.seq
mu <- link( m4.3 , data=data.frame(weight=weight.seq) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
str(mu)
## num [1:1000, 1:46] 23.1 23.7 20.3 21.3 23 ...
## R code 4.55
# use type="n" to hide raw data
plot( height ~ weight , d2 , type="n" )
# loop over samples and plot each mu value
for ( i in 1:100 )
points( weight.seq , mu[i,] , pch=16 , col=col.alpha(rangi2,0.1) )
## R code 4.56
# summarize the distribution of mu
mu.mean <- apply( mu , 2 , mean )
mu.HPDI <- apply( mu , 2 , HPDI , prob=0.89 )
## R code 4.57
# plot raw data
# fading out points to make line and interval more visible
plot( height ~ weight , data=d2 , col=col.alpha(rangi2,0.5) )
# plot the MAP line, aka the mean mu for each weight
lines( weight.seq , mu.mean )
# plot a shaded region for 89% HPDI
shade( mu.HPDI , weight.seq )
## R code 4.58
post <- extract.samples(m4.3)
mu.link <- function(weight) post$a + post$b*weight
weight.seq <- seq( from=25 , to=70 , by=1 )
mu <- sapply( weight.seq , mu.link )
mu.mean <- apply( mu , 2 , mean )
mu.HPDI <- apply( mu , 2 , HPDI , prob=0.89 )
## R code 4.59
sim.height <- sim( m4.3 , data=list(weight=weight.seq) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
str(sim.height)
## num [1:1000, 1:46] 136 142 133 139 149 ...
## R code 4.60
height.PI <- apply( sim.height , 2 , PI , prob=0.89 )
## R code 4.61
# plot raw data
plot( height ~ weight , d2 , col=col.alpha(rangi2,0.5) )
# draw MAP line
lines( weight.seq , mu.mean )
# draw HPDI region for line
shade( mu.HPDI , weight.seq )
# draw PI region for simulated heights
shade( height.PI , weight.seq )
## R code 4.62
sim.height <- sim( m4.3 , data=list(weight=weight.seq) , n=1e4 )
## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]
height.PI <- apply( sim.height , 2 , PI , prob=0.89 )
## R code 4.63
post <- extract.samples(m4.3)
weight.seq <- 25:70
sim.height <- sapply( weight.seq , function(weight)
rnorm(
n=nrow(post) ,
mean=post$a + post$b*weight ,
sd=post$sigma ) )
height.PI <- apply( sim.height , 2 , PI , prob=0.89 )
## R code 4.64
library(rethinking)
data(Howell1)
d <- Howell1
str(d)
## 'data.frame': 544 obs. of 4 variables:
## $ height: num 152 140 137 157 145 ...
## $ weight: num 47.8 36.5 31.9 53 41.3 ...
## $ age : num 63 63 65 41 51 35 32 27 19 54 ...
## $ male : int 1 0 0 1 0 1 0 1 0 1 ...
## R code 4.65
d$weight.s <- ( d$weight - mean(d$weight) )/sd(d$weight)
## R code 4.66
d$weight.s2 <- d$weight.s^2
m4.5 <- map(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + b1*weight.s + b2*weight.s2 ,
a ~ dnorm( 178 , 100 ) ,
b1 ~ dnorm( 0 , 10 ) ,
b2 ~ dnorm( 0 , 10 ) ,
sigma ~ dunif( 0 , 50 )
) ,
data=d )
## R code 4.67
precis( m4.5 )
## Mean StdDev 5.5% 94.5%
## a 146.66 0.37 146.07 147.26
## b1 21.40 0.29 20.94 21.86
## b2 -8.42 0.28 -8.86 -7.97
## sigma 5.75 0.17 5.47 6.03
## R code 4.68
weight.seq <- seq( from=-2.2 , to=2 , length.out=30 )
pred_dat <- list( weight.s=weight.seq , weight.s2=weight.seq^2 )
mu <- link( m4.5 , data=pred_dat )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.mean <- apply( mu , 2 , mean )
mu.PI <- apply( mu , 2 , PI , prob=0.89 )
sim.height <- sim( m4.5 , data=pred_dat )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
height.PI <- apply( sim.height , 2 , PI , prob=0.89 )
## R code 4.69
plot( height ~ weight.s , d , col=col.alpha(rangi2,0.5) )
lines( weight.seq , mu.mean )
shade( mu.PI , weight.seq )
shade( height.PI , weight.seq )
## R code 4.70
d$weight.s3 <- d$weight.s^3
m4.6 <- map(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + b1*weight.s + b2*weight.s2 + b3*weight.s3 ,
a ~ dnorm( 178 , 100 ) ,
b1 ~ dnorm( 0 , 10 ) ,
b2 ~ dnorm( 0 , 10 ) ,
b3 ~ dnorm( 0 , 10 ) ,
sigma ~ dunif( 0 , 50 )
) ,
data=d )
## R code 4.71
plot( height ~ weight.s , d , col=col.alpha(rangi2,0.5) , xaxt="n" )
## R code 4.72
at <- c(-2,-1,0,1,2)
labels <- at*sd(d$weight) + mean(d$weight)
axis( side=1 , at=at , labels=round(labels,1) )
## R code 4.73
plot( height ~ weight , data=Howell1 ,
col=col.alpha(rangi2,0.4) )
