【例題でわかりやすく!】ベイズ推定でポアソン分布が更新される過程をエクセルで可視化

2021年8月31日

物流センターXでは返品を受け付けることになりました。

荷主の過去の経験から、毎日平均3±2パレット分くらいの物量がありそうです。

実際に業務を開始したら、下記のような返品がありました。

返品パレット数
1 3
2 3
3 0
4 1
5 4
6 3
7 4
8 0
9 1
10 1
11 2
12 3
13 1
14 4
15 1

 

荷主の過去の経験を事前確率分布として、7日後、15日後の事後確率分布を作ってみましょう。

 

事前確率分布を作る

返品数はポアソン分布に従うと仮定します。

ベイズ推定

事後確率=事前確率×尤度

で表現されますが、尤度はポアソン分布の式になるということです。

 

ポアソン分布の式はこうでした。

f(x) = λx e–λ / x!

λ:ある期間内に起きる平均回数

x:ある期間内に起きる知りたい回数

>> 【具体例でわかりやすく!】ポアソン分布は何に使えるの?|死に筋商品への応用を実演

 

この式を事前分布の式と掛け合わせるのです。

計算が難しそうですね。

でも、この式は指数関数です。

指数関数同士の掛け算なら

ea ×eb=ea+b

というように、指数同士の足し算にできるので簡単です。

従って、事前分布も同じ指数分布の形にしてやれば、計算が簡単になります。

このように尤度関数との計算を簡単にできる事前分布を共役事前分布といいます。

 

ポアソン分布の共役事前分布はガンマ関数です。

ガンマ関数は次式で表せます。

f(x) = λα xα-1eλx / Γ(α)

x:サイクル数

λ:サイクル内で事象が起こる平均回数

α:事象が起こる回数(確率を知りたい回数)

Γ(α):確率密度関数であるf(x)を積分して1になるようにするための係数

>> 【具体例でわかりやすく!】ガンマ分布は何に使えるの?|ヒヤリハットへの応用を実演

 

そして、ガンマ分布の平均(期待値)と分散

平均=α/λ

分散=α/λ2

になることが知られています。

 

これを使って荷主の経験を式で表現してみましょう。

荷主は平均3と言っていましたので、

α/λ=3

ですね。

それと、プラスマイナス2くらい振れると言っていましたので、標準偏差は2とします。

標準偏差の二乗が分散ですので、

α/λ2=4

です。

これを解くと、

α=2.25

λ=0.75

になります。

 

従って事前分布は、

f(x) = 0.752.25 x(2.25-1)e-0.75x / Γ(2.25)

となりますが、0.752.25とΓ(2.25)は定数なので、まとめてAとすると、

f(x) = A x1.25e-0.75x

と簡単にできます。

これが事前確率分布になります。

 

グラフで描くと次のようになります。

 

ここで注意です。

平均は3のはずなのに、このグラフの頂点は1.6くらいになっていますね。

これはガンマ分布に対称性がないため仕方ないのです。

グラフが左に寄っているため平均がどこか分かりにくいのですが、ちゃんと平均は3標準偏差は2になっています。

 

ベイズ推定する

それではこの事前確率分布を15個のデータでベイズ推定していきましょう。

1日目には3パレット返品されてきました。

ポアソン分布の式にそのまま当てはめると、

f(3) = λ3 e–λ / 3!

となりますね。

しかし、これはxを変数、λを定数と見た場合です。

λはポアソン分布の形を特徴づけるパラメータですが、今回はこれが未知です。

 

逆にxは15個のデータが得られています。

ですので、今回はλを変数、xを定数とします。

すると先ほどの式は

f(λ) = λ3 e–λ / 3!

のようにλの関数として書くことができます。

これが1日目の尤度になります。

 

事前分布は

f(x) = A x1.25e-0.75x

でした。

これは平均3、標準偏差2の確率分布の式です。

今回知りたいのはλなので、確率変数をλに変えます。

f(λ) = A λ1.25e-0.75λ

これが事前分布になります。

 

よって、1日目の事後分布は次のようになります。

事前分布×尤度

= A λ1.25e-0.75λ λ3 e–λ / 3!

= A1 λ(1.25+3)e-(0.75+1)λ

A/3!を新たな定数A1にしています。

 

同様に、2日目の事後分布は次のようになります。

A1 λ(1.25+3)e-(0.75+1)λ λ3 e–λ / 3!

= A2 λ(1.25+3+3)e-(0.75+1+1)λ

A1/3!を新たな定数A2にしています。

 

同様に計算すると、7日目の事後分布は次のようになります。

A λ(1.25+3+3+0+1+4+3+4)e-(0.75+1+1+1+1+1+1+1)λ

= A7 λ19.25e-7.75λ

 

λの指数部分は1.25+7日間の返品数の合計、eの指数部分は0.75+データの数(=7になっていることが分かります。

 

グラフで描くと次のようになります。

 

事前分布のグラフと比べてみて下さい。

だいぶ山が急になって、λの推定値が絞れてきていることが分かります。

また7日間の平均返品数は2.6なのですが、山の中心が2.6に近づいたことも分かります。

 

では同様にして15日目の事後分布も求めてみましょう。

A15 λ32.25e-15.75λ

になりますね。

グラフにするとこうなります。

 

7日目より更に山が急になって、λの推定値の精度が上がっていることが分かりますね。

 

3つのグラフの描き方ですが、まずはエクセルシート上で下図のように計算します。

後は、これを折れ線グラフにするだけです。

 

最後に、事前分布、7日目の事後分布、15日目の事後分布をまとめて描くとこうなります。

 

新しいデータで更新されるたびに、確信度が上がっていく様子が見てとれます。