棄却サンプリング法とは？アルゴリズムをExcelに実装して試してみた。
ここでいうサンプリングとは、予め式が分かっている確率分布から、その分布に従う数字の羅列を得ることです。
マーケティング調査や音楽のサンプリングについて検索していた人は、これ以上読んでもムダですので退出して下さい。＾＾
サンプリングとは任意の確率分布に従う乱数を作ること
例えば、ある商品の毎日の売上数量が、平均100、標準偏差50の正規分布に従うものとしましょう。
この分布から10個のデータをサンプリングすると、
41、117、131、122、97、83、120、150、79、68
とか
10、85、76、148、118、83、134、142、82、54
といったような数字の羅列になります。
「それが何の役に立つのか？」
こういうサンプリングができると、例えば適正在庫シミュレーションなどで色々な売上パターンによるシミュレーションを行うことができるようになります。
「乱数と何が違うのか？」
同じことです。
ただ、
【サンプルファイル付き】ポアソン分布の乱数をエクセルで生成する2通りの方法
でポアソン分布に従う乱数の作り方を解説しましたが、このように有名な確率分布でなくても、サンプリング法を使えばどんな分布でも乱数を作ることができます。
本記事ではなかなかイメージの掴みにくいサンプリングについて、その入門編である棄却サンプリング法を使い慣れたExcelを使って解説します。
サンプリングの緩い説明
平均2、標準偏差0.5の正規分布を想定して、その分布からサンプリングしてみましょう。
横軸は確率変数zで、0から4までの目盛りでグラフを描いています。
このグラフから分かるように、2付近の数字が出る確率が高くて、0や4付近の数字が出る確率はほとんどゼロです。
ということは、この分布からサンプリングしたデータには2付近の数が多く含まれ、0や4付近の数はほとんど含まれないはずです。
サンプリングしたデータが正しいか間違っているかは、このようにして大体の当たりをつけることができます。
棄却サンプリングのアルゴリズム
それでは、このようなデータをサンプリングしてくるにはどうすれば良いでしょうか？
良い方法が知られています。
下の図を見て下さい。
まず0から4の数字を適当に選びます。
0から4までの間で一つの乱数を作るということです。
0.9が選ばれたとしましょう。
0.9に対応する確率密度は計算できます。
これは平均2、標準偏差0.5の正規分布のグラフなので、正規分布の式に代入して
確率密度p(z)
= 1/√(2π0.52)exp(-(0.9-2)2/0.52)
= 0.071
となります。
グラフの最高点は0.8くらいですが、ここを天井としましょう。
0.071というのは天井までの高さの約9%です。
これを0.9をサンプリングデータに入れる確率とすれば、都合が良さそうですね。
逆に言うと、91%の確率で0.9がサンプリングデータに入らないように棄却するということです。
次に2.4を考えてみましょう。
2.4に対応する確率密度は計算してみると0.58です。
このことから、2.4は58%の確率でサンプリングデータに入れる、逆に言うと42%の確率で棄却します。
このように生成された乱数を判断すれば、例えば2付近の数は、ほぼいつでもサンプリングデータに採用されますので、沢山の2付近の数がサンプリングデータに入ることになります。
逆に0付近や4付近の数は、ほとんどサンプリングデータに入りません。
願ったり叶ったりですね。
これをまとめて計算するには次のようにアルゴリズムを組めば良いでしょう。
- 0から4の乱数zを作る
- zに対応する確率密度p(z)を計算する
- p(z)／8を計算する
- 0から8の乱数qを作る
- q < p(z)／8ならzを採用、さもなければ棄却する
Excelに棄却サンプリング法を実装してみる
これをExcelでやってみましょう。
C列だけ正規分布の式を入れるので面倒ですが、他は簡単ですね。
F列にサンプリングされたデータが入ります。
1,000回、5,000回、10,000回、50,000回で試したところ、次のような結果になりました。
段々とサンプリングの精度が上がって、きれいな正規分布に近づいていることが分かります。
いずれの試行でも、採用率は30%くらいでした。
つまり、70%のデータはムダになって捨てられたわけです。
効率が悪いですね。
2付近のデータは高い確率で採用されますが、1以下や3以上のデータは高い確率で棄却されることが効いています。
今回は0から4から一様乱数を作りましたが、2付近が多く出る乱数を作れば採用率を上げれることが容易に想像できます。
そのために色々な改良版が開発されています。
重点サンプリングによるモンテカルロ積分をExcelで試して効果を見てみた。