——————————–
model611.stan <- '
data {
int<lower=0> N; //20
real<lower=0> x[N];
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
for(n in 1:N) {
x[n] ~ normal(mu,sigma);
}
}
generated quantities{ //発生する量
real<lower=0,upper=1> mu_over;
real<lower=0,upper=1> mu_over2;
real es;
real<lower=0,upper=1> es_over;
mu_over <- step(mu-2500);#前期の平均値（既知）
mu_over2 <- step(mu-3000); #boder line
es <-(mu-2500)/sigma;
es_over <- step(es-0.8);
}
'
#--------------------------------
data <-list(N=N, x=x)
data
## $N
## [1] 20
##
## $x
## [1] 3060 2840 1780 3280 3550 2450 2200 3070 2100 4100 3630 3060 3280 1870
## [15] 2980 3120 2150 3830 4300 1880
par<-c("mu","sigma","mu_over","mu_over2","es","es_over")
war<-3000 #バーンイン期間
ite<-10000 #サンプル数
see<-12345 #シード
dig<-3 #有効数字
cha<-2 #チェーンの数
fit <- stan(model_code = model611.stan ,
data = data,
iter=ite,
seed=see,
warmup=war,
pars=par,
chains=cha)
## COMPILING THE C++ CODE FOR MODEL 'model611.stan' NOW.
##
## SAMPLING FOR MODEL 'model611.stan' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1, Iteration: 3001 / 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: 0.187 seconds (Warm-up)
## # 0.158 seconds (Sampling)
## # 0.345 seconds (Total)
##
##
## SAMPLING FOR MODEL 'model611.stan' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2, Iteration: 3001 / 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: 0.203 seconds (Warm-up)
## # 0.172 seconds (Sampling)
## # 0.375 seconds (Total)
print(fit,pars=par,digits_summary=dig)
## Inference for Stan model: model611.stan.
## 2 chains, each with iter=10000; warmup=3000; thin=1;
## post-warmup draws per chain=7000, total post-warmup draws=14000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## mu 2923.231 2.382 188.294 2536.264 2805.806 2924.213 3043.545
## sigma 812.655 1.751 141.267 588.236 712.066 793.353 892.705
## mu_over 0.983 0.002 0.128 1.000 1.000 1.000 1.000
## mu_over2 0.334 0.005 0.472 0.000 0.000 0.000 1.000
## es 0.536 0.003 0.245 0.041 0.375 0.539 0.698
## es_over 0.139 0.004 0.346 0.000 0.000 0.000 0.000
## 97.5% n_eff Rhat
## mu 3290.744 6251 1
## sigma 1140.648 6512 1
## mu_over 1.000 5584 1
## mu_over2 1.000 8514 1
## es 1.016 6599 1
## es_over 1.000 7701 1
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 15 12:18:29 2015.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
traceplot(fit,inc_warmup=F)
#traceplot(fit,inc_warmup=T)
plot(fit)
print(fit, pars=par,
probs=c(.1,.5,.9))
## Inference for Stan model: model611.stan.
## 2 chains, each with iter=10000; warmup=3000; thin=1;
## post-warmup draws per chain=7000, total post-warmup draws=14000.
##
## mean se_mean sd 10% 50% 90% n_eff Rhat
## mu 2923.23 2.38 188.29 2688.89 2924.21 3158.82 6251 1
## sigma 812.66 1.75 141.27 650.83 793.35 998.61 6512 1
## mu_over 0.98 0.00 0.13 1.00 1.00 1.00 5584 1
## mu_over2 0.33 0.01 0.47 0.00 0.00 1.00 8514 1
## es 0.54 0.00 0.24 0.23 0.54 0.84 6599 1
## es_over 0.14 0.00 0.35 0.00 0.00 1.00 7701 1
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 15 12:18:29 2015.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
pairs(fit, include = FALSE, las = 1)
#pairs(fit, include = TRUE, las = 1)
### code chunk number 11: rstan_vignette.Rnw:614-619
###################################################
# all chains combined
summary(do.call(rbind, args = get_sampler_params(
fit, inc_warmup = TRUE)),
digits = 2)
## accept_stat__ stepsize__ treedepth__ n_leapfrog__
## Min. :0.00 Min. :6.1e-05 Min. : 1 Min. : 1.0
## 1st Qu.:0.86 1st Qu.:7.3e-01 1st Qu.: 2 1st Qu.: 3.0
## Median :0.95 Median :8.0e-01 Median : 2 Median : 3.0
## Mean :0.88 Mean :8.4e-01 Mean : 2 Mean : 6.1
## 3rd Qu.:0.99 3rd Qu.:8.0e-01 3rd Qu.: 2 3rd Qu.: 3.0
## Max. :1.00 Max. :2.6e+01 Max. :11 Max. :2047.0
## n_divergent__
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.0016
## 3rd Qu.:0.0000
## Max. :1.0000
# each chain separately
lapply(get_sampler_params(fit, inc_warmup = TRUE),
summary, digits = 2)
## [[1]]
## accept_stat__ stepsize__ treedepth__ n_leapfrog__
## Min. :0.00 Min. :6.1e-05 Min. : 1 Min. : 1.0
## 1st Qu.:0.84 1st Qu.:8.0e-01 1st Qu.: 2 1st Qu.: 3.0
## Median :0.95 Median :8.0e-01 Median : 2 Median : 3.0
## Mean :0.87 Mean :8.6e-01 Mean : 2 Mean : 5.9
## 3rd Qu.:1.00 3rd Qu.:8.0e-01 3rd Qu.: 2 3rd Qu.: 3.0
## Max. :1.00 Max. :2.6e+01 Max. :11 Max. :2047.0
## n_divergent__
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.0018
## 3rd Qu.:0.0000
## Max. :1.0000
##
## [[2]]
## accept_stat__ stepsize__ treedepth__ n_leapfrog__
## Min. :0.00 Min. :8.9e-05 Min. : 1.0 Min. : 1.0
## 1st Qu.:0.87 1st Qu.:7.3e-01 1st Qu.: 2.0 1st Qu.: 3.0
## Median :0.96 Median :7.3e-01 Median : 2.0 Median : 3.0
## Mean :0.89 Mean :8.1e-01 Mean : 2.1 Mean : 6.3
## 3rd Qu.:0.99 3rd Qu.:7.3e-01 3rd Qu.: 2.0 3rd Qu.: 3.0
## Max. :1.00 Max. :2.5e+01 Max. :11.0 Max. :2047.0
## n_divergent__
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.0015
## 3rd Qu.:0.0000
## Max. :1.0000
s <- extract(fit, pars = par, permuted = TRUE)
names(s)
## [1] "mu" "sigma" "mu_over" "mu_over2" "es" "es_over"
dim(s$mu)
## [1] 14000
dim(s$mu_over2)
## [1] 14000
s2 <- extract(fit, pars = "mu", permuted = FALSE)
dim(s2)
## [1] 7000 2 1
dimnames(s2)
## $iterations
## NULL
##
## $chains
## [1] "chain:1" "chain:2"
##
## $parameters
## [1] "mu"