アヒル本4章問題のこり

前回に続き、松浦健太郎先生著の「StanとRでベイズ統計モデリング (Wonderful R)」（以下アヒル本）の練習問題をやります。

4章（５）が残っています。

Y1は、平均0、標準偏差5の正規分布に従う標本数30のベクトル、

Y2は、平均1、標準偏差4の正規分布に従う標本数30のベクトルのt検定をせよだそうです。

データ設定

以下コード

set.seed(123)
N1<- 30
N2<- 20
Y1 <-rnorm(n=N1, mean = 0, sd =5)
Y2 <- rnorm (n=N2, mean = 1, sd = 4)

 以上

以上をexercise4-data.Rの名前で保存します。

stanコードを書きます。

全く同じでは、流石につまらないので、ベクトル化してみます。

以下コード

data {
int N1;
int N2;
vector[N1] Y1;
vector[N2] Y2;
}

parameters {
real mu1;
real mu2;
real<lower = 0> sigma1;
real<lower = 0> sigma2;
}

model {
Y1 ~ normal(mu1, sigma1);
Y2 ~ normal(mu2, sigma2);
}

以上をex5v.stanと名付けて保存します。

データセクションの

real Y1[N1];

が実数列の書き方、

ベクトル化すると

vector[N1] Y1;

になります。処理が長い場合にベクトル化のほうがやや高速化するとか。

今回は単純なのであまり変わりません。

デメリットは可読性が下がり、バグが入りやすくなるそうです。

なので最初にforでモデルを書いて動くことを確認したのちに、

ベクトル化してchainを増やしたり、サンプリングを増やしたりするそうです。

つづいてRのコード

library(rstan)

source('~/Rstan/chap4/exercise4-data.R')

data <- list(N1 = N1, N2 = N2, Y1 = Y1, Y2 =Y2)
fit <- stan(file='~/Rstan/chap4/ex5v.stan', data = data, seed = 1234)

save.image('result-ex5v.RData')

ms5 <- rstan::extract(fit)
prob5 <- mean(ms5$mu1 < ms5$mu2)
prob5

以上をex5v.Rと名付けて保存しました

prob5は 0.9335でした

	mean	se_mean	sd	X2.5.	X25.	X50.	X75.	X97.5.	n_eff	Rhat
mu1	-0.24	0.02	0.98	-2.21	-0.86	-0.23	0.42	1.67	3297.23	1.00
mu2	1.59	0.01	0.81	-0.01	1.07	1.61	2.11	3.18	3662.62	1.00
sigma1	5.15	0.01	0.71	3.93	4.65	5.08	5.57	6.75	3217.31	1.00
sigma2	3.62	0.01	0.63	2.62	3.17	3.55	3.98	5.07	3279.02	1.00
lp__	-95.37	0.04	1.51	-99.20	-96.11	-95.02	-94.25	-93.52	1765.64	1.00

結果も収束していますが、問題（４）と同じくmu2が大きいです。

収束をしらべます。

以下Rコード
library(ggmcmc)

load('result-ex5v.RData')

write.table(data.frame(summary(fit)$summary),
file='~/Rstan/chap4/fit-summary5v.txt', sep='\t', quote=FALSE, col.names=NA)

p <- ggs_traceplot(ggs(fit, inc_warmup=TRUE, stan_include_auxiliar=TRUE))
p <- p + theme_bw(base_size=18)
p <- p + labs(color='Chain')
ggsave(p, file='20181026XYZ.png', dpi=300, w=7, h=10)

以上

収束結果を図示します。

いい感じに自己相関なくサンプリングできている気がします。

best regard