【サンプルファイル付き】ポアソン分布の乱数をエクセルで生成する2通りの方法
ポアソン分布の乱数生成は難しい
適正在庫のシミュレーションでは普通、過去の需要実績データを使います。
しかし、これは過去の特定の需要パターンでシミュレーションしたにすぎません。
需要パターンが変わった場合にも適正在庫にコントロールできることまでは示せません。
これを示すためにモンテカルロシミュレーションを行います。
各商品アイテムの需要の平均と標準偏差だけを指定して、それに従う正規分布乱数を生成し、それを需要データとしてシミュレーションするのです。
乱数は何個でも生成できるので、何百回でも何千回でも繰り返すことができ、シミュレーション結果の信頼性が高まります。
正規分布乱数を生成する方法については、
エクセルを使って正規分布の【疑似乱数】を生成する方法を具体例で実演！
で紹介しました。
で乱数を生成することができました。
ところが、ポアソン分布についてはこのような逆関数がエクセルにありません。
また、逆関数を式で解析的に求めることも困難です。
少しトリッキーなやり方が必要になりますので、2通りの方法を紹介します。
ポアソン分布については、下記の記事で具体的な使い方を解説しています。
【ポアソン分布の使い方】在庫管理への適用方法を具体例で解説します。
【バックスクリーン３連発】ポアソン分布で伝説が生まれた確率を計算してみた
ポアソン分布の累積分布関数からIf文で逆関数を求める方法
累積分布関数はエクセル関数にある
エクセルにはポアソン分布の累積分布確率を求める関数はあります。
=POISSON.DIST(起こる回数、λ、TRUE)
λ：単位期間内に起こる平均回数
λ=2の場合は次のように計算できます。
グラフにするとこうなります。
単位期間を1日とすると、λ=2ということは1日に平均2回起こる事象ということです。
例えば、1日に平均2個売れる商品を思い浮かべて下さい。
この商品が1日に1個も売れない確率は0.14、1個以下だけ売れる確率は0.4、2個以下だけ売れる確率は0.68、、、ということです。
累積分布関数の逆関数が解析的に計算できない
ここで縦軸は累積確率で0から1の値です。
ですので、0から1の一様乱数を生成してこれを縦軸の値とし、対応する横軸の値が分かればポアソン分布の乱数になります。
このように縦軸から横軸を求める関数は逆関数といいますが、ポアソン分布の場合、この逆関数を解析的に求めることが難しいのです。
If文で強制的に求める
でも幸いなことにポアソン分布は離散確率分布です。
離散ということは、1、2、3、、、というようにx軸は整数の値しか取らないのです。
従って力業が使えます。
次の図を見て下さい。
0から1の乱数を生成して、その値が
0から0.14の時・・・0
0.14から0.4の時・・・1
0.4から0.68の時・・・2
0.68から0.86の時・・・3
というように新しい数字を作れば、それはポアソン分布に従います。
従って、先ほど作ったシートを使って、次のようにポアソン乱数を生成できます。
できた1,000個のポアソン乱数をグラフにしてみると次のようになり、確かにλ=2のポアソン分布に従います。
指数分布の乱数からポアソン分布に変換する方法
ポアソン分布と指数分布は表裏の関係
もう一つのやり方は、ポアソン分布が指数分布と表裏の関係にあることを利用する方法です。
指数分布の逆関数は解析的に求められますので、乱数も簡単に生成できます。
できた指数分布の乱数を姉妹関係にあるポアソン分布の乱数に変換するのです。
指数分布とは次に起こるまでの時間の分布
【指数分布】具体例でわかりやすく解説します。｜ロングテール商品への応用
で解説したように、サイクル期間内で一定回数起こる事象があった時、それが最初に起こるまでのサイクル数の確率分布を指数分布といいます。
f(x) = λe–λx
x：サイクル数
λ：サイクル内で事象が起こる平均回数
もっと噛み砕いていうと、次に起こるまでの時間の分布です。
先ほどのポアソン分布と同じように、1日に2個しか売れない商品を思い浮かべて下さい。
今、商品が売れたとしたら、次に売れるまでの時間はどれくらいでしょうか？
エクセルで計算すると次のようになります。
指数分布の式で計算されるのは確率密度ですので、確率と累積確率も計算しています。
これによると、次に売れるまでの時間が0.5日、つまり12時間後になる確率は6.7%です。
随分低い確率のように見えますが、これは0.4～0.5日の間にだけ売れる確率なので、滅多に起きないことなのです。
累積確率で見ると、12時間以内に売れる確率は70%であることが分かります。
ポアソン分布とは単位期間内に何回起こるかの分布
一方で、ポアソン分布を噛み砕いていうと、単位期間内に何回起こるかを表しているので、指数関数と同じことを違う角度から見ていることになります。
どういうことかを具体的に説明します。
例えば今が日付が変わった午前0時で、1個目が売れたのが8時で、その6時間後に2個目、更にその12時間後に3個目が売れたとしましょう。
すると、1/3日（8時間）、1/4日（6時間）、1/2日（12時間）、、が指数分布に従います。
一方、2個、、、がポアソン分布に従います。
なぜポアソン分布の最初のデータが2かというと、3個目が売れたのは次の日ですので、1日目に売れたのは2個だからです。
このようにポアソン分布は指数分布を一日に売れた個数でみているだけですので、指数分布をポアソン分布に変換できます。
指数分布の数字の列を順に見ていって、足して1以内になる最大のデータ数がポアソン分布になります。
指数分布の逆関数は簡単に求められる
更に、指数分布は次式のように累積分布関数の逆関数が求められます。
x = -Log(確率)/λ
ですから「確率」に0から1の一様乱数を入れれば、xは指数分布に従います。
今回はλ=2ですので、
x = -Log(0～1の乱数)/2
で指数分布の乱数が次のように生成できます。
これでD列に指数分布の乱数が生成されました。
指数分布をポアソン分布に変換する
これをポアソン分布に変換します。
先ほど説明したように、指数分布の数字の列を順に見ていって、足して1以内になる最大のデータ数を調べます。
次のようにします。
指数分布の乱数を1,000個生成したとすると、これを手作業でやるのは大変です。
そこでVBA（マクロ）の力を借ります。
詳細は割愛しますが、次のようにプログラムすれば、E列にポアソン分布に従う乱数が自動的に入力されます。
Sub Macro1()
'
Dim SS As Double : Dim CNT As Long : Dim I As Long : Dim J As Long
'
SS = 0
CNT = 0
For I = 4 To 1004
If CNT <= 1 Then
If Range(Cells(I, 4), Cells(I, 4)) >= 1 Then
Range(Cells(I, 5), Cells(I, 5)) = 0
CNT = 0
Else
For J = 1 To 10
If J = 1 Then
SS = Range(Cells(I, 4), Cells(I, 4))
CNT = 0
End If
SS = SS + Range(Cells(I + J, 4), Cells(I + J, 4))
CNT = CNT + 1
If SS >= 1 Then
Range(Cells(I, 5), Cells(I, 5)) = CNT
SS = 0
Exit For
End If
Next J
End If
Else
CNT = CNT – 1
End If
Next I
'
End Sub
E列に作られた乱数をヒストグラムにすると、このようにλ=2のポアソン分布に従っていることが分かります。
このマクロが入ったサンプルファイルはこちらからダウンロードできますので、自由にご利用下さい。