## R code 14.1
# simulate a pancake and return randomly ordered sides
sim_pancake <- function() {
pancake <- sample(1:3,1)
sides <- matrix(c(1,1,1,0,0,0),2,3)[,pancake]
sample(sides)
}
# sim 10,000 pancakes
pancakes <- replicate( 1e4 , sim_pancake() )
up <- pancakes[1,]
down <- pancakes[2,]
# compute proportion 1/1 (BB) out of all 1/1 and 1/0
num_11_10 <- sum( up==1 )
num_11 <- sum( up==1 & down==1 )
num_11/num_11_10
## [1] 0.6660606
## R code 14.2
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)
data(WaffleDivorce)
d <- WaffleDivorce
# points
plot( d$Divorce ~ d$MedianAgeMarriage , ylim=c(4,15) ,
xlab="Median age marriage" , ylab="Divorce rate" )
# standard errors
for ( i in 1:nrow(d) ) {
ci <- d$Divorce[i] + c(-1,1)*d$Divorce.SE[i]
x <- d$MedianAgeMarriage[i]
lines( c(x,x) , ci )
}
## R code 14.3
dlist <- list(
div_obs=d$Divorce,
div_sd=d$Divorce.SE,
R=d$Marriage,
A=d$MedianAgeMarriage
)
m14.1 <- map2stan(
alist(
div_est ~ dnorm(mu,sigma),
mu <- a + bA*A + bR*R,
div_obs ~ dnorm(div_est,div_sd),
a ~ dnorm(0,10),
bA ~ dnorm(0,10),
bR ~ dnorm(0,10),
sigma ~ dcauchy(0,2.5)
) ,
data=dlist ,
start=list(div_est=dlist$div_obs) ,
WAIC=FALSE , iter=5000 , warmup=1000 , chains=2 , cores=2 ,
control=list(adapt_delta=0.95) )
## Precompiling model...
## the number of chains is less than 1; sampling not done
## Dispatching 2 chains to 2 cores.
##
## SAMPLING FOR MODEL 'div_est ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## R code 14.4
precis( m14.1 , depth=2 )
## Warning: 使われていないコネクション 6 (<-kuma:11429) を閉じます
## Warning: 使われていないコネクション 5 (<-kuma:11429) を閉じます
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## div_est[1] 11.78 0.67 10.70 12.84 8000 1
## div_est[2] 11.19 1.02 9.65 12.90 8000 1
## div_est[3] 10.48 0.62 9.50 11.48 8000 1
## div_est[4] 12.30 0.87 10.94 13.70 8000 1
## div_est[5] 8.05 0.24 7.66 8.41 8000 1
## div_est[6] 11.01 0.74 9.83 12.18 8000 1
## div_est[7] 7.24 0.64 6.22 8.25 8000 1
## div_est[8] 9.36 0.90 7.94 10.81 8000 1
## div_est[9] 7.01 1.09 5.22 8.72 8000 1
## div_est[10] 8.54 0.31 8.03 9.03 8000 1
## div_est[11] 11.15 0.53 10.31 11.99 8000 1
## div_est[12] 9.09 0.87 7.67 10.43 8000 1
## div_est[13] 9.69 0.90 8.25 11.13 4288 1
## div_est[14] 8.10 0.42 7.43 8.76 8000 1
## div_est[15] 10.68 0.55 9.80 11.57 8000 1
## div_est[16] 10.17 0.70 9.04 11.29 8000 1
## div_est[17] 10.50 0.80 9.21 11.74 8000 1
## div_est[18] 11.94 0.64 10.89 12.93 8000 1
## div_est[19] 10.50 0.70 9.34 11.57 8000 1
## div_est[20] 10.18 1.02 8.58 11.82 5259 1
## div_est[21] 8.76 0.59 7.79 9.67 8000 1
## div_est[22] 7.77 0.48 6.99 8.52 8000 1
## div_est[23] 9.15 0.48 8.38 9.91 8000 1
## div_est[24] 7.75 0.54 6.88 8.59 8000 1
## div_est[25] 10.43 0.76 9.26 11.67 8000 1
## div_est[26] 9.53 0.57 8.59 10.44 8000 1
## div_est[27] 9.43 0.94 7.91 10.88 8000 1
## div_est[28] 9.26 0.74 8.04 10.40 8000 1
## div_est[29] 9.18 0.94 7.67 10.65 8000 1
## div_est[30] 6.39 0.44 5.68 7.08 8000 1
## div_est[31] 9.98 0.80 8.67 11.19 8000 1
## div_est[32] 6.70 0.30 6.23 7.19 8000 1
## div_est[33] 9.89 0.43 9.19 10.55 8000 1
## div_est[34] 9.75 0.95 8.29 11.31 8000 1
## div_est[35] 9.43 0.41 8.79 10.11 8000 1
## div_est[36] 11.95 0.77 10.75 13.21 8000 1
## div_est[37] 10.06 0.64 9.10 11.16 8000 1
## div_est[38] 7.80 0.41 7.16 8.46 8000 1
## div_est[39] 8.22 1.01 6.58 9.82 8000 1
## div_est[40] 8.41 0.60 7.48 9.39 8000 1
## div_est[41] 10.00 1.05 8.32 11.65 8000 1
## div_est[42] 10.95 0.64 9.96 11.97 8000 1
## div_est[43] 10.02 0.33 9.48 10.53 8000 1
## div_est[44] 11.06 0.80 9.81 12.36 8000 1
## div_est[45] 8.91 0.98 7.28 10.42 8000 1
## div_est[46] 9.01 0.47 8.26 9.77 7614 1
## div_est[47] 9.95 0.56 9.07 10.85 8000 1
## div_est[48] 10.63 0.88 9.20 11.97 8000 1
## div_est[49] 8.48 0.51 7.68 9.30 8000 1
## div_est[50] 11.49 1.10 9.72 13.19 8000 1
## a 21.42 6.51 10.86 31.51 2382 1
## bA -0.55 0.21 -0.89 -0.22 2509 1
## bR 0.12 0.08 0.01 0.25 2744 1
## sigma 1.12 0.21 0.79 1.44 2385 1
## R code 14.5
dlist <- list(
div_obs=d$Divorce,
div_sd=d$Divorce.SE,
mar_obs=d$Marriage,
mar_sd=d$Marriage.SE,
A=d$MedianAgeMarriage )
m14.2 <- map2stan(
alist(
div_est ~ dnorm(mu,sigma),
mu <- a + bA*A + bR*mar_est[i],
div_obs ~ dnorm(div_est,div_sd),
mar_obs ~ dnorm(mar_est,mar_sd),
a ~ dnorm(0,10),
bA ~ dnorm(0,10),
bR ~ dnorm(0,10),
sigma ~ dcauchy(0,2.5)
) ,
data=dlist ,
start=list(div_est=dlist$div_obs,mar_est=dlist$mar_obs) ,
WAIC=FALSE , iter=5000 , warmup=1000 , chains=3 , cores=3 ,
control=list(adapt_delta=0.95) )
## Precompiling model...
## the number of chains is less than 1; sampling not done
## Dispatching 3 chains to 3 cores.
##
## SAMPLING FOR MODEL 'div_est ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Warning in map2stan(alist(div_est ~ dnorm(mu, sigma), mu <- a + bA * A + : There were 23 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## R code 14.6
library(rethinking)
data(milk)
d <- milk
d$neocortex.prop <- d$neocortex.perc / 100
d$logmass <- log(d$mass)
head(d)
## clade species kcal.per.g perc.fat perc.protein
## 1 Strepsirrhine Eulemur fulvus 0.49 16.60 15.42
## 2 Strepsirrhine E macaco 0.51 19.27 16.91
## 3 Strepsirrhine E mongoz 0.46 14.11 16.85
## 4 Strepsirrhine E rubriventer 0.48 14.91 13.18
## 5 Strepsirrhine Lemur catta 0.60 27.28 19.50
## 6 New World Monkey Alouatta seniculus 0.47 21.22 23.58
## perc.lactose mass neocortex.perc neocortex.prop logmass
## 1 67.98 1.95 55.16 0.5516 0.6678294
## 2 63.82 2.09 NA NA 0.7371641
## 3 69.04 2.51 NA NA 0.9202828
## 4 71.91 1.62 NA NA 0.4824261
## 5 53.22 2.19 NA NA 0.7839015
## 6 55.20 5.25 64.54 0.6454 1.6582281
## R code 14.7
# prep data
data_list <- list(
kcal = d$kcal.per.g,
neocortex = d$neocortex.prop,
logmass = d$logmass )
# fit model
m14.3 <- map2stan(
alist(
kcal ~ dnorm(mu,sigma),
mu <- a + bN*neocortex + bM*logmass,
neocortex ~ dnorm(nu,sigma_N),#!!!!
a ~ dnorm(0,100),
c(bN,bM) ~ dnorm(0,10),
nu ~ dnorm(0.5,1),
sigma_N ~ dcauchy(0,1),
sigma ~ dcauchy(0,1)
) ,
data=data_list , iter=1e4 , chains=2 )
## Imputing 12 missing values (NA) in variable 'neocortex'.
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 4.914 seconds (Warm-up)
## # 5.629 seconds (Sampling)
## # 10.543 seconds (Total)
##
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 6.885 seconds (Warm-up)
## # 6.803 seconds (Sampling)
## # 13.688 seconds (Total)
##
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## Warning: 使われていないコネクション 8 (<-kuma:11429) を閉じます
## Warning: 使われていないコネクション 6 (<-kuma:11429) を閉じます
## Warning: 使われていないコネクション 5 (<-kuma:11429) を閉じます
## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]
## R code 14.8
precis(m14.3,depth=2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## neocortex_impute[1] 0.63 0.05 0.56 0.72 7541 1
## neocortex_impute[2] 0.63 0.05 0.54 0.70 6159 1
## neocortex_impute[3] 0.62 0.05 0.54 0.71 6697 1
## neocortex_impute[4] 0.65 0.05 0.57 0.73 8433 1
## neocortex_impute[5] 0.70 0.05 0.63 0.79 8293 1
## neocortex_impute[6] 0.66 0.05 0.58 0.73 10000 1
## neocortex_impute[7] 0.69 0.05 0.61 0.76 8817 1
## neocortex_impute[8] 0.70 0.05 0.62 0.77 9290 1
## neocortex_impute[9] 0.71 0.05 0.64 0.79 8301 1
## neocortex_impute[10] 0.65 0.05 0.57 0.72 7615 1
## neocortex_impute[11] 0.66 0.05 0.58 0.73 10000 1
## neocortex_impute[12] 0.70 0.05 0.62 0.78 7618 1
## a -0.54 0.48 -1.31 0.22 1749 1
## bN 1.91 0.76 0.74 3.13 1715 1
## bM -0.07 0.02 -0.11 -0.03 2239 1
## nu 0.67 0.01 0.65 0.69 5477 1
## sigma_N 0.06 0.01 0.05 0.08 3748 1
## sigma 0.13 0.02 0.10 0.17 3316 1
## R code 14.9
# prep data
dcc <- d[ complete.cases(d$neocortex.prop) , ]
data_list_cc <- list(
kcal = dcc$kcal.per.g,
neocortex = dcc$neocortex.prop,
logmass = dcc$logmass )
# fit model
m14.3cc <- map2stan(
alist(
kcal ~ dnorm(mu,sigma),
mu <- a + bN*neocortex + bM*logmass,
a ~ dnorm(0,100),
c(bN,bM) ~ dnorm(0,10),
sigma ~ dcauchy(0,1)
) ,
data=data_list_cc , iter=1e4 , chains=2 )
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 1.583 seconds (Warm-up)
## # 1.85 seconds (Sampling)
## # 3.433 seconds (Total)
##
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 1.64 seconds (Warm-up)
## # 3.147 seconds (Sampling)
## # 4.787 seconds (Total)
##
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]
precis(m14.3cc)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -1.07 0.58 -2.01 -0.17 2218 1
## bN 2.76 0.90 1.37 4.22 2189 1
## bM -0.10 0.03 -0.14 -0.05 2535 1
## sigma 0.14 0.03 0.09 0.18 2553 1
## R code 14.10
m14.4 <- map2stan(
alist(
kcal ~ dnorm(mu,sigma),
mu <- a + bN*neocortex + bM*logmass,
neocortex ~ dnorm(nu,sigma_N),
nu <- a_N + gM*logmass,
a ~ dnorm(0,100),
c(bN,bM,gM) ~ dnorm(0,10),
a_N ~ dnorm(0.5,1),
sigma_N ~ dcauchy(0,1),
sigma ~ dcauchy(0,1)
) ,
data=data_list , iter=1e4 , chains=2 )
## Imputing 12 missing values (NA) in variable 'neocortex'.
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 11.188 seconds (Warm-up)
## # 11.416 seconds (Sampling)
## # 22.604 seconds (Total)
##
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 11.396 seconds (Warm-up)
## # 11.41 seconds (Sampling)
## # 22.806 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 2"
## count
## Exception thrown at line 37: normal_log: Scale parameter is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'kcal ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]
precis(m14.4,depth=2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## neocortex_impute[1] 0.63 0.04 0.58 0.69 10000 1
## neocortex_impute[2] 0.63 0.04 0.57 0.69 10000 1
## neocortex_impute[3] 0.62 0.04 0.56 0.68 10000 1
## neocortex_impute[4] 0.65 0.03 0.59 0.70 10000 1
## neocortex_impute[5] 0.66 0.04 0.61 0.72 10000 1
## neocortex_impute[6] 0.63 0.04 0.57 0.69 10000 1
## neocortex_impute[7] 0.68 0.03 0.62 0.73 10000 1
## neocortex_impute[8] 0.70 0.03 0.64 0.75 10000 1
## neocortex_impute[9] 0.71 0.03 0.66 0.77 10000 1
## neocortex_impute[10] 0.66 0.04 0.61 0.72 10000 1
## neocortex_impute[11] 0.68 0.03 0.62 0.73 10000 1
## neocortex_impute[12] 0.74 0.04 0.68 0.80 10000 1
## a -0.87 0.49 -1.65 -0.11 2876 1
## bN 2.44 0.76 1.24 3.64 2853 1
## bM -0.09 0.02 -0.13 -0.05 3731 1
## gM 0.02 0.01 0.01 0.03 7127 1
## a_N 0.64 0.01 0.62 0.66 6145 1
## sigma_N 0.04 0.01 0.03 0.05 4679 1
## sigma 0.13 0.02 0.09 0.16 5076 1
## R code 14.11
nc_missing <- ifelse( is.na(d$neocortex.prop) , 1 , 0 )
nc_missing <- sapply( 1:length(nc_missing) ,
function(n) nc_missing[n]*sum(nc_missing[1:n]) )
nc_missing
## [1] 0 1 2 3 4 0 0 0 5 0 0 0 0 6 7 0 8 0 9 0 10 0 11
## [24] 0 0 12 0 0 0
## R code 14.12
nc <- ifelse( is.na(d$neocortex.prop) , -1 , d$neocortex.prop )
## R code 14.13
model_code <- '
data{
int N;
int nc_num_missing;
vector[N] kcal;
real neocortex[N];
vector[N] logmass;
int nc_missing[N];
}
parameters{
real alpha;
real<lower=0> sigma;
real bN;
real bM;
vector[nc_num_missing] nc_impute;
real mu_nc;
real<lower=0> sigma_nc;
}
model{
vector[N] mu;
vector[N] nc_merged;
alpha ~ normal(0,10);
bN ~ normal(0,10);
bM ~ normal(0,10);
mu_nc ~ normal(0.5,1);
sigma ~ cauchy(0,1);
sigma_nc ~ cauchy(0,1);
// merge missing and observed
for ( i in 1:N ) {
nc_merged[i] <- neocortex[i];
if ( nc_missing[i] > 0 ) nc_merged[i] <- nc_impute[nc_missing[i]];
}
// imputation
nc_merged ~ normal( mu_nc , sigma_nc );
// regression
mu <- alpha + bN*nc_merged + bM*logmass;
kcal ~ normal( mu , sigma );
}'
## R code 14.14
data_list <- list(
N = nrow(d),
kcal = d$kcal.per.g,
neocortex = nc,
logmass = d$logmass,
nc_missing = nc_missing,
nc_num_missing = max(nc_missing)
)
start <- list(
alpha=mean(d$kcal.per.g), sigma=sd(d$kcal.per.g),
bN=0, bM=0, mu_nc=0.68, sigma_nc=0.06,
nc_impute=rep( 0.5 , max(nc_missing) )
)
library(rstan)
m14.3stan <- stan( model_code=model_code , data=data_list , init=list(start) ,
iter=1e4 , chains=1 )
##
## SAMPLING FOR MODEL 'b04c58eddd0abbeb3a901686231b923e' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 4.865 seconds (Warm-up)
## # 6.94 seconds (Sampling)
## # 11.805 seconds (Total)
## R code 14.15
set.seed(100)
x <- c( rnorm(10) , NA )
y <- c( rnorm(10,x) , 100 )
d <- list(x=x,y=y)