Chapter5 練習問題 Part1

(1)model5-3.stan(P.57)をRから実行して得られるMCMCサンプルを使って、それぞれのnにおけるε[n]=Y[n]-μ[n]のMCMCサンプルを算出するRコードを書け

インプットするデータは以下のようなもの（冒頭の１０行ほど）

data-attendance-1.txt

A,Score,Y
0,69,0.286
1,145,0.196
0,125,0.261
1,86,0.109
1,158,0.23
0,133,0.35
0,111,0.33
1,147,0.194
0,146,0.413

まず、model5-3.stanの元となるモデル式は・・・

model5-3.stanは以下のとおり。

data{
  int N;
  int<lower=0,upper=1> A[N];
  real<lower=0,upper=1> Score[N];
  real<lower=0,upper=1> Y[N];
}

parameters{
  real b1;
  real b2;
  real b3;
  real<lower=0> sigma;
}

transformed parameters{
  real myu[N];
  for(n in 1:N)
    myu[n] = b1 + b2*A[n] + b3*Score[n];
}

model{
  for(n in 1:N)
    Y[n] ~ normal(myu[n], sigma);
}

generated quantities{
  real y_pred[N];
  for(n in 1:N)
    y_pred[n] = normal_rng(myu[n],sigma);
}

これをRで実行する。

#データ読み込み
d <- read.csv(file="data-attendance-1.txt")
#データリストの指定
#スコアは２００で割ってスケーリングし、０－１の値にする
data <- list(N=nrow(d), A=d$A, Score=d$Score/200, Y=d$Y)

#モデル作成
fit <- stan(file="model5-3.stan", data=data, seed=1234)

ようやく本題。ε[n]を求めるRコードは・・・

ms <- rstan::extract(fit) #MCMCデータの取り出し
e <- d$Y - ms$myu

これでいいのかな、と思って回答を見に行ったら・・・

noise <- t(-t(ms$mu) + d$Y)

なんだこれは・・・よくわからないので、ひとまず置いておいて先に進む。

(2)model5-3.stanを修正し、Stan内でε[n]を算出せよ

generated quantitiesのところを修正すればいいはず

data{
  int N;
  int<lower=0,upper=1> A[N];
  real<lower=0,upper=1> Score[N];
  real<lower=0,upper=1> Y[N];
}

parameters{
  real b1;
  real b2;
  real b3;
  real<lower=0> sigma;
}

transformed parameters{
  real myu[N];
  for(n in 1:N)
    myu[n] = b1 + b2*A[n] + b3*Score[n];
}

model{
  for(n in 1:N)
    Y[n] ~ normal(myu[n], sigma);
}

generated quantities{
  real y_pred[N];
  real e[N];
  for(n in 1:N){
    y_pred[n] = normal_rng(myu[n],sigma);
    e[n] = Y[n] - myu[n];
  }
}

これを実行して、(1)の２つの計算方法で算出したε[n]の値と、Stan内で算出したε[n]を比較してみる。

fit2 <- stan(file="model5-3-2.stan", data=data, seed=1234)
ms2 <- rstan::extract(fit2)
e2 <- ms2$e
e <- d$Y - ms2$myu
noise <- t(-t(ms2$myu)+d$Y)

結果は

#差の合計でチェック
> sum(e-e2)
[1] 4.551914e-15
> sum(noise-e2)
[1] 0

t関数を組み入れた方がやはり一致。なぜだかよく理解できていないが、こういうものと覚えておくしかない。。。

＜9/5追記＞
ようやくわかった。
t関数は行列の行と列を入れ替える関数（要するに転置行列を作る関数）
今回扱っているデータ件数は50件のもの。
d$Yは行の方向にデータが入っている（50行）けど、ms2$myuは列の方向にデータが入っている（50列）。
ms2$myuの行方向にはMCMCサンプル数分のデータが入っている（ここでは4000行）
よって引き算するために、ms2$myuの転置行列を作らないといけない（1つ目のt関数）
しかし、そのままだと結果が50行×4000列の行列となるので、もう１回t関数をかまして、
4000行×50列の行列にしている。
解釈はこれであってるはず。

(3)以降はまた後日。