#Metropolis algorithm
## R code 8.1
num_weeks <- 1e5
positions <- rep(0,num_weeks)
current <- 10
for ( i in 1:num_weeks ) {
# record current position
positions[i] <- current
# flip coin to generate proposal
proposal <- current + sample( c(-1,1) , size=1 )
# now make sure he loops around the archipelago
if ( proposal < 1 ) proposal <- 10
if ( proposal > 10 ) proposal <- 1
# move?
prob_move <- proposal/current
current <- ifelse( runif(1) < prob_move , proposal , current )
}
## R code 8.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(rugged)
d <- rugged
d$log_gdp <- log(d$rgdppc_2000)
dd <- d[ complete.cases(d$rgdppc_2000) , ]
head(dd)
## isocode isonum country rugged rugged_popw rugged_slope
## 3 AGO 24 Angola 0.858 0.714 2.274
## 5 ALB 8 Albania 3.427 1.597 10.451
## 8 ARE 784 United Arab Emirates 0.769 0.316 2.112
## 9 ARG 32 Argentina 0.775 0.220 2.268
## 10 ARM 51 Armenia 2.688 0.934 8.178
## 12 ATG 28 Antigua and Barbuda 0.006 0.003 0.012
## rugged_lsd rugged_pc land_area lat lon soil desert tropical
## 3 0.228 4.906 124670 -12.299 17.551 26.676 0.425 44.346
## 5 1.006 62.133 2740 41.143 20.070 68.088 0.000 0.000
## 8 0.191 6.142 8360 23.913 54.331 0.000 77.280 0.000
## 9 0.226 9.407 273669 -35.396 -65.170 35.678 0.000 0.000
## 10 0.799 50.556 2820 40.294 44.938 30.148 0.000 0.000
## 12 0.003 0.000 44 17.271 -61.800 100.000 0.000 100.000
## dist_coast near_coast gemstones rgdppc_2000 rgdppc_1950_m rgdppc_1975_m
## 3 0.428 13.1587 47756 1794.729 1051.822 1073.036
## 5 0.048 94.6919 0 3703.113 1001.339 2289.472
## 8 0.065 75.7464 0 20604.460 15797.558 25465.002
## 9 0.352 13.0167 0 12173.680 4986.725 8122.497
## 10 0.348 0.0000 0 2421.985 NA NA
## 12 0.001 100.0000 0 10022.030 NA NA
## rgdppc_2000_m rgdppc_1950_2000_m q_rule_law cont_africa cont_asia
## 3 765.215 1106.763 -1.567 1 0
## 5 2741.420 1931.784 -0.820 0 0
## 8 17567.883 20119.992 0.913 0 1
## 9 8543.558 6926.810 0.033 0 0
## 10 4565.035 NA -0.453 0 1
## 12 NA NA 0.990 0 0
## cont_europe cont_oceania cont_north_america cont_south_america
## 3 0 0 0 0
## 5 1 0 0 0
## 8 0 0 0 0
## 9 0 0 0 1
## 10 0 0 0 0
## 12 0 0 1 0
## legor_gbr legor_fra legor_soc legor_deu legor_sca colony_esp colony_gbr
## 3 0 1 0 0 0 0 0
## 5 0 0 1 0 0 0 0
## 8 1 0 0 0 0 0 1
## 9 0 1 0 0 0 1 0
## 10 0 0 1 0 0 0 0
## 12 1 0 0 0 0 0 1
## colony_fra colony_prt colony_oeu africa_region_n africa_region_s
## 3 0 1 0 0 0
## 5 0 0 0 0 0
## 8 0 0 0 0 0
## 9 0 0 0 0 0
## 10 0 0 0 0 0
## 12 0 0 0 0 0
## africa_region_w africa_region_e africa_region_c slave_exports
## 3 0 0 1 3610000
## 5 0 0 0 0
## 8 0 0 0 0
## 9 0 0 0 0
## 10 0 0 0 0
## 12 0 0 0 0
## dist_slavemkt_atlantic dist_slavemkt_indian dist_slavemkt_saharan
## 3 5.669 6.981 4.926
## 5 NA NA NA
## 8 NA NA NA
## 9 NA NA NA
## 10 NA NA NA
## 12 NA NA NA
## dist_slavemkt_redsea pop_1400 european_descent log_gdp
## 3 3.872 1223208 2.000 7.492609
## 5 NA 200000 100.000 8.216929
## 8 NA 19200 0.000 9.933263
## 9 NA 276632 89.889 9.407032
## 10 NA 105743 0.500 7.792343
## 12 NA 747 NA 9.212541
## R code 8.3
m8.1 <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged + bA*cont_africa +
bAR*rugged*cont_africa ,
a ~ dnorm(0,100),
bR ~ dnorm(0,10),
bA ~ dnorm(0,10),
bAR ~ dnorm(0,10),
sigma ~ dunif(0,10)
) ,
data=dd )
precis(m8.1)
## Mean StdDev 5.5% 94.5%
## a 9.22 0.14 9.00 9.44
## bR -0.20 0.08 -0.32 -0.08
## bA -1.95 0.22 -2.31 -1.59
## bAR 0.39 0.13 0.19 0.60
## sigma 0.93 0.05 0.85 1.01
## R code 8.4
dd.trim <- dd[ , c("log_gdp","rugged","cont_africa") ]
str(dd.trim)
## 'data.frame': 170 obs. of 3 variables:
## $ log_gdp : num 7.49 8.22 9.93 9.41 7.79 ...
## $ rugged : num 0.858 3.427 0.769 0.775 2.688 ...
## $ cont_africa: int 1 0 0 0 0 0 0 0 0 1 ...
## R code 8.5
m8.1stan <- map2stan(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged + bA*cont_africa + bAR*rugged*cont_africa ,
a ~ dnorm(0,100),
bR ~ dnorm(0,10),
bA ~ dnorm(0,10),
bAR ~ dnorm(0,10),
sigma ~ dcauchy(0,2)
) ,
data=dd.trim )
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 0.578 seconds (Warm-up)
## # 0.562 seconds (Sampling)
## # 1.14 seconds (Total)
##
##
## SAMPLING FOR MODEL 'log_gdp ~ 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
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## R code 8.6
precis(m8.1stan)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.23 0.13 9.02 9.43 223 1.00
## bR -0.21 0.07 -0.33 -0.08 293 1.00
## bA -1.95 0.21 -2.28 -1.61 231 1.01
## bAR 0.39 0.13 0.20 0.60 190 1.00
## sigma 0.95 0.05 0.87 1.04 511 1.00
## R code 8.7
m8.1stan_4chains <- map2stan( m8.1stan , chains=4 , cores=4 )
## Dispatching 4 chains to 4 cores.
##
## SAMPLING FOR MODEL 'log_gdp ~ 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
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
precis(m8.1stan_4chains)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.23 0.14 9.00 9.44 1041 1
## bR -0.20 0.08 -0.33 -0.08 1015 1
## bA -1.95 0.23 -2.34 -1.62 718 1
## bAR 0.40 0.13 0.18 0.60 823 1
## sigma 0.95 0.05 0.86 1.03 1649 1
## R code 8.8
post <- extract.samples( m8.1stan )
str(post)
## List of 5
## $ a : num [1:1000(1d)] 9.2 9.4 9.48 9.41 9.29 ...
## $ bR : num [1:1000(1d)] -0.172 -0.308 -0.289 -0.233 -0.24 ...
## $ bA : num [1:1000(1d)] -1.77 -2.18 -2.23 -1.95 -2.31 ...
## $ bAR : num [1:1000(1d)] 0.418 0.499 0.541 0.37 0.547 ...
## $ sigma: num [1:1000(1d)] 0.97 0.936 1.074 0.981 0.937 ...
## R code 8.9
pairs(post)
## R code 8.10
pairs(m8.1stan)
## R code 8.11
show(m8.1stan)
## map2stan model fit
## 1000 samples from 1 chain
##
## Formula:
## log_gdp ~ dnorm(mu, sigma)
## mu <- a + bR * rugged + bA * cont_africa + bAR * rugged * cont_africa
## a ~ dnorm(0, 100)
## bR ~ dnorm(0, 10)
## bA ~ dnorm(0, 10)
## bAR ~ dnorm(0, 10)
## sigma ~ dcauchy(0, 2)
##
## Log-likelihood at expected values: -229.45
## Deviance: 458.89
## DIC: 468.33
## Effective number of parameters (pD): 4.72
##
## WAIC (SE): 468.75 (14.7)
## pWAIC: 4.82
## R code 8.12
plot(m8.1stan)
## R code 8.13
y <- c(-1,1)
m8.2 <- map2stan(
alist(
y ~ dnorm( mu , sigma ) ,
mu <- alpha
) ,
data=list(y=y) , start=list(alpha=0,sigma=1) ,
chains=2 , iter=4000 , warmup=1000 )
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 1, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 1, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 1, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 1, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 1, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 1, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 1, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 0.218 seconds (Warm-up)
## # 0.343 seconds (Sampling)
## # 0.561 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 2, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 2, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 2, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 2, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 2, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 2, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 2, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 0.047 seconds (Warm-up)
## # 0.502 seconds (Sampling)
## # 0.549 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ 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
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
## Warning in map2stan(alist(y ~ dnorm(mu, sigma), mu <- alpha), data = list(y = y), : There were 769 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## R code 8.14
precis(m8.2)
## Warning in precis(m8.2): There were 769 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 85946438 325586539 -60711981.15 157804787 43 1.1
## sigma 981741780 27253368293 926.39 449191597 2106 1.0
## R code 8.15
m8.3 <- map2stan(
alist(
y ~ dnorm( mu , sigma ) ,
mu <- alpha ,
alpha ~ dnorm( 1 , 10 ) ,
sigma ~ dcauchy( 0 , 1 )
) ,
data=list(y=y) , start=list(alpha=0,sigma=1) ,
chains=2 , iter=4000 , warmup=1000 )
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 1, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 1, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 1, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 1, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 1, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 1, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 1, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 0.016 seconds (Warm-up)
## # 0.078 seconds (Sampling)
## # 0.094 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 2, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 2, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 2, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 2, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 2, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 2, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 2, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 0.016 seconds (Warm-up)
## # 0.062 seconds (Sampling)
## # 0.078 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ 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
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
## Warning: 使われていないコネクション 10 (<-kuma:11432) を閉じます
## Warning: 使われていないコネクション 9 (<-kuma:11432) を閉じます
## Warning: 使われていないコネクション 8 (<-kuma:11432) を閉じます
## Warning: 使われていないコネクション 7 (<-kuma:11432) を閉じます
## [ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
precis(m8.3)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 0.06 2.08 -2.66 2.34 502 1
## sigma 2.16 2.66 0.47 3.67 594 1
## R code 8.16
y <- rcauchy(1e4,0,5)
mu <- sapply( 1:length(y) , function(i) sum(y[1:i])/i )
plot(mu,type="l")
## R code 8.17
y <- rnorm( 100 , mean=0 , sd=1 )
## R code 8.18
m8.4 <- map2stan(
alist(
y ~ dnorm( mu , sigma ) ,
mu <- a1 + a2 ,
sigma ~ dcauchy( 0 , 1 )
) ,
data=list(y=y) , start=list(a1=0,a2=0,sigma=1) ,
chains=2 , iter=4000 , warmup=1000 )
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 1, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 1, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 1, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 1, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 1, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 1, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 1, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 9.814 seconds (Warm-up)
## # 35.061 seconds (Sampling)
## # 44.875 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 1"
## count
## Exception thrown at line 16: 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 'y ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 2, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 2, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 2, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 2, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 2, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 2, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 2, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 9.31 seconds (Warm-up)
## # 35.211 seconds (Sampling)
## # 44.521 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ 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
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
precis(m8.4)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a1 -834.85 1660.14 -3138.49 1217.30 1 3.84
## a2 835.00 1660.14 -1217.01 3138.68 1 3.84
## sigma 0.88 0.06 0.79 0.97 40 1.05
## R code 8.19
m8.5 <- map2stan(
alist(
y ~ dnorm( mu , sigma ) ,
mu <- a1 + a2 ,
a1 ~ dnorm( 0 , 10 ) ,
a2 ~ dnorm( 0 , 10 ) ,
sigma ~ dcauchy( 0 , 1 )
) ,
data=list(y=y) , start=list(a1=0,a2=0,sigma=1) ,
chains=2 , iter=4000 , warmup=1000 )
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 1, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 1, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 1, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 1, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 1, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 1, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 1, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 1.563 seconds (Warm-up)
## # 4.227 seconds (Sampling)
## # 5.79 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 1"
## count
## Exception thrown at line 18: 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 'y ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 2, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 2, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 2, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 2, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 2, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 2, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 2, Iteration: 4000 / 4000 [100%] (Sampling)
## # Elapsed Time: 1.628 seconds (Warm-up)
## # 5.309 seconds (Sampling)
## # 6.937 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ 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
## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]
precis(m8.5)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a1 0.09 7.02 -11.07 11.15 956 1
## a2 0.07 7.02 -11.12 11.12 956 1
## sigma 0.89 0.06 0.79 0.99 1267 1
## R code 8.20
mp <- map2stan(
alist(
a ~ dnorm(0,1),
b ~ dcauchy(0,1)
),
data=list(y=1),
start=list(a=0,b=0),
iter=1e4, warmup=100 , WAIC=FALSE )
##
## SAMPLING FOR MODEL 'a ~ dnorm(0, 1)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 101 / 10000 [ 1%] (Sampling)
## Chain 1, Iteration: 1100 / 10000 [ 11%] (Sampling)
## Chain 1, Iteration: 2100 / 10000 [ 21%] (Sampling)
## Chain 1, Iteration: 3100 / 10000 [ 31%] (Sampling)
## Chain 1, Iteration: 4100 / 10000 [ 41%] (Sampling)
## Chain 1, Iteration: 5100 / 10000 [ 51%] (Sampling)
## Chain 1, Iteration: 6100 / 10000 [ 61%] (Sampling)
## Chain 1, Iteration: 7100 / 10000 [ 71%] (Sampling)
## Chain 1, Iteration: 8100 / 10000 [ 81%] (Sampling)
## Chain 1, Iteration: 9100 / 10000 [ 91%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0.208 seconds (Sampling)
## # 0.208 seconds (Total)
##
##
## SAMPLING FOR MODEL 'a ~ dnorm(0, 1)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## R code 8.21
N <- 100 # number of individuals
height <- rnorm(N,10,2) # sim total height of each
leg_prop <- runif(N,0.4,0.5) # leg as proportion of height
leg_left <- leg_prop*height + # sim left leg as proportion + error
rnorm( N , 0 , 0.02 )
leg_right <- leg_prop*height + # sim right leg as proportion + error
rnorm( N , 0 , 0.02 )
# combine into data frame
d <- data.frame(height,leg_left,leg_right)
head(d)
## height leg_left leg_right
## 1 10.502628 4.374049 4.345735
## 2 7.742391 3.554417 3.569487
## 3 11.111753 5.361891 5.355601
## 4 9.719459 4.740331 4.767011
## 5 9.936417 4.803539 4.789984
## 6 9.963384 4.771129 4.806627
## R code 8.22
m5.8s <- map2stan(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + bl*leg_left + br*leg_right ,
a ~ dnorm( 10 , 100 ) ,
bl ~ dnorm( 2 , 10 ) ,
br ~ dnorm( 2 , 10 ) ,
sigma ~ dcauchy( 0 , 1 )
) ,
data=d, chains=4,
start=list(a=10,bl=0,br=0,sigma=1) )
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 10.599 seconds (Warm-up)
## # 14.015 seconds (Sampling)
## # 24.614 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 11.021 seconds (Warm-up)
## # 14.285 seconds (Sampling)
## # 25.306 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 11.689 seconds (Warm-up)
## # 13.47 seconds (Sampling)
## # 25.159 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 3"
## count
## Exception thrown at line 22: 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 'height ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 11.224 seconds (Warm-up)
## # 12.131 seconds (Sampling)
## # 23.355 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ 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
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 8.23
m5.8s2 <- map2stan(
alist(
height ~ dnorm( mu , sigma ) ,
mu <- a + bl*leg_left + br*leg_right ,
a ~ dnorm( 10 , 100 ) ,
bl ~ dnorm( 2 , 10 ) ,
br ~ dnorm( 2 , 10 ) & T[0,] ,
sigma ~ dcauchy( 0 , 1 )
) ,
data=d, chains=4,
start=list(a=10,bl=0,br=0,sigma=1) )
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 5.388 seconds (Warm-up)
## # 5.884 seconds (Sampling)
## # 11.272 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 5.116 seconds (Warm-up)
## # 5.873 seconds (Sampling)
## # 10.989 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 5.392 seconds (Warm-up)
## # 5.966 seconds (Sampling)
## # 11.358 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
## # Elapsed Time: 5.917 seconds (Warm-up)
## # 5.178 seconds (Sampling)
## # 11.095 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 4"
## count
## Exception thrown at line 22: 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 'height ~ 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
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Warning in map2stan(alist(height ~ dnorm(mu, sigma), mu <- a + bl * leg_left + : There were 1423 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.