## R code 13.1
a <- 3.5 # average morning wait time
b <- (-1) # average difference afternoon wait time
sigma_a <- 1 # std dev in intercepts
sigma_b <- 0.5 # std dev in slopes
rho <- (-0.7) # correlation between intercepts and slopes
## R code 13.2
Mu <- c( a , b )
## R code 13.3
cov_ab <- sigma_a*sigma_b*rho
Sigma <- matrix( c(sigma_a^2,cov_ab,cov_ab,sigma_b^2) , ncol=2 )
## R code 13.4
matrix( c(1,2,3,4) , nrow=2 , ncol=2 )
## [,1] [,2]
## [1,] 1 3
## [2,] 2 4
## R code 13.5
sigmas <- c(sigma_a,sigma_b) # standard deviations
Rho <- matrix( c(1,rho,rho,1) , nrow=2 ) # correlation matrix
# now matrix multiply to get covariance matrix
Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas)
## R code 13.6
N_cafes <- 20
## R code 13.7
library(MASS)
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)
library(ellipse)
set.seed(5) # used to replicate example
vary_effects <- mvrnorm( N_cafes , Mu , Sigma )
## R code 13.8
a_cafe <- vary_effects[,1]
b_cafe <- vary_effects[,2]
## R code 13.9
par(mfrow=c(1,1))
plot( a_cafe , b_cafe , col=rangi2 ,
xlab="intercepts (a_cafe)" , ylab="slopes (b_cafe)" )
# overlay population distribution
#library(ellipse)
for ( l in c(0.1,0.3,0.5,0.8,0.99) )
lines(ellipse(Sigma,centre=Mu,level=l),col=col.alpha("black",0.2))
## R code 13.10
N_visits <- 10
afternoon <- rep(0:1,N_visits*N_cafes/2)
cafe_id <- rep( 1:N_cafes , each=N_visits )
mu <- a_cafe[cafe_id] + b_cafe[cafe_id]*afternoon
sigma <- 0.5 # std dev within cafes
wait <- rnorm( N_visits*N_cafes , mu , sigma )
d <- data.frame( cafe=cafe_id , afternoon=afternoon , wait=wait )
## R code 13.11
R <- rlkjcorr( 1e4 , K=2 , eta=2 )
dens( R[,1,2] , xlab="correlation" )
## R code 13.12
m13.1 <- map2stan(
alist(
wait ~ dnorm( mu , sigma ),
mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
c(a_cafe,b_cafe)[cafe] ~ dmvnorm2(c(a,b),sigma_cafe,Rho),
a ~ dnorm(0,10),
b ~ dnorm(0,10),
sigma_cafe ~ dcauchy(0,2),
sigma ~ dcauchy(0,2),
Rho ~ dlkjcorr(2)
) ,
data=d ,
iter=5000 , warmup=2000 , chains=2 )
##
## SAMPLING FOR MODEL 'wait ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 1, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 1, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1, Iteration: 5000 / 5000 [100%] (Sampling)
## # Elapsed Time: 2.972 seconds (Warm-up)
## # 4.93 seconds (Sampling)
## # 7.902 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 1"
## count
## -0 0 1
## 0 -0 1
## validate transformed params: SRS_sigma_cafeRho is not positive definite: 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 'wait ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2, Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2, Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2, Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2, Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2, Iteration: 2001 / 5000 [ 40%] (Sampling)
## Chain 2, Iteration: 2500 / 5000 [ 50%] (Sampling)
## Chain 2, Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2, Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2, Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2, Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2, Iteration: 5000 / 5000 [100%] (Sampling)
## # Elapsed Time: 3.295 seconds (Warm-up)
## # 4.271 seconds (Sampling)
## # 7.566 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 2"
## count
## validate transformed params: SRS_sigma_cafeRho is not symmetric. SRS_sigma_cafeRho[1,2] = 1.#INF, bu 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 'wait ~ 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(wait ~ dnorm(mu, sigma), mu <- a_cafe[cafe] + : There were 109 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## R code 13.13
post <- extract.samples(m13.1)
dens( post$Rho[,1,2] )
## R code 13.14
# compute unpooled estimates directly from data
a1 <- sapply( 1:N_cafes ,
function(i) mean(wait[cafe_id==i & afternoon==0]) )
b1 <- sapply( 1:N_cafes ,
function(i) mean(wait[cafe_id==i & afternoon==1]) ) - a1
# extract posterior means of partially pooled estimates
post <- extract.samples(m13.1)
a2 <- apply( post$a_cafe , 2 , mean )
b2 <- apply( post$b_cafe , 2 , mean )
# plot both and connect with lines
plot( a1 , b1 , xlab="intercept" , ylab="slope" ,
pch=16 , col=rangi2 , ylim=c( min(b1)-0.1 , max(b1)+0.1 ) ,
xlim=c( min(a1)-0.1 , max(a1)+0.1 ) )
points( a2 , b2 , pch=1 )
for ( i in 1:N_cafes ) lines( c(a1[i],a2[i]) , c(b1[i],b2[i]) )
## R code 13.15
# compute posterior mean bivariate Gaussian
Mu_est <- c( mean(post$a) , mean(post$b) )
rho_est <- mean( post$Rho[,1,2] )
sa_est <- mean( post$sigma_cafe[,1] )
sb_est <- mean( post$sigma_cafe[,2] )
cov_ab <- sa_est*sb_est*rho_est
Sigma_est <- matrix( c(sa_est^2,cov_ab,cov_ab,sb_est^2) , ncol=2 )
# draw contours
#library(ellipse)
for ( l in c(0.1,0.3,0.5,0.8,0.99) )
lines(ellipse(Sigma_est,centre=Mu_est,level=l),
col=col.alpha("black",0.2))
## R code 13.16
# convert varying effects to waiting times
wait_morning_1 <- (a1)
wait_afternoon_1 <- (a1 + b1)
wait_morning_2 <- (a2)
wait_afternoon_2 <- (a2 + b2)
## R code 13.17
library(rethinking)
data(UCBadmit)
d <- UCBadmit
d$male <- ifelse( d$applicant.gender=="male" , 1 , 0 )
d$dept_id <- coerce_index( d$dept )
## R code 13.18
m13.2 <- map2stan(
alist(
admit ~ dbinom( applications , p ),
logit(p) <- a_dept[dept_id] + bm*male,
a_dept[dept_id] ~ dnorm( a , sigma_dept ),
a ~ dnorm(0,10),
bm ~ dnorm(0,1),
sigma_dept ~ dcauchy(0,2)
) ,
data=d , warmup=500 , iter=4500 , chains=3 )
## Warning: Variable 'applicant.gender' contains dots '.'.
## Will attempt to remove dots internally.
## Warning in FUN(X[[i]], ...): data with name dept is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name applicant_gender is not numeric
## and not used
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4500 [ 0%] (Warmup)
## Chain 1, Iteration: 450 / 4500 [ 10%] (Warmup)
## Chain 1, Iteration: 501 / 4500 [ 11%] (Sampling)
## Chain 1, Iteration: 950 / 4500 [ 21%] (Sampling)
## Chain 1, Iteration: 1400 / 4500 [ 31%] (Sampling)
## Chain 1, Iteration: 1850 / 4500 [ 41%] (Sampling)
## Chain 1, Iteration: 2300 / 4500 [ 51%] (Sampling)
## Chain 1, Iteration: 2750 / 4500 [ 61%] (Sampling)
## Chain 1, Iteration: 3200 / 4500 [ 71%] (Sampling)
## Chain 1, Iteration: 3650 / 4500 [ 81%] (Sampling)
## Chain 1, Iteration: 4100 / 4500 [ 91%] (Sampling)
## Chain 1, Iteration: 4500 / 4500 [100%] (Sampling)
## # Elapsed Time: 0.052 seconds (Warm-up)
## # 0.343 seconds (Sampling)
## # 0.395 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4500 [ 0%] (Warmup)
## Chain 2, Iteration: 450 / 4500 [ 10%] (Warmup)
## Chain 2, Iteration: 501 / 4500 [ 11%] (Sampling)
## Chain 2, Iteration: 950 / 4500 [ 21%] (Sampling)
## Chain 2, Iteration: 1400 / 4500 [ 31%] (Sampling)
## Chain 2, Iteration: 1850 / 4500 [ 41%] (Sampling)
## Chain 2, Iteration: 2300 / 4500 [ 51%] (Sampling)
## Chain 2, Iteration: 2750 / 4500 [ 61%] (Sampling)
## Chain 2, Iteration: 3200 / 4500 [ 71%] (Sampling)
## Chain 2, Iteration: 3650 / 4500 [ 81%] (Sampling)
## Chain 2, Iteration: 4100 / 4500 [ 91%] (Sampling)
## Chain 2, Iteration: 4500 / 4500 [100%] (Sampling)
## # Elapsed Time: 0.06 seconds (Warm-up)
## # 0.41 seconds (Sampling)
## # 0.47 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 4500 [ 0%] (Warmup)
## Chain 3, Iteration: 450 / 4500 [ 10%] (Warmup)
## Chain 3, Iteration: 501 / 4500 [ 11%] (Sampling)
## Chain 3, Iteration: 950 / 4500 [ 21%] (Sampling)
## Chain 3, Iteration: 1400 / 4500 [ 31%] (Sampling)
## Chain 3, Iteration: 1850 / 4500 [ 41%] (Sampling)
## Chain 3, Iteration: 2300 / 4500 [ 51%] (Sampling)
## Chain 3, Iteration: 2750 / 4500 [ 61%] (Sampling)
## Chain 3, Iteration: 3200 / 4500 [ 71%] (Sampling)
## Chain 3, Iteration: 3650 / 4500 [ 81%] (Sampling)
## Chain 3, Iteration: 4100 / 4500 [ 91%] (Sampling)
## Chain 3, Iteration: 4500 / 4500 [100%] (Sampling)
## # Elapsed Time: 0.06 seconds (Warm-up)
## # 0.32 seconds (Sampling)
## # 0.38 seconds (Total)
## Warning in FUN(X[[i]], ...): data with name dept is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name applicant_gender is not numeric
## and not used
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' 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
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
precis( m13.2 , depth=2 ) # depth=2 to display vector parameters
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_dept[1] 0.67 0.10 0.51 0.83 4426 1
## a_dept[2] 0.63 0.12 0.45 0.82 4849 1
## a_dept[3] -0.58 0.08 -0.69 -0.45 7232 1
## a_dept[4] -0.62 0.09 -0.76 -0.48 5772 1
## a_dept[5] -1.06 0.10 -1.21 -0.90 8026 1
## a_dept[6] -2.61 0.16 -2.85 -2.35 8763 1
## a -0.59 0.67 -1.54 0.47 4331 1
## bm -0.10 0.08 -0.23 0.03 3576 1
## sigma_dept 1.49 0.59 0.71 2.19 4419 1
## R code 13.19
m13.3 <- map2stan(
alist(
admit ~ dbinom( applications , p ),
logit(p) <- a_dept[dept_id] +
bm_dept[dept_id]*male,
c(a_dept,bm_dept)[dept_id] ~ dmvnorm2( c(a,bm) , sigma_dept , Rho ),
a ~ dnorm(0,10),
bm ~ dnorm(0,1),
sigma_dept ~ dcauchy(0,2),
Rho ~ dlkjcorr(2)
) ,
data=d , warmup=1000 , iter=5000 , chains=4 , cores=3 )
## Warning: Variable 'applicant.gender' contains dots '.'.
## Will attempt to remove dots internally.
## Precompiling model...
## Warning in FUN(X[[i]], ...): data with name dept is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name applicant_gender is not numeric
## and not used
## the number of chains is less than 1; sampling not done
## Dispatching 4 chains to 3 cores.
## Warning in FUN(X[[i]], ...): data with name dept is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name applicant_gender is not numeric
## and not used
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' 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
## [ 1600 / 16000 ]
[ 3200 / 16000 ]
[ 4800 / 16000 ]
[ 6400 / 16000 ]
[ 8000 / 16000 ]
[ 9600 / 16000 ]
[ 11200 / 16000 ]
[ 12800 / 16000 ]
[ 14400 / 16000 ]
[ 16000 / 16000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
## Warning in map2stan(alist(admit ~ dbinom(applications, p), logit(p) <- a_dept[dept_id] + : There were 6 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## R code 13.20
precis(m13.3,pars=c("a_dept","bm_dept"),depth=2)
## Warning in precis(m13.3, pars = c("a_dept", "bm_dept"), depth = 2): There were 6 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
## bm_dept[1] -0.79 0.27 -1.22 -0.37 5577 1
## bm_dept[2] -0.21 0.33 -0.73 0.31 6454 1
## bm_dept[3] 0.08 0.14 -0.14 0.31 9509 1
## bm_dept[4] -0.09 0.14 -0.32 0.12 8414 1
## bm_dept[5] 0.12 0.19 -0.19 0.41 9592 1
## bm_dept[6] -0.13 0.27 -0.54 0.33 7422 1
## a_dept[1] 1.30 0.25 0.91 1.72 5819 1
## a_dept[2] 0.74 0.33 0.20 1.24 6511 1
## a_dept[3] -0.65 0.09 -0.79 -0.52 9929 1
## a_dept[4] -0.62 0.10 -0.79 -0.45 8498 1
## a_dept[5] -1.13 0.11 -1.31 -0.95 9899 1
## a_dept[6] -2.60 0.20 -2.92 -2.27 8653 1
#plot( precis(m13.3,pars=c("a_dept","bm_dept"),depth=2) )
## R code 13.21
m13.4 <- map2stan(
alist(
admit ~ dbinom( applications , p ),
logit(p) <- a_dept[dept_id],
a_dept[dept_id] ~ dnorm( a , sigma_dept ),
a ~ dnorm(0,10),
sigma_dept ~ dcauchy(0,2)
) ,
data=d , warmup=500 , iter=4500 , chains=3 )
## Warning: Variable 'applicant.gender' contains dots '.'.
## Will attempt to remove dots internally.
## Warning in FUN(X[[i]], ...): data with name dept is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name applicant_gender is not numeric
## and not used
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4500 [ 0%] (Warmup)
## Chain 1, Iteration: 450 / 4500 [ 10%] (Warmup)
## Chain 1, Iteration: 501 / 4500 [ 11%] (Sampling)
## Chain 1, Iteration: 950 / 4500 [ 21%] (Sampling)
## Chain 1, Iteration: 1400 / 4500 [ 31%] (Sampling)
## Chain 1, Iteration: 1850 / 4500 [ 41%] (Sampling)
## Chain 1, Iteration: 2300 / 4500 [ 51%] (Sampling)
## Chain 1, Iteration: 2750 / 4500 [ 61%] (Sampling)
## Chain 1, Iteration: 3200 / 4500 [ 71%] (Sampling)
## Chain 1, Iteration: 3650 / 4500 [ 81%] (Sampling)
## Chain 1, Iteration: 4100 / 4500 [ 91%] (Sampling)
## Chain 1, Iteration: 4500 / 4500 [100%] (Sampling)
## # Elapsed Time: 0.04 seconds (Warm-up)
## # 0.27 seconds (Sampling)
## # 0.31 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4500 [ 0%] (Warmup)
## Chain 2, Iteration: 450 / 4500 [ 10%] (Warmup)
## Chain 2, Iteration: 501 / 4500 [ 11%] (Sampling)
## Chain 2, Iteration: 950 / 4500 [ 21%] (Sampling)
## Chain 2, Iteration: 1400 / 4500 [ 31%] (Sampling)
## Chain 2, Iteration: 1850 / 4500 [ 41%] (Sampling)
## Chain 2, Iteration: 2300 / 4500 [ 51%] (Sampling)
## Chain 2, Iteration: 2750 / 4500 [ 61%] (Sampling)
## Chain 2, Iteration: 3200 / 4500 [ 71%] (Sampling)
## Chain 2, Iteration: 3650 / 4500 [ 81%] (Sampling)
## Chain 2, Iteration: 4100 / 4500 [ 91%] (Sampling)
## Chain 2, Iteration: 4500 / 4500 [100%] (Sampling)
## # Elapsed Time: 0.03 seconds (Warm-up)
## # 0.261 seconds (Sampling)
## # 0.291 seconds (Total)
##
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 4500 [ 0%] (Warmup)
## Chain 3, Iteration: 450 / 4500 [ 10%] (Warmup)
## Chain 3, Iteration: 501 / 4500 [ 11%] (Sampling)
## Chain 3, Iteration: 950 / 4500 [ 21%] (Sampling)
## Chain 3, Iteration: 1400 / 4500 [ 31%] (Sampling)
## Chain 3, Iteration: 1850 / 4500 [ 41%] (Sampling)
## Chain 3, Iteration: 2300 / 4500 [ 51%] (Sampling)
## Chain 3, Iteration: 2750 / 4500 [ 61%] (Sampling)
## Chain 3, Iteration: 3200 / 4500 [ 71%] (Sampling)
## Chain 3, Iteration: 3650 / 4500 [ 81%] (Sampling)
## Chain 3, Iteration: 4100 / 4500 [ 91%] (Sampling)
## Chain 3, Iteration: 4500 / 4500 [100%] (Sampling)
## # Elapsed Time: 0.032 seconds (Warm-up)
## # 0.26 seconds (Sampling)
## # 0.292 seconds (Total)
## Warning in FUN(X[[i]], ...): data with name dept is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name applicant_gender is not numeric
## and not used
##
## SAMPLING FOR MODEL 'admit ~ dbinom(applications, p)' 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
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
## Warning: 使われていないコネクション 14 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 13 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 12 (<-kuma:11085) を閉じます
## [ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
compare( m13.2 , m13.3 , m13.4 )
## WAIC pWAIC dWAIC weight SE dSE
## m13.3 5191.2 11.2 0.0 0.99 57.26 NA
## m13.4 5201.1 6.0 9.9 0.01 56.88 6.79
## m13.2 5201.6 7.0 10.4 0.01 56.92 6.48
## R code 13.22
library(rethinking)
data(chimpanzees)
d <- chimpanzees
d$recipient <- NULL
d$block_id <- d$block
head(d)
## actor condition block trial prosoc_left chose_prosoc pulled_left
## 1 1 0 1 2 0 1 0
## 2 1 0 1 4 0 0 1
## 3 1 0 1 6 1 0 0
## 4 1 0 1 8 0 1 0
## 5 1 0 1 10 1 1 1
## 6 1 0 1 12 1 1 1
## block_id
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
m13.6 <- map2stan(
alist(
# likeliood
pulled_left ~ dbinom(1,p),
# linear models
logit(p) <- A + (BP + BPC*condition)*prosoc_left,
A <- a + a_actor[actor] + a_block[block_id],
BP <- bp + bp_actor[actor] + bp_block[block_id],
BPC <- bpc + bpc_actor[actor] + bpc_block[block_id],
# adaptive priors
c(a_actor,bp_actor,bpc_actor)[actor] ~
dmvnorm2(0,sigma_actor,Rho_actor),
c(a_block,bp_block,bpc_block)[block_id] ~
dmvnorm2(0,sigma_block,Rho_block),
# fixed priors
c(a,bp,bpc) ~ dnorm(0,1),
sigma_actor ~ dcauchy(0,2),
sigma_block ~ dcauchy(0,2),
Rho_actor ~ dlkjcorr(4),
Rho_block ~ dlkjcorr(4)
) , data=d , iter=5000 , warmup=1000 , chains=3 , cores=3 )
## Precompiling model...
## the number of chains is less than 1; sampling not done
## Dispatching 3 chains to 3 cores.
##
## SAMPLING FOR MODEL 'pulled_left ~ dbinom(1, p)' 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
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
## Warning in map2stan(alist(pulled_left ~ dbinom(1, p), logit(p) <- A + (BP + : There were 347 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## R code 13.23
m13.6NC <- map2stan(
alist(
pulled_left ~ dbinom(1,p),
logit(p) <- A + (BP + BPC*condition)*prosoc_left,
A <- a + a_actor[actor] + a_block[block_id],
BP <- bp + bp_actor[actor] + bp_block[block_id],
BPC <- bpc + bpc_actor[actor] + bpc_block[block_id],
# adaptive NON-CENTERED priors
c(a_actor,bp_actor,bpc_actor)[actor] ~
dmvnormNC(sigma_actor,Rho_actor),
c(a_block,bp_block,bpc_block)[block_id] ~
dmvnormNC(sigma_block,Rho_block),
c(a,bp,bpc) ~ dnorm(0,1),
sigma_actor ~ dcauchy(0,2),
sigma_block ~ dcauchy(0,2),
Rho_actor ~ dlkjcorr(4),
Rho_block ~ dlkjcorr(4)
) , data=d , iter=5000 , warmup=1000 , chains=3 , cores=3 )
## Precompiling model...
## the number of chains is less than 1; sampling not done
## Dispatching 3 chains to 3 cores.
## Warning: 使われていないコネクション 14 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 13 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 12 (<-kuma:11085) を閉じます
##
## SAMPLING FOR MODEL 'pulled_left ~ dbinom(1, p)' 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
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
## R code 13.24
# extract n_eff values for each model
neff_c <- precis(m13.6,2)@output$n_eff
## Warning in precis(m13.6, 2): There were 347 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
neff_nc <- precis(m13.6NC,2)@output$n_eff
## Warning: 使われていないコネクション 23 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 22 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 21 (<-kuma:11085) を閉じます
# plot distributions
boxplot( list( 'm13.6'=neff_c , 'm13.6NC'=neff_nc ) ,
ylab="effective samples" , xlab="model" )
## R code 13.25
precis( m13.6NC , depth=2 , pars=c("sigma_actor","sigma_block") )
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## sigma_actor[1] 2.34 0.88 1.17 3.51 2605 1
## sigma_actor[2] 0.46 0.36 0.00 0.90 5542 1
## sigma_actor[3] 0.51 0.47 0.00 1.07 4490 1
## sigma_block[1] 0.23 0.22 0.00 0.46 3629 1
## sigma_block[2] 0.56 0.41 0.00 1.03 3809 1
## sigma_block[3] 0.52 0.42 0.00 1.03 5435 1
## R code 13.26
p <- link(m13.6NC)
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
str(p)
## List of 4
## $ p : num [1:1000, 1:504] 0.272 0.238 0.257 0.205 0.4 ...
## $ A : num [1:1000, 1:504] -0.984 -1.161 -1.064 -1.355 -0.406 ...
## $ BP : num [1:1000, 1:504] 0.176 0.735 1.029 1.524 0.145 ...
## $ BPC: num [1:1000, 1:504] -0.254 -1.165 -0.836 -0.78 0.483 ...
## R code 13.27
#compare( m13.6NC)# , m12.5 )
## R code 13.28
m13.6nc1 <- map2stan(
alist(
pulled_left ~ dbinom(1,p),
# linear models
logit(p) <- A + (BP + BPC*condition)*prosoc_left,
A <- a + za_actor[actor]*sigma_actor[1] +
za_block[block_id]*sigma_block[1],
BP <- bp + zbp_actor[actor]*sigma_actor[2] +
zbp_block[block_id]*sigma_block[2],
BPC <- bpc + zbpc_actor[actor]*sigma_actor[3] +
zbpc_block[block_id]*sigma_block[3],
# adaptive priors
c(za_actor,zbp_actor,zbpc_actor)[actor] ~ dmvnorm(0,Rho_actor),
c(za_block,zbp_block,zbpc_block)[block_id] ~ dmvnorm(0,Rho_block),
# fixed priors
c(a,bp,bpc) ~ dnorm(0,1),
sigma_actor ~ dcauchy(0,2),
sigma_block ~ dcauchy(0,2),
Rho_actor ~ dlkjcorr(4),
Rho_block ~ dlkjcorr(4)
) ,
data=d ,
start=list( sigma_actor=c(1,1,1), sigma_block=c(1,1,1) ),
constraints=list( sigma_actor="lower=0", sigma_block="lower=0" ),
types=list( Rho_actor="corr_matrix", Rho_block="corr_matrix" ),
iter=5000 , warmup=1000 , chains=3 , cores=3 )
## Precompiling model...
## the number of chains is less than 1; sampling not done
## Dispatching 3 chains to 3 cores.
##
## SAMPLING FOR MODEL 'pulled_left ~ dbinom(1, p)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0.01 seconds (Sampling)
## # 0.01 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
## Warning in map2stan(alist(pulled_left ~ dbinom(1, p), logit(p) <- A + (BP + : There were 1 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## R code 13.29
# load the distance matrix
library(rethinking)
data(islandsDistMatrix)
#head(islandsDistMatrix)
# display short column names, so fits on screen
Dmat <- islandsDistMatrix
colnames(Dmat) <- c("Ml","Ti","SC","Ya","Fi","Tr","Ch","Mn","To","Ha")
round(Dmat,1)
## Ml Ti SC Ya Fi Tr Ch Mn To Ha
## Malekula 0.0 0.5 0.6 4.4 1.2 2.0 3.2 2.8 1.9 5.7
## Tikopia 0.5 0.0 0.3 4.2 1.2 2.0 2.9 2.7 2.0 5.3
## Santa Cruz 0.6 0.3 0.0 3.9 1.6 1.7 2.6 2.4 2.3 5.4
## Yap 4.4 4.2 3.9 0.0 5.4 2.5 1.6 1.6 6.1 7.2
## Lau Fiji 1.2 1.2 1.6 5.4 0.0 3.2 4.0 3.9 0.8 4.9
## Trobriand 2.0 2.0 1.7 2.5 3.2 0.0 1.8 0.8 3.9 6.7
## Chuuk 3.2 2.9 2.6 1.6 4.0 1.8 0.0 1.2 4.8 5.8
## Manus 2.8 2.7 2.4 1.6 3.9 0.8 1.2 0.0 4.6 6.7
## Tonga 1.9 2.0 2.3 6.1 0.8 3.9 4.8 4.6 0.0 5.0
## Hawaii 5.7 5.3 5.4 7.2 4.9 6.7 5.8 6.7 5.0 0.0
head(Dmat)
## Ml Ti SC Ya Fi Tr Ch Mn To Ha
## Malekula 0.000 0.475 0.631 4.363 1.234 2.036 3.178 2.794 1.860 5.678
## Tikopia 0.475 0.000 0.315 4.173 1.236 2.007 2.877 2.670 1.965 5.283
## Santa Cruz 0.631 0.315 0.000 3.859 1.550 1.708 2.588 2.356 2.279 5.401
## Yap 4.363 4.173 3.859 0.000 5.391 2.462 1.555 1.616 6.136 7.178
## Lau Fiji 1.234 1.236 1.550 5.391 0.000 3.219 4.027 3.906 0.763 4.884
## Trobriand 2.036 2.007 1.708 2.462 3.219 0.000 1.801 0.850 3.893 6.653
## R code 13.30
# linear
curve( exp(-1*x) , from=0 , to=4 , lty=2 ,
xlab="distance" , ylab="correlation" )
# squared
curve( exp(-1*x^2) , add=TRUE )
## R code 13.31
data(Kline2) # load the ordinary data, now with coordinates
d <- Kline2
d$society <- 1:10 # index observations
head(d)
## culture population contact total_tools mean_TU lat lon lon2
## 1 Malekula 1100 low 13 3.2 -16.3 167.5 -12.5
## 2 Tikopia 1500 low 22 4.7 -12.3 168.8 -11.2
## 3 Santa Cruz 3600 low 24 4.0 -10.7 166.0 -14.0
## 4 Yap 4791 high 43 5.0 9.5 138.1 -41.9
## 5 Lau Fiji 7400 high 33 5.0 -17.7 178.1 -1.9
## 6 Trobriand 8000 high 19 4.0 -8.7 150.9 -29.1
## logpop society
## 1 7.003065 1
## 2 7.313220 2
## 3 8.188689 3
## 4 8.474494 4
## 5 8.909235 5
## 6 8.987197 6
m13.7 <- map2stan(
alist(
total_tools ~ dpois(lambda),
log(lambda) <- a + g[society] + bp*logpop,
g[society] ~ GPL2( Dmat , etasq , rhosq , 0.01 ),
a ~ dnorm(0,10),
bp ~ dnorm(0,1),
etasq ~ dcauchy(0,1),
rhosq ~ dcauchy(0,1)
),
data=list(
total_tools=d$total_tools,
logpop=d$logpop,
society=d$society,
Dmat=islandsDistMatrix),
warmup=2000 , iter=1e4 , chains=4 )
##
## SAMPLING FOR MODEL 'total_tools ~ dpois(lambda)' 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: 2001 / 10000 [ 20%] (Sampling)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Sampling)
## Chain 1, Iteration: 4000 / 10000 [ 40%] (Sampling)
## Chain 1, Iteration: 5000 / 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: 6.995 seconds (Warm-up)
## # 27.852 seconds (Sampling)
## # 34.847 seconds (Total)
##
##
## SAMPLING FOR MODEL 'total_tools ~ dpois(lambda)' 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: 2001 / 10000 [ 20%] (Sampling)
## Chain 2, Iteration: 3000 / 10000 [ 30%] (Sampling)
## Chain 2, Iteration: 4000 / 10000 [ 40%] (Sampling)
## Chain 2, Iteration: 5000 / 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.509 seconds (Warm-up)
## # 29.03 seconds (Sampling)
## # 35.539 seconds (Total)
##
##
## SAMPLING FOR MODEL 'total_tools ~ dpois(lambda)' NOW (CHAIN 3).
##
## Chain 3, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 3, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 3, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 3, Iteration: 2001 / 10000 [ 20%] (Sampling)
## Chain 3, Iteration: 3000 / 10000 [ 30%] (Sampling)
## Chain 3, Iteration: 4000 / 10000 [ 40%] (Sampling)
## Chain 3, Iteration: 5000 / 10000 [ 50%] (Sampling)
## Chain 3, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 3, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 3, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 3, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 3, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 6.51 seconds (Warm-up)
## # 29.824 seconds (Sampling)
## # 36.334 seconds (Total)
##
##
## SAMPLING FOR MODEL 'total_tools ~ dpois(lambda)' NOW (CHAIN 4).
##
## Chain 4, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 4, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 4, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 4, Iteration: 2001 / 10000 [ 20%] (Sampling)
## Chain 4, Iteration: 3000 / 10000 [ 30%] (Sampling)
## Chain 4, Iteration: 4000 / 10000 [ 40%] (Sampling)
## Chain 4, Iteration: 5000 / 10000 [ 50%] (Sampling)
## Chain 4, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 4, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 4, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 4, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 4, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 6.971 seconds (Warm-up)
## # 30.105 seconds (Sampling)
## # 37.076 seconds (Total)
##
##
## SAMPLING FOR MODEL 'total_tools ~ dpois(lambda)' 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
## [ 3200 / 32000 ]
## Warning: 使われていないコネクション 21 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 14 (<-kuma:11085) を閉じます
## Warning: 使われていないコネクション 13 (<-kuma:11085) を閉じます
## [ 6400 / 32000 ]
[ 9600 / 32000 ]
[ 12800 / 32000 ]
[ 16000 / 32000 ]
[ 19200 / 32000 ]
[ 22400 / 32000 ]
[ 25600 / 32000 ]
[ 28800 / 32000 ]
[ 32000 / 32000 ]
## R code 13.32
precis(m13.7,depth=2)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## g[1] -0.27 0.46 -0.98 0.43 2263 1
## g[2] -0.12 0.45 -0.84 0.52 2140 1
## g[3] -0.17 0.44 -0.83 0.47 2151 1
## g[4] 0.29 0.39 -0.27 0.88 2141 1
## g[5] 0.02 0.39 -0.55 0.57 2119 1
## g[6] -0.46 0.39 -1.03 0.09 2279 1
## g[7] 0.09 0.38 -0.47 0.63 2153 1
## g[8] -0.26 0.38 -0.83 0.26 2233 1
## g[9] 0.23 0.36 -0.28 0.74 2204 1
## g[10] -0.12 0.46 -0.82 0.58 3839 1
## a 1.31 1.19 -0.56 3.16 3205 1
## bp 0.25 0.12 0.07 0.43 4296 1
## etasq 0.36 1.01 0.00 0.72 8312 1
## rhosq 1.69 14.05 0.01 2.17 9608 1
## R code 13.33
post <- extract.samples(m13.7)
# plot the posterior median covariance function
curve( median(post$etasq)*exp(-median(post$rhosq)*x^2) , from=0 , to=10 ,
xlab="distance (thousand km)" , ylab="covariance" , ylim=c(0,1) ,
yaxp=c(0,1,4) , lwd=2 )
# plot 100 functions sampled from posterior
for ( i in 1:100 )
curve( post$etasq[i]*exp(-post$rhosq[i]*x^2) , add=TRUE ,
col=col.alpha("black",0.2) )
## R code 13.34
# compute posterior median covariance among societies
K <- matrix(0,nrow=10,ncol=10)
for ( i in 1:10 )
for ( j in 1:10 )
K[i,j] <- median(post$etasq) *
exp( -median(post$rhosq) * islandsDistMatrix[i,j]^2 )
diag(K) <- median(post$etasq) + 0.01
## R code 13.35
# convert to correlation matrix
Rho <- round( cov2cor(K) , 2 )
# add row/col names for convenience
colnames(Rho) <- c("Ml","Ti","SC","Ya","Fi","Tr","Ch","Mn","To","Ha")
rownames(Rho) <- colnames(Rho)
Rho
## Ml Ti SC Ya Fi Tr Ch Mn To Ha
## Ml 1.00 0.87 0.81 0.00 0.52 0.18 0.02 0.04 0.24 0
## Ti 0.87 1.00 0.92 0.00 0.52 0.19 0.04 0.06 0.21 0
## SC 0.81 0.92 1.00 0.00 0.37 0.30 0.07 0.11 0.12 0
## Ya 0.00 0.00 0.00 1.00 0.00 0.09 0.37 0.34 0.00 0
## Fi 0.52 0.52 0.37 0.00 1.00 0.02 0.00 0.00 0.76 0
## Tr 0.18 0.19 0.30 0.09 0.02 1.00 0.26 0.72 0.00 0
## Ch 0.02 0.04 0.07 0.37 0.00 0.26 1.00 0.53 0.00 0
## Mn 0.04 0.06 0.11 0.34 0.00 0.72 0.53 1.00 0.00 0
## To 0.24 0.21 0.12 0.00 0.76 0.00 0.00 0.00 1.00 0
## Ha 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1
## R code 13.36
# scale point size to logpop
psize <- d$logpop / max(d$logpop)
psize <- exp(psize*1.5)-2
# plot raw data and labels
plot( d$lon2 , d$lat , xlab="longitude" , ylab="latitude" ,
col=rangi2 , cex=psize , pch=16 , xlim=c(-50,30) )
labels <- as.character(d$culture)
text( d$lon2 , d$lat , labels=labels , cex=0.7 , pos=c(2,4,3,3,4,1,3,2,4,2) )
# overlay lines shaded by Rho
for( i in 1:10 )
for ( j in 1:10 )
if ( i < j )
lines( c( d$lon2[i],d$lon2[j] ) , c( d$lat[i],d$lat[j] ) ,
lwd=2 , col=col.alpha("black",Rho[i,j]^2) )
## R code 13.37
# compute posterior median relationship, ignoring distance
logpop.seq <- seq( from=6 , to=14 , length.out=30 )
lambda <- sapply( logpop.seq , function(lp) exp( post$a + post$bp*lp ) )
lambda.median <- apply( lambda , 2 , median )
lambda.PI80 <- apply( lambda , 2 , PI , prob=0.8 )
# plot raw data and labels
plot( d$logpop , d$total_tools , col=rangi2 , cex=psize , pch=16 ,
xlab="log population" , ylab="total tools" )
text( d$logpop , d$total_tools , labels=labels , cex=0.7 ,
pos=c(4,3,4,2,2,1,4,4,4,2) )
# display posterior predictions
lines( logpop.seq , lambda.median , lty=2 )
lines( logpop.seq , lambda.PI80[1,] , lty=2 )
lines( logpop.seq , lambda.PI80[2,] , lty=2 )
# overlay correlations
for( i in 1:10 )
for ( j in 1:10 )
if ( i < j )
lines( c( d$logpop[i],d$logpop[j] ) ,
c( d$total_tools[i],d$total_tools[j] ) ,
lwd=2 , col=col.alpha("black",Rho[i,j]^2) )
## R code 13.38
#S <- matrix( c( sa^2 , sa*sb*rho , sa*sb*rho , sb^2 ) , nrow=2 )