ベイズ推定でポアソン分布のパラメータを更新する方法をExcelで実演！
予測返品数をポアソン分布でベイズ推定する事例
物流センター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になっています。
ベイズ推定で返品予測数を更新する
1日目のベイズ推定
それではこの事前確率分布を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にしています。
7日目のベイズ推定
同様に、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日目のベイズ推定
では同様にして15日目の事後分布も求めてみましょう。
A15 λ32.25e-15.75λ
になりますね。
グラフにするとこうなります。
7日目より更に山が急になって、λの推定値の精度が上がっていることが分かりますね。
更新回数が増えるにつれて確信度が上がる
3つのグラフの描き方ですが、まずはエクセルシート上で下図のように計算します。
後は、これを折れ線グラフにするだけです。
最後に、事前分布、7日目の事後分布、15日目の事後分布をまとめて描くとこうなります。
新しいデータで更新されるたびに、確信度が上がっていく様子が見てとれます。