【サンプルファイル付き】ポアソン分布の乱数をエクセルで生成する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の一様乱数を生成してこれを縦軸の値とし、対応する横軸の値が分かればポアソン分布の乱数になります。

 

このように縦軸から横軸を求める関数は逆関数といいますが、ポアソン分布の場合、この逆関数を解析的に求めることが難しいのです。

でも幸いなことにポアソン分布は離散確率分布です。

離散ということは、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のポアソン分布に従っていることが分かります。

 

このマクロが入ったサンプルファイルはこちらからダウンロードできますので、自由にご利用下さい。