Chapter8
Easy
8E1
- Proposal distribution needs to be symmetric.
8E2
Distribution of proposed parameter values is adjusted to currect parameter values. This is done using particular combinations of priors and likelihood functions that can be analytically solved.
8E3
HMC needs continuous parameters, because it “glides” to the lowest multidimensional point over all the parameters. Discrete parameters cannot be used because there is no slope (like a bal will not roll of an even table)
8E4
Samples here means the number of iterations of the Markov chains, not data points. Markov chains are typically autocorrelated, making sequential samples not independent. n_eff is the number of independent samples.
8E5
Rhat should approach 1 from the top. However, even 1.01 is suspicious.
8E6
Should be uncorrelated, ie wildly going up and down, within reasonable estimates of the parameters. Something like:
plot(y=rnorm(1000,0,1),x=1:1000, type="l")
Medium
8M1
## From the book:
## R code 8.2
library(rethinking)
data(rugged)
d <- rugged
d$log_gdp <- log(d$rgdppc_2000)
dd <- d[ complete.cases(d$rgdppc_2000) , ]
## 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, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc6ba3357.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]pairs(m8.1stan)
precis(m8.1stan)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.22 0.14 8.99 9.44 4443 1
## bR -0.20 0.08 -0.33 -0.08 4466 1
## bA -1.95 0.23 -2.31 -1.58 4568 1
## bAR 0.39 0.13 0.18 0.61 4395 1
## sigma 0.95 0.05 0.87 1.03 5850 1With a uniform prior:
m8.1stan.unif <- 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 ~ dunif(0,10)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc64046010.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]precis(m8.1stan.unif)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.22 0.14 9.00 9.45 3639 1
## bR -0.20 0.08 -0.33 -0.08 3540 1
## bA -1.95 0.23 -2.33 -1.60 3721 1
## bAR 0.39 0.13 0.17 0.60 3831 1
## sigma 0.95 0.05 0.86 1.03 6596 1With a dexp prior:
m8.1stan.exp <- 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 ~ dexp(1)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc11b6b88.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0.001 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]precis(m8.1stan.exp)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.22 0.14 9.00 9.45 3499 1
## bR -0.20 0.08 -0.32 -0.08 3819 1
## bA -1.95 0.23 -2.33 -1.60 4063 1
## bAR 0.39 0.13 0.19 0.62 3950 1
## sigma 0.95 0.05 0.86 1.03 6732 1There doesn’t seem to be a difference, the data overwhelms the prior.
8M2
Cauchy
m8.1stan.cauchy1 <- 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,1)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc11625745.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]post.cauchy1<-extract.samples(m8.1stan.cauchy1)m8.1stan.cauchy01 <- 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,0.1)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc7761184.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]post.cauchy01<-extract.samples(m8.1stan.cauchy01)m8.1stan.cauchy001 <- 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,0.01)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc460526f4.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]post.cauchy001<-extract.samples(m8.1stan.cauchy01)dens(post.cauchy1$a, main="cauchy")
dens(post.cauchy01$a, add=T, col="red")
dens(post.cauchy001$a, add=T, col="orange")
Exponential
m8.1stan.exp1 <- 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 ~ dexp(1)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc49da1d48.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]post.exp1<-extract.samples(m8.1stan.exp1)
m8.1stan.exp01 <- 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 ~ dexp(0.1)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc3d5f43b3.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]post.exp01<-extract.samples(m8.1stan.exp01)
m8.1stan.exp001 <- 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 ~ dexp(0.01)
) ,
data=dd.trim, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc19fa46.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]post.exp001<-extract.samples(m8.1stan.exp01)
dens(post.exp1$a, main="exponential")
dens(post.exp01$a, add=T, col="red")
dens(post.exp001$a, add=T, col="orange")
It makes no difference in either case
8M3
m8.1stan.10 <- 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, warmup=10, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc5264e6a.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1197 / 11970 ]
[ 2394 / 11970 ]
[ 3591 / 11970 ]
[ 4788 / 11970 ]
[ 5985 / 11970 ]
[ 7182 / 11970 ]
[ 8379 / 11970 ]
[ 9576 / 11970 ]
[ 10773 / 11970 ]
[ 11970 / 11970 ]m8.1stan.20 <- 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, warmup=20, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc7e73264f.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1194 / 11940 ]
[ 2388 / 11940 ]
[ 3582 / 11940 ]
[ 4776 / 11940 ]
[ 5970 / 11940 ]
[ 7164 / 11940 ]
[ 8358 / 11940 ]
[ 9552 / 11940 ]
[ 10746 / 11940 ]
[ 11940 / 11940 ]m8.1stan.50 <- 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, warmup=50, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc37c94cca.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1185 / 11850 ]
[ 2370 / 11850 ]
[ 3555 / 11850 ]
[ 4740 / 11850 ]
[ 5925 / 11850 ]
[ 7110 / 11850 ]
[ 8295 / 11850 ]
[ 9480 / 11850 ]
[ 10665 / 11850 ]
[ 11850 / 11850 ]m8.1stan.100 <- 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, warmup=100, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc71575d3d.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1170 / 11700 ]
[ 2340 / 11700 ]
[ 3510 / 11700 ]
[ 4680 / 11700 ]
[ 5850 / 11700 ]
[ 7020 / 11700 ]
[ 8190 / 11700 ]
[ 9360 / 11700 ]
[ 10530 / 11700 ]
[ 11700 / 11700 ]m8.1stan.200 <- 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, warmup=200, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc7d6b151.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1140 / 11400 ]
[ 2280 / 11400 ]
[ 3420 / 11400 ]
[ 4560 / 11400 ]
[ 5700 / 11400 ]
[ 6840 / 11400 ]
[ 7980 / 11400 ]
[ 9120 / 11400 ]
[ 10260 / 11400 ]
[ 11400 / 11400 ]m8.1stan.500 <- 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, warmup=500, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc26343f14.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1050 / 10500 ]
[ 2100 / 10500 ]
[ 3150 / 10500 ]
[ 4200 / 10500 ]
[ 5250 / 10500 ]
[ 6300 / 10500 ]
[ 7350 / 10500 ]
[ 8400 / 10500 ]
[ 9450 / 10500 ]
[ 10500 / 10500 ]m8.1stan.1000 <- 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, warmup=1000, iter = 4000, chains=3, cores = 3 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc3dc8e41.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 900 / 9000 ]
[ 1800 / 9000 ]
[ 2700 / 9000 ]
[ 3600 / 9000 ]
[ 4500 / 9000 ]
[ 5400 / 9000 ]
[ 6300 / 9000 ]
[ 7200 / 9000 ]
[ 8100 / 9000 ]
[ 9000 / 9000 ]neff <- data.frame(precis(m8.1stan.10)@output[,"n_eff"],
precis(m8.1stan.20)@output[,"n_eff"],
precis(m8.1stan.50)@output[,"n_eff"],
precis(m8.1stan.100)@output[,"n_eff"],
precis(m8.1stan.200)@output[,"n_eff"],
precis(m8.1stan.500)@output[,"n_eff"],
precis(m8.1stan.1000)@output[,"n_eff"])
plot(y=neff[1,], x=c(10,20,50,100,200,500,1000), type="l", col=1, ylim=c(min(neff), max(neff)), log="x")
for(i in 2:nrow(neff)){lines(y=neff[i,], x=c(10,20,50,100,200,500,1000), type="l", col=i)}
abline(h=4000)
legend("topleft", legend=row.names(precis(m8.1stan.10)@output), fill=1:nrow(neff))
Even low warmup gives a good nr of effective samples. 50 warmup even gives almost the same posterior distributions as 1000 warmups.
precis(m8.1stan.50)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.22 0.14 8.99 9.44 4182 1
## bR -0.20 0.08 -0.33 -0.08 4938 1
## bA -1.94 0.23 -2.30 -1.55 2836 1
## bAR 0.39 0.13 0.19 0.61 3978 1
## sigma 0.95 0.05 0.86 1.03 11850 1precis(m8.1stan.1000)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.22 0.14 9.00 9.45 3560 1
## bR -0.20 0.08 -0.33 -0.08 3397 1
## bA -1.95 0.23 -2.32 -1.60 3948 1
## bAR 0.39 0.13 0.19 0.60 4037 1
## sigma 0.95 0.05 0.87 1.03 6453 1Hard
8H1
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 )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dcf74701c.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'a ~ dnorm(0, 1)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: There aren't enough warmup iterations to fit the
## three stages of adaptation as currently configured.
## Reducing each adaptation stage to 15%/75%/10% of
## the given number of warmup iterations:
## init_buffer = 15
## adapt_window = 75
## term_buffer = 10
##
## Iteration: 1 / 10000 [ 0%] (Warmup)
## Iteration: 101 / 10000 [ 1%] (Sampling)
## Iteration: 1100 / 10000 [ 11%] (Sampling)
## Iteration: 2100 / 10000 [ 21%] (Sampling)
## Iteration: 3100 / 10000 [ 31%] (Sampling)
## Iteration: 4100 / 10000 [ 41%] (Sampling)
## Iteration: 5100 / 10000 [ 51%] (Sampling)
## Iteration: 6100 / 10000 [ 61%] (Sampling)
## Iteration: 7100 / 10000 [ 71%] (Sampling)
## Iteration: 8100 / 10000 [ 81%] (Sampling)
## Iteration: 9100 / 10000 [ 91%] (Sampling)
## Iteration: 10000 / 10000 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.503 seconds (Sampling)
## 0.503 seconds (Total)
##
##
## SAMPLING FOR MODEL 'a ~ dnorm(0, 1)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)precis(mp)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 0.00 1.00 -1.53 1.63 9900 1
## b -0.34 21.14 -5.54 5.66 2303 1show(mp)## map2stan model fit
## 9900 samples from 1 chain
##
## Formula:
## a ~ dnorm(0, 1)
## b ~ dcauchy(0, 1)
##
## Log-likelihood at expected values: 0
## Deviance: 0
## DIC: 0
## Effective number of parameters (pD): 0plot(mp)
pairs(mp)
essentially, this model just finds back the priors. The Cauchy distribution flattens out more, giving a higher density to extreme values?
8H2
Original models:
library(rethinking)
data(WaffleDivorce)
d <- WaffleDivorce
# standardize predictor
d$MedianAgeMarriage.s <- (d$MedianAgeMarriage-mean(d$MedianAgeMarriage))/
sd(d$MedianAgeMarriage)
d$Marriage.s <- (d$Marriage - mean(d$Marriage))/sd(d$Marriage)
m5.1 <- map(
alist(
Divorce ~ dnorm( mu , sigma ) ,
mu <- a + bA * MedianAgeMarriage.s ,
a ~ dnorm( 10 , 10 ) ,
bA ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) , data = d )
m5.2 <- map(
alist(
Divorce ~ dnorm( mu , sigma ) ,
mu <- a + bR * Marriage.s ,
a ~ dnorm( 10 , 10 ) ,
bR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) , data = d )
m5.3 <- map(
alist(
Divorce ~ dnorm( mu , sigma ) ,
mu <- a + bR*Marriage.s + bA*MedianAgeMarriage.s ,
a ~ dnorm( 10 , 10 ) ,
bR ~ dnorm( 0 , 1 ) ,
bA ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data = d )Using MCMC:
m5.1.mcmc <- map2stan(
alist(
Divorce ~ dnorm( mu , sigma ) ,
mu <- a + bA * MedianAgeMarriage.s ,
a ~ dnorm( 10 , 10 ) ,
bA ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) , data = d[,c("Divorce","MedianAgeMarriage.s")] , warmup=500, iter = 4000, chains=3, cores = 3)## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc10fe4b0f.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'Divorce ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1050 / 10500 ]
[ 2100 / 10500 ]
[ 3150 / 10500 ]
[ 4200 / 10500 ]
[ 5250 / 10500 ]
[ 6300 / 10500 ]
[ 7350 / 10500 ]
[ 8400 / 10500 ]
[ 9450 / 10500 ]
[ 10500 / 10500 ]m5.2.mcmc <- map2stan(
alist(
Divorce ~ dnorm( mu , sigma ) ,
mu <- a + bR * Marriage.s ,
a ~ dnorm( 10 , 10 ) ,
bR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) , data = d[,c("Divorce","Marriage.s")] , warmup=500, iter = 4000, chains=3, cores = 3)## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc7f0561f.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'Divorce ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1050 / 10500 ]
[ 2100 / 10500 ]
[ 3150 / 10500 ]
[ 4200 / 10500 ]
[ 5250 / 10500 ]
[ 6300 / 10500 ]
[ 7350 / 10500 ]
[ 8400 / 10500 ]
[ 9450 / 10500 ]
[ 10500 / 10500 ]m5.3.mcmc <- map2stan(
alist(
Divorce ~ dnorm( mu , sigma ) ,
mu <- a + bR*Marriage.s + bA*MedianAgeMarriage.s ,
a ~ dnorm( 10 , 10 ) ,
bR ~ dnorm( 0 , 1 ) ,
bA ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data = d[,c("Divorce","Marriage.s","MedianAgeMarriage.s")] , warmup=500, iter = 4000, chains=3, cores = 3)## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc6c12ae4.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'Divorce ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 1050 / 10500 ]
[ 2100 / 10500 ]
[ 3150 / 10500 ]
[ 4200 / 10500 ]
[ 5250 / 10500 ]
[ 6300 / 10500 ]
[ 7350 / 10500 ]
[ 8400 / 10500 ]
[ 9450 / 10500 ]
[ 10500 / 10500 ]note: the warning about dots refers to the variable name, not the actual values.
compare(m5.1,m5.2,m5.3,m5.1.mcmc,m5.2.mcmc,m5.3.mcmc)## WAIC pWAIC dWAIC weight SE dSE
## m5.1.mcmc 186.3 3.9 0.0 0.44 12.50 NA
## m5.1 187.2 4.5 0.9 0.28 13.91 1.44
## m5.3.mcmc 188.0 4.9 1.8 0.18 12.50 0.79
## m5.3 189.4 5.9 3.1 0.09 14.49 2.22
## m5.2.mcmc 199.9 3.2 13.7 0.00 9.65 9.13
## m5.2 200.6 3.8 14.3 0.00 10.72 9.30So, compare shouldn’t be used this way? In any case, it seems the stan versions of the models perform better than the map versions.
8H3
## 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)
## 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) )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc2bec50f8.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 7.452 seconds (Warm-up)
## 10.372 seconds (Sampling)
## 17.824 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 7.429 seconds (Warm-up)
## 9.943 seconds (Sampling)
## 17.372 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 6.956 seconds (Warm-up)
## 11.092 seconds (Sampling)
## 18.048 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 7.862 seconds (Warm-up)
## 10.304 seconds (Sampling)
## 18.166 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]pairs(m5.8s)
## 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) )## In file included from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/paule/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file28dc106549ba.cpp:8:
## C:/Users/paule/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 2.281 seconds (Warm-up)
## 2.503 seconds (Sampling)
## 4.784 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 1.861 seconds (Warm-up)
## 2.254 seconds (Sampling)
## 4.115 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 3).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 1.836 seconds (Warm-up)
## 2.021 seconds (Sampling)
## 3.857 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 4).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 1.867 seconds (Warm-up)
## 2.565 seconds (Sampling)
## 4.432 seconds (Total)
##
##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
## Iteration: 1 / 1 [100%] (Sampling)
##
## Elapsed Time: 0 seconds (Warm-up)
## 0 seconds (Sampling)
## 0 seconds (Total)
##
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]pairs(m5.8s2)
To maintain the anticorrelation, Bl now has a mirrored distribution compared to Br. Their sum remains constant (around 2).
8H4
compare(m5.8s, m5.8s2)## WAIC pWAIC dWAIC weight SE dSE
## m5.8s 174.7 3.4 0.0 0.68 10.63 NA
## m5.8s2 176.3 2.4 1.5 0.32 10.41 3.92Effective number of parameters is pWAIC (or pD in DIC). In the second model, the effective number of parameters is lower because the priors are more constraining.
8H5
num_weeks <- 1e5
positions <- rep(0,num_weeks)
set.seed(1)
islands <- data.frame(id=1:10,pop=sample(1:10,10,replace=F))
current <- 10
for ( i in 1:num_weeks ) {
# record current position
positions[i] <- current
# flip coin to generate proposal, a random island
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?
# calculate prop with the populations, not the IDs
prob_move <- islands$pop[proposal]/islands$pop[current]
current <- ifelse( runif(1) < prob_move , proposal , current )
}
plot(y=positions[1:100],x=1:100, type="l")
hist(positions, breaks=0:10)
table(positions)/islands$pop## positions
## 1 2 3 4 5 6 7 8
## 1743.000 1705.500 1711.000 1728.571 1794.500 1877.125 1863.333 1859.667
## 9 10
## 1885.200 1908.000