ベイズ統計_13_混合分布モデル

ベイズ統計 13 混合分布モデル 分寺 杏介 神戸大学大学院 経営学研究科  [email protected] ※本スライドは，クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス（CC BY-NC 4.0）に従って利用が可能です。

• https://creativecommons.org/licenses/by-nc/4.0/

１ 分布を混ぜる より柔軟な確率分布の表現 13 混合分布モデル 2

いきなり本題 問 あなたはあるスーパーのマーケターとして，グミの購買について調べることにしました。 その第一歩として，そもそもどれくらいグミが購入されているかを確認したいと思います。 あくまでも仮の状況なので，「どういう人がグミを よく買うんだろう」みたいな話には展開させません。 【データの読み込み】 この名前はなんでもOK "data_pos.csv"をワーキングディレクトリに配置して dat <- read.csv("data_pos.csv") ▌ 中身はこんなデータ ID: POSデータ内でのID（特に使わない） n_buy: 各顧客が先週購入したグミの個数 こいつの分布を調べたい！ ※データは適当に作ったので、実際とは異なります。 13 混合分布モデル 3

とりあえず分布を見てみると… 問 あなたはあるスーパーのマーケターとして，グミの購買について調べることにしました。 その第一歩として，そもそもどれくらいグミが購入されているかを確認したいと思います。 あくまでも仮の状況なので，「どういう人がグミを よく買うんだろう」みたいな話には展開させません。 ▌とりあえず分布の形状を見てみると 【棒グラフ】 【カーネル密度】 barplot(table(dat$n_buy)) plot(density(dat$n_buy)) これは何分布で 表したら良さそう？ 13 混合分布モデル 4

先に種明かし ▌基本的にはポアソン分布でよい 「何個買ったか」というカウントデータなので ▌ただし複数のポアソン分布からデータを混ぜています ▶ イメージとしては グミのヘビーユーザー グミのライトユーザー 2つのデータが 混ざっている状態 ▌誰がどのグループに所属しているかは分からない 今回の話は顧客等のセグメンテーション にも繋がっていく話です たぶん ライト 混ざってる？ たぶん ヘビー …が，データからある程度推測はできる▶ 13 混合分布モデル 5

問題の整理 問 あなたはあるスーパーのマーケターとして，グミの購買について調べることにしました。 その第一歩として，そもそもどれくらいグミが購入されているかを確認したいと思います。 あくまでも仮の状況なので，「どういう人がグミを よく買うんだろう」みたいな話には展開させません。 ▌分布から，どうやら2つのクラスターがありそうだ ▶ 各クラスターに所属する人の平均購買数を調べよう ▌定式化 ライト／ヘビーのクラスターはそれぞれ添字𝐿/𝐻で表します 各顧客𝑖がどちらのクラスターに 所属しているかを表す 潜在変数𝑧𝑖 を用意する 0 𝑧𝑖 = ቊ 1 if 𝑖 ∈ 𝐿 if 𝑖 ∈ 𝐻 この潜在変数𝑧𝑖 の値に応じて 個人が従う確率分布が異なる 𝑦𝑖 ∼ ቊ 𝑃𝑜𝑖𝑠(𝜆𝐿 ) if 𝑧𝑖 = 0 𝑃𝑜𝑖𝑠(𝜆𝐻 ) if 𝑧𝑖 = 1 𝜆𝐿 , 𝜆𝐻 そして各個人の𝑧𝑖 を それぞれ推定したい！ n_buy 13 混合分布モデル 6

あとはstanコードを書くだけ
▌ポアソン分布のパラメータ推定モデルをもとに修正していきます
data {

資料05
p. 46より

int N;
array[N] int N_BUY;
}

parameters {

【完成予定のplate notation】

real <lower=0> lambda;
}

model {

各個人がどっちの𝜆になるかは
まだ分からないので
両方とも矢印をさしておきます

lambda ~ gamma(0.1,0.01);

𝜆𝐿
𝑧𝑖

N_BUY ~ poisson(lambda);

𝜆𝐻
𝑥𝑖
𝑖 data

}
13 混合分布モデル

7

parametersブロックの修正 ▌ポアソン分布のパラメータ推定モデルをもとに修正していきます 1. 2つのlambdaが必要 2. 潜在変数 𝑧𝑖 がデータの数だけ必要 parameters { parameters { real <lower=0> lambda; vector<lower=0> [2] lambda; } array[N] int z; } 【完成予定のplate notation】 𝜆𝐿 ポイント1に関しては，以下の書き方でも良いのですが いくつかの理由からまとめて扱ったほうが何かと都合が良いのです real<lower=0> lambda_L; 𝑧𝑖 real<lower=0> lambda_H; 𝜆𝐻 𝑥𝑖 𝑖 data 13 混合分布モデル 8

modelブロックの修正｜尤度 ▌直感的に実装するなら… 訳あって，一旦コードを書かずにしばらく見ていてください 1. まず潜在変数 𝑧𝑖 をサンプリング 0 𝑧𝑖 = ቊ 1 if 𝑖 ∈ 𝐿 if 𝑖 ∈ 𝐻 ▶ 「0か1」しか取らないので，ベルヌーイ分布が使えそう 各個人のクラスター所属確率 とみなすことができる ▶ そうなると，ベルヌーイ分布の確率を表すパラメータが必要！ 𝑧𝑖 ∼ 𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝜋𝑖 ) 2. 潜在変数 𝑧𝑖 の値に応じてn_buyの尤度を変える ▶ stanにもif構文があるので，それで対応できそう model { // 事前分布はあとで書く X ~ poisson(lambda); } model { // 事前分布はあとでかく z ~ bernoulli(pi); for(i in 1:N){ if(z[i] == 0){ N_BUY[i] ~ poisson(lambda[1]); } else { N_BUY[i] ~ poisson(lambda[2]); } } } 13 混合分布モデル 【完成予定のplate notation】 𝜋𝑖 𝑧𝑖 𝜆𝐿 𝜆𝐻 𝑥𝑖 𝑖 data 9

stanコードの完成？
▌あとは対応する事前分布を書き足せば

まだ見ていてください

model_mixture.stan
data {

model {

int N;

lambda ~ cauchy(0, 10);

array[N] int N_BUY;

pi ~ beta(1,1);

}

z ~ bernoulli(pi);

parameters {

for(i in 1:N){

vector<lower=0> [2] lambda;

if(z[i] == 0){

vector<lower=0, upper=1> [N] pi;

N_BUY[i] ~ poisson(lambda[1]);

array[N] int z;

【plate notation】

} else {

}

N_BUY[i] ~ poisson(lambda[2]);
}

𝜋𝑖

𝜆𝐿

𝜆𝐻

}

𝑧𝑖

}

𝑥𝑖
𝑖 data

13 混合分布モデル

10

stanコードの完成？
▌早速コンパイルしてMCMCじゃい！
library(cmdstanr)
model <- cmdstan_model("model_mixture.stan")

ちゃんと見ていてくれましたか？

Error in stanc(filename, allow_undefined = TRUE) : 0
Semantic error in 'string', line 7, column 2 to column 35:
------------------------------------------------5: parameters {
6:
vector<lower=0> [2] lambda;
7:
array[N] int<lower=0, upper=1> z;
^
8: }
9: model {
------------------------------------------------(Transformed) Parameters cannot be integers.
stanでは，整数（離散変数）がパラメータになることはできない
13 混合分布モデル

これが今日の本題です

11

なぜ離散パラメータが使えない？ ▌HMCのしくみを考えれば当然のこと 2変数のケースを例に考えてみると… 連続変数であれば，（対数）尤度関数は 全域が連続しているお椀のような形 ▶ どう転がしても𝑛秒後の位置は計算可能 離散変数があると，（対数）尤度関数は 完全に分離してしまう ▶ 離散変数をまたいだ計算は不可能 13 混合分布モデル 12

stanで離散パラメータを扱うためには ▌離散パラメータを周辺化する 𝑃 𝑦𝑖 z𝑖 𝑃𝑜𝑖𝑠 𝑦𝑖 𝜆𝐿 =ቊ 𝑃𝑜𝑖𝑠 𝑦𝑖 𝜆𝐻 if 𝑧𝑖 = 0 if 𝑧𝑖 = 1 𝑧𝑖 ∼ 𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝜋𝑖 ) 周辺化 1 𝑃 𝑦𝑖 = ෍ 𝑃 𝑦𝑖 𝑧𝑖 𝑃(𝑧𝑖 = 𝑘) 𝑘=0 𝑃 𝑧𝑖 = 1 = 𝜋𝑖 【尤度以外の部分の修正】 ▶ もう 𝑧 が不要，というだけ 𝑃 𝑧𝑖 = 0 = 1 − 𝜋𝑖 𝑃 𝑦𝑖 = 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐻 𝜋𝑖 + 𝑃𝑜𝑖𝑠 𝑦𝑖 𝜆𝐿 (1 − 𝜋𝑖 ) 周辺化によって離散パラメータ 𝑧𝑖 を 迂回したサンプリングが可能となる！ parameters { vector<lower=0> [2] lambda; vector<lower=0, upper=1> [N] pi; array[N] int z; } データ 𝑦𝑖 の尤度さえ 書けたら良いのです model { lambda ~ cauchy(0, 10); pi ~ beta(1,1); z ~ bernoulli(pi); // 尤度はどう書けば良い？ 【周辺化後のplate notation】 } 13 混合分布モデル 𝜋𝑖 𝑧𝑖 𝜆𝐿 𝜆𝐻 𝑥𝑖 𝑖 data 13

周辺化後の尤度の書き方 ▌stanでは対数スケールで計算していることを意識しながら 𝑃 𝑦𝑖 = 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐻 𝜋𝑖 + 𝑃𝑜𝑖𝑠 𝑦𝑖 𝜆𝐿 (1 − 𝜋𝑖 ) まず対数に log 𝑃 𝑦𝑖 = log 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐻 𝜋𝑖 + 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐿 1 − 𝜋𝑖 ▌stanコードと実際の計算の関係を思い出すと… target += poisson_lpmf(y|lambda) y ~ poisson(lambda) いずれも log𝑃 𝑦𝑖 = log 𝑃𝑜𝑖𝑠(𝑦|𝜆)を計算していた stanには 𝑃𝑜𝑖𝑠(𝑦|𝜆)を直接計算する関数が無い 𝑃𝑜𝑖𝑠 𝑦 𝜆 = exp log 𝑃𝑜𝑖𝑠 𝑦 𝜆 で求める必要がある！ ※ 混合分布の場合sampling記法ではなく，target記法で直接対数尤度を示す必要がある 【できる】 target += exp(poisson_lpmf(y|lambda)) 【できない】y ~ exp(poisson(lambda)) 13 混合分布モデル 14

結局周辺尤度の書き方は？ log 𝑃 𝑦𝑖 = log 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐻 𝜋𝑖 + 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐿 1 − 𝜋𝑖 = log exp log 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐻 𝜋𝑖 + exp log 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐿 1 − 𝜋𝑖 = log exp log 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐻 + log 𝜋𝑖 + exp log 𝑃𝑜𝑖𝑠 𝑦𝑖 |𝜆𝐿 + log 1 − 𝜋𝑖 poisson_lpmf(N_BUY[i]|lambda[2]) poisson_lpmf(N_BUY[i]|lambda[1]) ▌これを素直に書けば… model { for(i in 1:N){ logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]); 1行が長くなってしまうので 分けて書いてみました logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]); target += log(exp(logp1) + exp(logp0)); } これでも良いのですが，まだ少し問題があるのです… } 13 混合分布モデル 15

expは結構怖い model { for(i in 1:N){ logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]); logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]); target += log(exp(logp1) + exp(logp0)); } } ▌場合によっては計算ができないことが 対数確率をexp(𝑥)で変換するということは，元の確率のスケールに戻すということ ▶ あまりに小さい確率の場合，桁あふれを起こす可能性がある 【例】log 𝑃(𝑥) = −700の場合，𝑃 𝑥 = exp log 𝑃(𝑥) ≃ 9.85 × 10−305 対数スケールなら3桁 【Rで試すと】 もとのスケールだと305桁以上 (0.0000…000985…) exp(-700) [1] 9.859677e-305 exp(-800) [1] 0 これが桁あふれ 13 混合分布モデル 305個のゼロ log 0 = −∞なので計算不能に 16

stanは準備万端です model { for(i in 1:N){ logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]); logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]); target += log(exp(logp1) + exp(logp0)); } } ▌不用意にexp⇄logのスケールを行き来しないですむような関数がすでにある log_sum_exp(a,b)= log exp 𝑎 + exp 𝑏 model { これを使えば for(i in 1:N){ logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]); logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]); target += log_sum_exp(logp1, logp0); } } 13 混合分布モデル 17

ようやく推定できるstanコードが完成
▌全体を整えてあげたら
model_mixture.stan
data {

model {

int N;

lambda ~ cauchy(0, 10);

array[N] int N_BUY;

pi ~ beta(1,1);

}

real logp0;
real logp1;

parameters {
vector<lower=0> [2] lambda;

計算途中の生成物は保存が不要なので
modelブロック内で宣言するのがおすすめ

for(i in 1:N){

vector<lower=0, upper=1> [N] pi;

logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]);

}

logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]);
target += log_sum_exp(logp1, logp0);
}

}

13 混合分布モデル

18

推定しましょう ▌推定 model <- cmdstan_model("model_mixture.stan") stan_data <- list(N=nrow(dat), N_BUY=dat$n_buy) 混合分布モデルもそこそこ時間がかかるので 並列化しておくのがおすすめ result <- model$sample(data=stan_data, parallel_chains = 4) ▌結果を見る result$summary() 人によっては（何回も繰り返すと）全然違う結果になっている可能性もあります 各クラスターの 平均購入数 各顧客が 「ヘビーユーザー」 に所属する確率 13 混合分布モデル 19

あまりにもうまくいっていない 各クラスターの 平均購入数 各顧客が 「ヘビーユーザー」 に所属する確率 ▌クラスターが全く分離できていない？ ▶ その原因はlambdaのトレースプロットに 𝜆𝐻 っぽい bayesplot::mcmc_trace(result$draws("lambda")) chainによってlambda[1]とlambda[2]の どちらが𝜆ℎ , 𝜆𝐿 にあたるかが異なっていた 𝜆𝐿 っぽい ラベルスイッチングの問題と呼ばれます 13 混合分布モデル 20

対策はかんたん
▌2つのlambdaに順序制約をつけてあげるだけ
model_mixture.stan
data {

model {

int N;

lambda ~ cauchy(0, 10);

array[N] int N_BUY;

pi ~ beta(1,1);

}

real logp0;
real logp1;

parameters {
ordered[2] lambda;

for(i in 1:N){

vector<lower=0, upper=1> [N] pi;
}

logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]);
logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]);

【ordered型】
基本的にはvector型だが

target += log_sum_exp(logp1, logp0);

中の成分が必ず小さい順になる
という制約が追加される
▶ lambda[1] < lambda[2]
【注意】lower, upperの制約を
同時にかけることはできない

}

}

13 混合分布モデル

21

改めて結果を見てみると 各クラスターの 平均購入数 各顧客が 「ヘビーユーザー」 に所属する確率 ▌クラスターがうまく分離できた！ 𝜆𝐻 っぽい lambdaのトレースプロットもみてみると bayesplot::mcmc_trace(result$draws("lambda")) 𝜆𝐿 っぽい 全てのchainで 𝜆ℎ , 𝜆𝐿 の割り当てが おなじになりました 13 混合分布モデル 22

事後予測チェックとかしたい
▌事後予測分布はどうなる？
データ発生メカニズムに従って考えると

𝑦𝑖 ∼ ቊ

𝑃𝑜𝑖𝑠(𝜆𝐿 ) if 𝑧𝑖 = 0
𝑃𝑜𝑖𝑠(𝜆𝐻 ) if 𝑧𝑖 = 1

𝑧𝑖 ∼ 𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝜋𝑖 )なので…

▶ まずzを乱数発生させて，その値に応じてpoissonのパラメータを変えたら良さそう！
data {
int N;
array[N] int N_BUY;
}
parameters {
ordered[2] lambda;
vector<lower=0, upper=1> [N] pi;
}
model {
lambda ~ cauchy(0, 10);
pi ~ beta(1,1);
real logp0;
real logp1;
// 右上に続く

}

// modelブロックの続き
for(i in 1:N){
logp1 = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]);
logp0 = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]);
target += log_sum_exp(logp1, logp0);
}

1

𝑧𝑖 ∼ 𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝜋𝑖 )
にしたがって乱数発生

generated quantities {
array[N] real N_BUY_pred;
for (i in 1:N){
N_BUY_pred[i] = poisson_rng(lambda[bernoulli_rng(pi[i]) ? 2 : 1]);
}
}

3 ②で決定したポアソン分布
から乱数を生成

13 混合分布モデル

2

三項演算子を用いて
• 乱数が1ならlambda[2]
• 乱数が0ならlambda[1]

23

再び推定して事後予測チェック ▌データの分布 # 推定が終わった後で N_BUY_pred <- result$draws("N_BUY_pred", format = "matrix")[seq(1,4000,by=500),] ppc_dens_overlay(y=dat$n_buy, y_pred= N_BUY_pred) このあたりの盛り上がりを うまく表現できていそう 13 混合分布モデル 【ちなみに】単一のポアソン分布で 推定→事後予測チェックした場合 24

ただ混合分布も完璧ではない ▌あくまでも「2つのポアソン分布」しか考えていない ▶ それ以上に細かな変化には対応しきれない # 推定が終わった後で ppc_scatter_avg(dat$n_buy, yrep=result$draws("N_BUY_pred",format="matrix")) ヘビーユーザーの分布 （𝜆𝐻 ）よりも大きな値 ライトユーザーの分布 （𝜆𝐿 ）よりも小さな値 13 混合分布モデル 25

混合分布モデルのその他もろもろ ▌もちろん，混ぜる確率分布はポアソン分布に限らない 特に，正規分布の混合分布モデルは「混合ガウスモデル (Gaussian mixture)」 ▌もちろん，混ぜる確率分布の個数は3個以上でも良い 𝐾個の分布の混合を考える場合も基本は同じ 𝐾 𝑃 𝑦𝑖 = ෍ 𝑃 𝑦𝑖 𝑧𝑖 𝑃(𝑧𝑖 = 𝑘) 𝑧𝑖 が1から𝐾になりうるので，categorical分布が使える 𝑘=1 事前分布はdirichletがおすすめ ▌混ざっている確率分布の個数がわからない場合は？ いくつかの個数で試してモデル比較する，といった方法が考えられます ▲ データから「最適な混合の個数」を予測させる 13 混合分布モデル 混合の個数をdataブロックで定義してRから与えてやることで 再コンパイルの必要なく様々な個数を試すことができます 26

２ 分布を分岐させる 実世界では案外良くあるかもしれない話 13 混合分布モデル 27

話のつづき… ▌データ発生メカニズムをもっと考えてみると… そもそもグミを全く買わない層っていうのがあるのでは？ グミを 全く買わない人 𝑃 𝑦𝑖 = 0 = 1 グミを買う人 𝑦𝑖 ∼ 𝑃𝑜𝑖𝑠(𝜆) 全ての顧客 問題を簡単にするために 「グミを買う人は全員同じクラスター」だと仮定しますが， もちろん先程の「ライト／ヘビーユーザー」という分類を 組み合わせることも可能です。 1 グミを 全く買わない人 全ての顧客 13 混合分布モデル =0 =1 1 ライトユーザー ( ) ヘビーユーザー ( グミを買う人 28 )

とはいえ基本的な考え方は同じです ▌データの値によって混合されるかどうかが変わる 1 − 𝜋𝑖 グミを 全く買わない人 𝑃 𝑦𝑖 = 0 = 1 𝜋𝑖 グミを買う人 𝑦𝑖 ∼ 𝑃𝑜𝑖𝑠(𝜆) 全ての顧客 先程の混合分布モデルとの違いは • 混合の有無がデータによって異なる • 一方の分布が確率分布ではなく定数 の2点です。 𝑦𝑖 ≠ 0の場合 𝑦𝑖 = 0の場合 購入個数0というデータには • もともとグミを買わない人 • グミを買うこともあるがたまたま今週は買わなかった人 が混ざっている，と考えている 13 混合分布モデル 29

コードを拡張していきます
▌p. 21のモデルをもとに
model_zip.stan
data {

model {

int N;

lambda ~ cauchy(0, 10);

array[N] int N_BUY;

pi ~ beta(1,1);

}
for(i in 1:N){

parameters {

// この中をどう書けば良い？

real<lower=0> lambda;
vector<lower=0, upper=1> [N] pi;

}

}

}

13 混合分布モデル

30

少しずつ丁寧に書いていくしか無いのですが ▌N_BUY[i] > 0の場合 グミを ▶ 全く買わない人 はいないので， ポアソン分布だけ考えたらよい の グミを買う人 // N_BUY[i] > 0の場合 1 target += poisson_lpmf(N_BUY[i]|lambda); ▶ と 混合分布を考える必要がある グミを買う人 // N_BUY[i] == 0の場合 logp_not_buy = log(1-pi[i]); =0 =1 グミを買う人 ( ) 全ての顧客 ▌N_BUY[i] == 0の場合 グミを 全く買わない人 グミを 全く買わない人 の2つの 0の場合 グミを全く買わないクラスターでは 𝑃 𝑦𝑖 = 0 = 1 なので，log 1 = 0は無視できる = 0の場合 logp_buy = poisson_lpmf(N_BUY[i]|lambda)+log(pi[i]); target += log_sum_exp(logp_not_buy, logp_buy); 13 混合分布モデル 31

コードを拡張していきます
▌p. 21のモデルをもとに
model_zip.stan
data {

// modelブロックの続き

int N;

1

for(i in 1:N){

array[N] int N_BUY;

if(N_BUY[i] == 0){

}

グミを
全く買わない人

=0 =1

グミを買う人

( )

全ての顧客

// N_BUY[i] == 0の場合
parameters {

logp_not_buy = log(1-pi[i]);

real<lower=0> lambda;

= 0の場合

logp_buy = poisson_lpmf(N_BUY[i]|lambda)+log(pi[i]);

vector<lower=0, upper=1> [N] pi;

target += log_sum_exp(logp_not_buy, logp_buy);

}

} else {
// N_BUY[i] > 0の場合

model {

target += poisson_lpmf(N_BUY[i]|lambda);

lambda ~ cauchy(0, 10);
}

pi ~ beta(1,1);
real logp_not_buy;
real logp_buy;
// 右上に続く

}

グミを買う人

}

( )
0の場合

13 混合分布モデル

32

事後予測チェックも同じように ▌事後予測分布はどうなる？ データ発生メカニズムに従って考えると 【𝑦𝑖 = 0の場合】 =0 if 𝑧𝑖 = 0 𝑦𝑖 ቊ ∼ 𝑃𝑜𝑖𝑠(𝜆) if 𝑧𝑖 = 1 【𝑦𝑖 > 0の場合】 𝑦𝑖 ∼ 𝑃𝑜𝑖𝑠(𝜆) 𝑧𝑖 ∼ 𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝜋𝑖 )なので… ▶ N_BUY[i]が0のときだけzを乱数発生させて，と処理したら良さそう！ model_zip.stanに追記 // modelブロックまでは前ページと同じなので省略 generated quantities { array[N] real N_BUY_pred; for (i in 1:N){ if(N_BUY[i] == 0) { 1 N_BUY[i] == 0なら まず𝑧𝑖 ∼ 𝑏𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝜋𝑖 ) にしたがって乱数発生 N_BUY_pred[i] = bernoulli_rng(pi[i]) ? poisson_rng(lambda) : 0; } else { N_BUY_pred[i] = poisson_rng(lambda); } } } 3 2 N_BUY[i]>0なら ポアソン分布から乱数を生成 13 混合分布モデル 三項演算子を用いて • 乱数が1ならポアソン分布から乱数を生成 • 乱数が0なら0を返す 33

推定して事後予測チェック ▌データの分布 result <- model$sample(data=stan_data, parallel_chains = 4) N_BUY_pred <- result$draws("N_BUY_pred", format = "matrix")[seq(1,4000,by=500),] ppc_dens_overlay(y=dat$n_buy, y_pred= N_BUY_pred) 単一のポアソン分布よりは 「0が多い」ことは 今回のデータに関しては 表現できているが… どうやら違ったようです… 13 混合分布モデル 【ちなみに】単一のポアソン分布で 推定→事後予測チェックした場合 34

今回紹介したモデル ▌ゼロ過剰ポアソン分布 (zero-inflated Poisson: ZIP) と呼ばれます 他の分布でも「ゼロ過剰」を考えることはできるはず （e.g., ゼロ過剰二項分布，ゼロ過剰負の二項分布，0-1過剰ベータ分布） ▌（混合分布も同じですが）混合比率に説明変数をいれることも可能 例えば「どういう属性にヘビーユーザーが多いか」などの検討が可能に ▲ ロジスティック回帰的なイメージで，ロジットを説明変数の線形結合で表現するなど ▌よく似たモデルにハードルモデル (Hurdle model) というのもあります 違いは𝑦𝑖 = 0のときに混合する分布 【ゼロ過剰モデルの場合】 =0 if 𝑧𝑖 = 0 𝑦𝑖 ቊ ∼ 𝑃𝑜𝑖𝑠(𝜆) if 𝑧𝑖 = 1 【ハードルモデルの場合】 =0 if 𝑧𝑖 = 0 𝑦𝑖 ቊ ∼ 𝑃𝑜𝑖𝑠𝑦>0 (𝜆) if 𝑧𝑖 = 1 13 混合分布モデル 𝑃𝑜𝑖𝑠𝑦>0 (𝜆)は切断ポアソン分布 ▶ 違いはポアソン分布に0を含むかどうかだけ 35

（おまけ）ゼロ過剰と混合分布を混ぜたコード例
data {
int N;
array[N] int N_BUY;
}
parameters {
ordered[2] lambda;
vector<lower=0, upper=1> [N] pi;
vector<lower=0, upper=1> [N] tau;
}
model {
lambda ~ cauchy(0, 10);
// 「ヘビーユーザー」である確率
pi ~ beta(1,1);
// 「買う人」である確率
tau ~ beta(1,1);
real logp_not_buy;
real logp_buy_L;
real logp_buy_H;
real logp_buy;
// 右上に続く

たぶん合っていますが
間違いに気づいたらこっそり教えて下さい…

model_zip_mixture.stan

// modelブロックの続き
for(i in 1:N){
if(N_BUY[i] == 0){
// N_BUY == 0の場合
logp_not_buy = log(1-tau[i]);
logp_buy_H = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]);
logp_buy_L = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]);
// まずライト・ヘビーの混合
// そこに「買う」の確率tau[i]の重みをつけておく
logp_buy = log_sum_exp(logp_buy_H, logp_buy_L) + log(tau[i]);
// 買わない・買うの混合
target += log_sum_exp(logp_not_buy, logp_buy);
} else {
// N_BUY > 0の場合
logp_buy_H = poisson_lpmf(N_BUY[i]|lambda[2])+log(pi[i]);
logp_buy_L = poisson_lpmf(N_BUY[i]|lambda[1])+log(1-pi[i]);
target += log_sum_exp(logp_buy_H, logp_buy_L);
}
}

}
generated quantities {
array[N] real N_BUY_pred;
for (i in 1:N){
if(N_BUY[i] == 0) {
N_BUY_pred[i] = bernoulli_rng(tau[i]) ? poisson_rng(lambda[bernoulli_rng(pi[i]) ? 2 :
1]) : 0;
} else {
N_BUY_pred[i] = poisson_rng(lambda[bernoulli_rng(pi[i]) ? 2 : 1]);
}
}
}

13 混合分布モデル

36

まとめと次回予告 【まとめ】 ▌確率分布を混ぜるモデルを紹介しました 実世界の複雑なデータに対しては，単一の確率分布では足りないことは多々ある 離散パラメータをstanで扱う場合には周辺化が必要 【次回予告】 ▌受講生の方に研究事例を紹介してもらいます 実際の研究でベイズ統計がどのように使われているのか ここまでついてこれた皆さんならもう理解できるだけの力はついているはずです 13 混合分布モデル 37