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

こんにちは！
今日はポアソン回帰の続きをやっていきます。

5.4.3　Stanで実装

↓モデル式の実装例
λ[n] = exp(b1+b2A[n]+b3Score[n]) n=1,...,N
M[n] ~ Poisson(λ[n]) n=1,...,N

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

parameters {
 real b[3];
}

transformed parameters {
 real lambda[N];
 for (n in 1:N)
   lambda[n] = exp(b[1] + b[2]*A[n] + b[3]*Score[n]);
}

model {
 for (n in 1:N)
   M[n] ~ poisson(lambda[n]);
}

generated quantities {
 int m_pred[N];
 for (n in 1:N)
   m_pred[n] = poisson_rng(lambda[n]);
}

二項ロジスティック回帰のコードから変更した部分を説明してくれています。

int<lower=0> M[N];

応答変数を授業回数Mに変更する。

lambda[n] = exp(b[1] + b[2]*A[n] + b[3]*Score[n]);

exp関数を用いて0以上に変換している。

M[n] ~ poisson(lambda[n]);

 m_pred[n] = poisson_rng(lambda[n]);

確率的な変動の部分をポアソン分布に変更している

ここで異なる方法が書かれていました！

Stanにはpoisson_log(x)という分布関数が用意されており、これはpoisson(exp(x))と等価であるそうです。

poisson_log_rng(x)を使った場合↓

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

parameters {
 real b[3];
}

transformed parameters {
 real lambda[N];
 for (n in 1:N)
   lambda[n] = b[1] + b[2]*A[n] + b[3]*Score[n];
}

model {
 for (n in 1:N)
   M[n] ~ poisson_log(lambda[n]);
}

generated quantities {
 int m_pred[N];
 for (n in 1:N)
   m_pred[n] = poisson_log_rng(lambda[n]);
}

今日はここまでです！

おそらく今日の夕刻ごろまたブログをアップするかと思われます。

またお会いしましょう〜