Rstan単回帰練習
初心者です。
伝統的統計は記述統計の発展であり、（手続きを理解できている場合は）データの直感的理解を促進します。
一方、機械学習は的確な予測を与えますが、カーネルや多層化など直感を超える処理を含む印象があります。
さて、ベイズモデルはパラメーターが分布で得られることから、モデリングの「理解」と「予測」を両立することが特徴的だと考えています。
さて、久保拓弥「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・ＭＣＭＣ」と豊田秀樹「基礎からのベイズ統計」を読んだので
松浦健太郎「StanとRでベイズ統計モデリング」（通称アヒル本）を学んでいます。
ゴールドスタンダードの順番です。
書評は多くの方々がされており、自身の実力も不足していることから割愛します。
右も左もわからないので、GitHubの写経を行い、一行ずつ理解していくつもりです。
アヒル本の第4章から行います。
GitHubのアドレス
（１）データ設定
以下コード
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の名前で保存します。
Y1は、平均0、標準偏差5の正規分布に従う標本数30のベクトル、
Y2は、平均1、標準偏差4の正規分布に従う標本数30のベクトルと設定できました。
以下作図コード
source('~/Rstan/chap4/excercise4-data.R')
library(ggplot2)
d1 <- data.frame(group = 1, Y = Y1)
d2 <- data.frame(group = 2, Y = Y2)
d <- rbind(d1, d2)
d$group <- as.factor(d$group)
p <- ggplot(data=d, aes(x = group, y = Y, group = group, col = group))
p <- p + geom_boxplot(outlier.size = 0)
p <- p + geom_point(position=position_jitter(w = 0.4, h = 0), size = 2)
ggsave('file = fig-ex1.png', plot = p, dpi =300, w = 4, h = 3)
以上
初心者なので今まではhist関数で大まかにしか作図してませんでしたが、ggplot2パッケージを進められています。
一つ目のpをプロット
二つ目のpをプロット
三つ目のpをプロット
pというオブジェクトに少しずつ情報を追加していくコンセプトのパッケージのようです。
（２）stanでStudentのt検定
まず、コンソールにlibrary(rstan)を打ち込みます。
このパッケージを読み込まないと、Rstudioがstanの入力を補助してくれませんでした。
以下stanコード
data {
int N1;
int N2;
real Y1[N1];
real Y2[N2];
}
parameters {
real mu1;
real mu2;
real<lower = 0> sigma;
}
model {
for (i in 1:N1)
Y1[i] ~ normal(mu1, sigma);
for (i in 1:N2)
Y2[i] ~ normal(mu2, sigma);
}
以上
上記をex3.stanで保存する。
Rではなくstanなので、コードの最後の行が改行だけであることが必要でした。
また、初心者なので;を書き忘れてrstanに怒られます。
HatenaBlogにコピペするとtabが反映されないので、読みにくいです。
以下R側で実行するコード
library(rstan)
source('~/Rstan/chap4/exercise4-data.R')
data <- list(N1 = N1, N2 = N2, Y1 = Y1, Y2 =Y2)
fit <- stan(file='~/Rstan/chap4/ex3.stan', data = data, seed = 1234)
save.image('result-ex3.RData')
以上
Rhatは収束しているように思えました。
Inference for Stan model: ex3.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu1 -0.23 0.01 0.83 -1.89 -0.79 -0.22 0.34 1.35 3610 1
mu2 1.61 0.02 1.00 -0.33 0.91 1.61 2.28 3.59 3391 1
sigma 4.49 0.01 0.47 3.68 4.16 4.44 4.77 5.53 3446 1
lp__ -97.76 0.03 1.30 -101.21 -98.35 -97.40 -96.83 -96.33 1666 1
どう考えてもExcelに張り付けて表にして張るべきです。
次回に活かします。
（４）mu1 < mu2の確率を求めよ
以下Rコード
library(rstan)
load('result-ex3.RData')
ms <- rstan::extract(fit)
prob <- mean(ms$mu1 < ms$mu2)
以上
GitHub通りではなく、extract関数は多様なパッケージで実装されているんので、パッケージを指定しています。
(豊田本(たしか)の方針に従っています。)
とりあえず、単回帰の練習問題はできました。
(5)はparametersの部分で
real<lower=0> sigama;
を
real<lower=0> sigama1;
real<lower=0> sigama2;
に変え、modelも変更するとできるため略します。
こちらは一応ggmcmcのパッケージでTrace Plotを作図しました。
写経すると、「できない」and「おろそかにしている」部分が明示されるので、身が引き締まります。
best regard
