6.2.1 独立な2群の平均値差に関する推測
setwd("~/new/Bayes/Beyes/rstan/rstan/new200/chap6")
source("data621.R")
#scr<-"model621.stan"
data <-list(N1=N1, N2=N2, x1=x1, x2=x2)
data
## $N1
## [1] 20
##
## $N2
## [1] 20
##
## $x1
## [1] 30.86 29.75 31.55 32.29 29.90 31.71 31.35 29.03 30.37 31.55 29.26
## [12] 32.29 29.90 30.18 30.72 32.28 30.72 29.90 31.55 31.55
##
## $x2
## [1] 31.36 33.34 33.16 31.36 36.19 29.80 31.11 35.23 31.36 31.27 31.63
## [12] 31.63 32.00 31.11 31.63 31.36 31.81 31.63 29.21 33.37
model621.stan <- '
data {
int<lower=0> N1;
int<lower=0> N2;
real<lower=0> x1[N1];
real<lower=0> x2[N2];
}
parameters {
real mu1;
real mu2;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
transformed parameters {
real<lower=0> sigma1sq;
real<lower=0> sigma2sq;
sigma1sq <- pow(sigma1,2);
sigma2sq <- pow(sigma2,2);
}
model {
x1 ~ normal(mu1,sigma1);
x2 ~ normal(mu2,sigma2);
}
generated quantities{
real delta;
real delta_over;
real delta_over1;
delta <- mu2 - mu1;
delta_over <- step(delta);
delta_over1 <- if_else(delta>1,1,0);
}
'
par<-c("mu1","mu2","sigma1","sigma2","delta","delta_over",
"delta_over1")
war<-3000
ite<-11000
see<-1234
dig<-3
cha<-2
library(rstan)
## Warning: package 'rstan' was built under R version 3.2.2
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 3.2.2
## Loading required package: inline
##
## Attaching package: 'inline'
##
## The following object is masked from 'package:Rcpp':
##
## registerPlugin
##
## rstan (Version 2.7.0-1, packaged: 2015-07-17 18:12:01 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())
fit <- stan(model_code= model621.stan,
data = data,
iter=ite,
seed=see,
warmup=war,
pars=par,chains=cha)
## COMPILING THE C++ CODE FOR MODEL 'model621.stan' NOW.
##
## SAMPLING FOR MODEL 'model621.stan' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 11000 [ 0%] (Warmup)
## Chain 1, Iteration: 1100 / 11000 [ 10%] (Warmup)
## Chain 1, Iteration: 2200 / 11000 [ 20%] (Warmup)
## Chain 1, Iteration: 3001 / 11000 [ 27%] (Sampling)
## Chain 1, Iteration: 4100 / 11000 [ 37%] (Sampling)
## Chain 1, Iteration: 5200 / 11000 [ 47%] (Sampling)
## Chain 1, Iteration: 6300 / 11000 [ 57%] (Sampling)
## Chain 1, Iteration: 7400 / 11000 [ 67%] (Sampling)
## Chain 1, Iteration: 8500 / 11000 [ 77%] (Sampling)
## Chain 1, Iteration: 9600 / 11000 [ 87%] (Sampling)
## Chain 1, Iteration: 10700 / 11000 [ 97%] (Sampling)
## Chain 1, Iteration: 11000 / 11000 [100%] (Sampling)
## # Elapsed Time: 0.063 seconds (Warm-up)
## # 0.156 seconds (Sampling)
## # 0.219 seconds (Total)
##
##
## SAMPLING FOR MODEL 'model621.stan' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 11000 [ 0%] (Warmup)
## Chain 2, Iteration: 1100 / 11000 [ 10%] (Warmup)
## Chain 2, Iteration: 2200 / 11000 [ 20%] (Warmup)
## Chain 2, Iteration: 3001 / 11000 [ 27%] (Sampling)
## Chain 2, Iteration: 4100 / 11000 [ 37%] (Sampling)
## Chain 2, Iteration: 5200 / 11000 [ 47%] (Sampling)
## Chain 2, Iteration: 6300 / 11000 [ 57%] (Sampling)
## Chain 2, Iteration: 7400 / 11000 [ 67%] (Sampling)
## Chain 2, Iteration: 8500 / 11000 [ 77%] (Sampling)
## Chain 2, Iteration: 9600 / 11000 [ 87%] (Sampling)
## Chain 2, Iteration: 10700 / 11000 [ 97%] (Sampling)
## Chain 2, Iteration: 11000 / 11000 [100%] (Sampling)
## # Elapsed Time: 0.062 seconds (Warm-up)
## # 0.188 seconds (Sampling)
## # 0.25 seconds (Total)
print(fit,pars=par,digits_summary=dig)
## Inference for Stan model: model621.stan.
## 2 chains, each with iter=11000; warmup=3000; thin=1;
## post-warmup draws per chain=8000, total post-warmup draws=16000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
## mu1 30.839 0.003 0.253 30.341 30.676 30.840 30.999 31.331 7986
## mu2 31.978 0.004 0.392 31.208 31.722 31.980 32.233 32.763 8465
## sigma1 1.091 0.002 0.196 0.788 0.952 1.062 1.196 1.550 7221
## sigma2 1.742 0.003 0.309 1.256 1.521 1.702 1.919 2.449 9040
## delta 1.139 0.005 0.466 0.224 0.837 1.142 1.440 2.075 7876
## delta_over 0.992 0.001 0.092 1.000 1.000 1.000 1.000 1.000 10292
## delta_over1 0.621 0.005 0.485 0.000 0.000 1.000 1.000 1.000 10164
## Rhat
## mu1 1.000
## mu2 1.000
## sigma1 1.000
## sigma2 1.001
## delta 1.000
## delta_over 1.000
## delta_over1 1.000
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 15 12:43:25 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)
plot(fit)
pairs(fit, include = FALSE, las = 1)
