6.2.2 対応のある2群の平均値差に関する推測
setwd("~/new/Bayes/Beyes/rstan/rstan/new200/chap6")
source("data622.R")
#scr<-"model622.stan"
data <-list(N=N, x=x)
data
## $N
## [1] 20
##
## $x
## [,1] [,2]
## [1,] 6 7
## [2,] 11 11
## [3,] 10 14
## [4,] 13 13
## [5,] 17 16
## [6,] 10 12
## [7,] 10 8
## [8,] 7 15
## [9,] 9 12
## [10,] 1 3
## [11,] 14 17
## [12,] 7 11
## [13,] 7 9
## [14,] 11 11
## [15,] 12 14
## [16,] 12 12
## [17,] 14 11
## [18,] 12 15
## [19,] 7 11
## [20,] 13 17
model622.stan <- '
data {
int<lower=0> N;
vector[2] x[N];
}
parameters {
vector[2] mu;
vector<lower=0>[2] sigma;
real<lower=-1,upper=1> rho;
}
transformed parameters {
vector<lower=0>[2] sigmasq;
matrix[2,2] Sigma;
sigmasq[1] <- pow(sigma[1],2);
sigmasq[2] <- pow(sigma[2],2);
Sigma[1,1] <- sigmasq[1];
Sigma[2,2] <- sigmasq[2];
Sigma[1,2] <- sigma[1]*sigma[2]*rho;
Sigma[2,1] <- sigma[1]*sigma[2]*rho;
}
model {
for(i in 1:N){
x[i] ~ multi_normal(mu,Sigma);
}
}
generated quantities{
real delta;
real delta_over;
real delta_over2;
real rho_over;
real rho_over05;
delta <- mu[2] - mu[1];
delta_over <- step(delta);
delta_over2 <- if_else(delta>2,1,0);
rho_over <- step(rho);
rho_over05 <- if_else(rho>0.5,1,0);
}
'
par<-c("mu","Sigma","rho","delta","delta_over",
"delta_over2","rho_over","rho_over05")
war<-4000
ite<-11000
see<-12345
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 = model622.stan,
data = data,
iter=ite,
seed=see,
warmup=war,
pars=par,chains=cha)
## COMPILING THE C++ CODE FOR MODEL 'model622.stan' NOW.
##
## SAMPLING FOR MODEL 'model622.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: 3300 / 11000 [ 30%] (Warmup)
## Chain 1, Iteration: 4001 / 11000 [ 36%] (Sampling)
## Chain 1, Iteration: 5100 / 11000 [ 46%] (Sampling)
## Chain 1, Iteration: 6200 / 11000 [ 56%] (Sampling)
## Chain 1, Iteration: 7300 / 11000 [ 66%] (Sampling)
## Chain 1, Iteration: 8400 / 11000 [ 76%] (Sampling)
## Chain 1, Iteration: 9500 / 11000 [ 86%] (Sampling)
## Chain 1, Iteration: 10600 / 11000 [ 96%] (Sampling)
## Chain 1, Iteration: 11000 / 11000 [100%] (Sampling)
## # Elapsed Time: 1.647 seconds (Warm-up)
## # 3.581 seconds (Sampling)
## # 5.228 seconds (Total)
##
##
## SAMPLING FOR MODEL 'model622.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: 3300 / 11000 [ 30%] (Warmup)
## Chain 2, Iteration: 4001 / 11000 [ 36%] (Sampling)
## Chain 2, Iteration: 5100 / 11000 [ 46%] (Sampling)
## Chain 2, Iteration: 6200 / 11000 [ 56%] (Sampling)
## Chain 2, Iteration: 7300 / 11000 [ 66%] (Sampling)
## Chain 2, Iteration: 8400 / 11000 [ 76%] (Sampling)
## Chain 2, Iteration: 9500 / 11000 [ 86%] (Sampling)
## Chain 2, Iteration: 10600 / 11000 [ 96%] (Sampling)
## Chain 2, Iteration: 11000 / 11000 [100%] (Sampling)
## # Elapsed Time: 1.616 seconds (Warm-up)
## # 2.896 seconds (Sampling)
## # 4.512 seconds (Total)
print(fit,pars=par,digits_summary=dig)
## Inference for Stan model: model622.stan.
## 2 chains, each with iter=11000; warmup=4000; thin=1;
## post-warmup draws per chain=7000, total post-warmup draws=14000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
## mu[1] 10.139 0.012 0.883 8.373 9.563 10.147 10.709 11.894 5385
## mu[2] 11.934 0.012 0.853 10.247 11.390 11.935 12.479 13.660 5417
## Sigma[1,1] 15.379 0.082 5.672 7.860 11.471 14.196 18.080 29.475 4769
## Sigma[1,2] 10.438 0.075 4.625 4.056 7.293 9.523 12.690 21.926 3809
## Sigma[2,1] 10.438 0.075 4.625 4.056 7.293 9.523 12.690 21.926 3809
## Sigma[2,2] 14.318 0.077 5.326 7.257 10.627 13.216 16.812 27.607 4755
## rho 0.696 0.002 0.124 0.402 0.626 0.715 0.785 0.880 5223
## delta 1.795 0.006 0.668 0.478 1.366 1.794 2.228 3.110 11338
## delta_over 0.993 0.001 0.082 1.000 1.000 1.000 1.000 1.000 7591
## delta_over2 0.376 0.004 0.484 0.000 0.000 0.000 1.000 1.000 14000
## rho_over 1.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 14000
## rho_over05 0.927 0.003 0.260 0.000 1.000 1.000 1.000 1.000 7236
## Rhat
## mu[1] 1.000
## mu[2] 1.000
## Sigma[1,1] 1.001
## Sigma[1,2] 1.001
## Sigma[2,1] 1.001
## Sigma[2,2] 1.000
## rho 1.000
## delta 1.000
## delta_over 1.000
## delta_over2 1.000
## rho_over NaN
## rho_over05 1.000
##
## Samples were drawn using NUTS(diag_e) at Tue Sep 15 13:21: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)
