ベイズ推定でポアソン分布のパラメータを更新する方法を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日目の事後分布をまとめて描くとこうなります。
新しいデータで更新されるたびに、確信度が上がっていく様子が見てとれます。
【Udemyの関連講座】
Excel(エクセル)による確率分布基礎講座【確率思考を身につけよう】
統計学の確率分布という優れた道具を理解し、実際にエクセルを使って確率分布を実践することで、確率的に考える思考プロセス（確率思考）を身につけましょう！
はじめての統計（推定・検定編） ～記述統計から推測統計へ！しっかり9時間、97レクチャーでデータ時代の入場券を手に入れる
データサイエンス時代にまず押さえるべきデータの扱い方・見方を扱った統計講座。データをどう要約し、分かりやすく伝えるのか（記述統計）から、そのデータから母集団について何が言えるのか（推測統計）まで、丁寧に統計的発想を身に付けます
いちばん理解できる統計学ベーシック講座その１【確率分布・推定・検定】
統計学の基礎を効率的に学べるベーシック講座です。統計学の入り口となる「確率分布・推定・検定」について豊富な図を用いて説明していきます。統計学や確率思考という一生モノのスキルを一緒に身につけましょう！