【2回目】100日後にベイズ統計ができるゆーみん【20日目】

おはようございます！

今日も少しずつ進めていきます。

5.2.5 Stanで実装

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

parameters {
 real b1;
 real b2;
 real b3;
}

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

model {
 for (n in 1:N)
   Y[n] ~ binomial(M[n], q[n]);
}

generated quantities {
 real y_pred[N];
 for (n in 1:N)
   y_pred[n] = binomial_rng(M[n], q[n]);
}

細かく見ていくと

int<lower=0> Y[N];

Yを出席率のreal型から出席回数のint型に変更します。

q[n] = inv_logit(b1 + b2*A[n] + b3*Score[n]);

inv_logit関数を用いて(0,1)の値に変換しています。

 Y[n] ~ binomial(M[n], q[n]);

y_pred[n] = binomial_rng(M[n], q[n]);

確率的な変動部分を正規分布から二項分布に変更しています。

Rコードの方はほとんど同じなので解説なしでした。

library(rstan)

d <- read.csv(file='input/data-attendance-2.txt')
data <- list(N=nrow(d), A=d$A, Score=d$Score/200, M=d$M, Y=d$Y)
fit <- stan(file='model/model5-4.stan', data=data, seed=1234)

save.image('output/result-model5-4.RData')

短いですが、今日はここまでです。
もうあまりにもゆっくり進めているので、
少し気合を入れ直して頑張ります！

では、また明日お会いしましょう〜