様々な確率分布とC++を使った乱数生成
<これはKogakuin Univ Advent Calendar 2020 10日目の記事です>
はじめに
どうも2年ぶりです、Linebrackと申します。
嫌な事件だったね…(去年のアドベントカレンダー)
まあ、気を取り直して今年も(?)書きます。
背景と目的
最近はPythonでベイズ推定について学んでいます。このような分野では確率分布という概念を積極的に扱いまして、乱数の生成とも密接な関係を持っています。中でも、乱数生成に興味を持っている人は多いのではないかと思います。そのため、今回は確率分布に従った乱数を生成しようと思います。具体的には、C/C++を使って、一様分布、正規分布、ベータ分布、ポアソン分布の疑似乱数を生成します。そして、pythonでヒストグラムを書いて実際の確率密度関数と比較して妥当か調べます。
理論
確率分布
まず、とは確率変数に対して、各々の値をとる確率を表したものである(Wikiより)。ここで、確率変数とは確率的に様々な値をとる変数であり、離散的なものと連続的なものがある。離散の例としては、サイコロの出目X∈{1,2,..,5,6}がわかりやすい。一方、連続なものは文字通り連続的な出力が得られるものである。もし、得られたデータから確率分布を見たい場合はヒストグラムを書いてみるといいだろう。また、代表的な確率分布には確率密度関数というものが決まっており、X=xのとき、確率p(X=x)を計算することができる。さて、代表的な確率密度関数を見てみよう。
一様分布
一様分布は、確率変数がどの値でも確率が同じ分布である。下図は連続な場合であるが、aからbの範囲の乱数が一様に得られることを示す。これは、C言語でいうrand()関数で生成できる。また、他の乱数生成アルゴリズムとして線形合同法、メルセンヌツイスタが利用できる。
正規分布
正規分布は、平均μ,標準偏差σからなる確率分布である。下図のような中心で確率が最も高く、離れるほど確率が低くなるような分布となっているが、これが実世界の多くの現象にあてはまる。例えば、身長が平均に近い人と身長が極端に高い人の比率を考えると、平均に近い人の方が多くいると考えられる。
ベータ分布
ベータ分布はベイズ推定によく使われる分布である。パラメータはaとbを持ち、それはa-1回当たり、b-1回外れるときの分布に対応している。具体的な説明は省くが、もっと詳しく知りたい方は
ベータ分布の意味と平均・分散の導出 | 高校数学の美しい物語
ポアソン分布
ポアソン分布は、ランダムなイベントが一定時間内に何回発生するかの分布である。パラメータはλで、イベントが単位時間内に平均λ回起きることに対応している。また、ポアソン分布の乱数は回数となるため離散であることに注意する。これも具体的な説明は省きますが、
ポアソン分布の意味と平均・分散 | 高校数学の美しい物語 参考にどうぞ。
乱数生成の実験
一様分布の乱数
ここでは、以下の方法を用いて0から1の連続な乱数をとることを考える。
線形合同法
CongruenceMethodクラスとして実装した。a,b,cの値を決めることで、合同式による漸化式の更新で乱数を生成することができる。
メルセンヌツイスタ
具体的な理論は良くわからないが、今回は線形合同法の比較としてメルセンヌツイスタを使った。
確率分布に従った乱数の生成
理論で紹介したどの分布も0から1の連続一様乱数からの変換で生成できる。具体的に、一様分布・正規分布・ポアソン分布は[1]のpdfを参考にした。 ベータ分布は[2]の基本アルゴリズムを用いた。
乱数の生成とpythonによる描画
C/C++で以上のアルゴリズムを用いて10000個づつ乱数を作り、csvファイルとして保存し、それをpythonで確率密度関数と共に描画した。
結果
線形合同法を用いた場合
一様分布
正規分布
ベータ分布
ポアソン分布
メルセンヌツイスタを用いた場合
一様分布
正規分布
ベータ分布
ポアソン分布
考察
どちらの方法でも、大体うまく疑似乱数が生成できた。しかしながら線形合同法を用いたベータ分布のサンプリングでは乱数がある周期でループしてしまい、疑似乱数の生成に失敗した。このことから、周期の長いメルセンヌツイスタの使用を推奨する。
おわりに
今回使ったコードはgithubで共有します！(https://github.com/Line-brack/KogakuinAdvent2020_12_10)
本当はc++で確率分布の乱数生成をやりたかっただけなのに、細かいアルゴリズムから実装したくなり、内容が盛りだくさんになってしまった...。どうかこの記事を見た方には、様々な分布があることだけでも知ってもらえば幸いです。
参考
[1]http://sys.ci.ritsumei.ac.jp/probability/2011/10-6.pdf
[2]
ベータ分布のサンプリングアルゴリズムを読みとく - Qiita
[3]C言語による乱数生成
