モデリングにも役立つ確率分布の性質（再生性と共役事前分布）

自分が分かりやすいように, 応用しやすいように疑似Stanコードで書きました。

再生性

再生性を使うとモデルをシンプルに書けることがあり、推定のスピードアップにつながります。 以下で用いられる2つの確率変数x1, x2は互いに独立とします。

正規分布

x1 ~ normal(mu1, s1);
x2 ~ normal(mu2, s2);
# ならば =>
y <- x1 + x2;
y ~ normal(mu1 + mu2, sqrt(square(s1)　+　square(s2)));

対数正規分布の場合は変数の積について再生性が成り立つことになります。

コーシー分布

x1 ~ cauchy(mu1, s1);
x2 ~ cauchy(mu2, s2);
# ならば =>
y <- x1 + x2;
y ~ cauchy(mu1 + mu2, s1　+　s2);

ガンマ分布

x1 ~ gamma(a1, b);
x2 ~ gamma(a2, b);
# ならば =>
y <- x1 + x2;
y ~ gamma(a1 + a2, b);

rate parameterであるbが異なる場合は成り立ちません。a1,a2が半整数である場合はカイ二乗分布に相当し再生性を持ちます。

ポアソン分布

x1 ~ poisson(lambda1);
x2 ~ poisson(lambda2);
# ならば =>
y <- x1 + x2;
y ~ poisson(lambda1 + lambda2);

二項分布

x1 ~ binomial(n1, p);
x2 ~ binomial(n2, p);
# ならば =>
y <- x1 + x2;
y ~ binomial(n1 + n2, p);

確率pが異なる場合は成り立ちません。

---

共役事前分布

共役事前分布を使うと事後分布の形が決まるためパラメータの推定の式がシンプルになることがあり、スピードアップにつながります。主要なところをいくつか紹介します。その他のものはwikipedia参照。

ポアソン分布

# data & prior
int<lower=0> Y[N_obs];  # data
real<lower=0> lambda;   # parameter

for (i in 1:N_obs)
   Y[i] ~ poisson(lambda);
lambda ~ gamma(a, b);   # prior

# ならば posterior =>
lambda ~ gamma(a + sum(Y), b + N_obs);

二項分布

# data & prior
int<lower=0> Y[N_obs];    # data
int<lower=1> N[N_obs];    # data
real<lower=0,upper=1> p;  # parameter

for (i in 1:N_obs)
   Y[i] ~ binomial(N[i], p);
p ~ beta(a, b);           # prior

# ならば posterior =>
p ~ beta(a + sum(Y), b + sum(N) - sum(Y));

N[i]がすべて1ならベルヌーイ分布に相当します。

離散分布（カテゴリカル分布）

# data & prior
int<lower=1,upper=K> Z[N_obs];  # data: category index
simplex[K] theta;               # parameter

for (i in 1:N_obs)
   Z[i] ~ categorical(theta);
theta ~ dirichlet(alpha);       # prior

# ならば posterior =>
for (k in 1:K)
   alpha_post[k] <- alpha[k];
for (i in 1:N_obs)
   alpha_post[Z[i]] <- alpha_post[Z[i]] + 1;
theta ~ dirichlet(alpha_post);

多項分布

# data & prior
int<lower=0> Y[N_obs, K];  # data
simplex[K] theta;          # parameter

for (i in 1:N_obs)
   Y[i] ~  multinomial(theta);
theta ~ dirichlet(alpha);  # prior

# ならば posterior =>
for (k in 1:K)
   alpha_post[k] <- alpha[k];
for (i in 1:N_obs)
   for (k in 1:K)
      alpha_post[k] <- alpha_post[k] + Y[i, k];
theta ~ dirichlet(alpha_post);

Stanのmultinomial()関数の引数はN[i]は省略し、Y[i,]の合計からN[i]が定まります。N[i]がすべて1なら離散分布（カテゴリカル分布）に相当します。