確率分布からの乱数の生成2 分布に従う乱数の発生,Box Muller法

ということで前回のエントリでは,下記の条件が既知の場合に変換後変数の確率密度分布を求める方法を考えました.

1.変換前の確率変数Xの確率密度関数がわかっている.

2.確率変数Xと確率変数Yの変換式Y=Φ(X)

で,結果的に変換後の確率変数が従う分布 g(y) と変換前の確率変数が従う分布 f(x) の間には下記の関係があることがわかりました.

f:id:rkoichi2001:20181217020009j:plain
確率変数変換に伴う分布関数の変化

一様分布に従う確率変数からの変換式

ここで,PCでは一様乱数を生成するのは簡単であることを述べました.なので,変換前の確率変数 X を 0 から 1 の確率変数に絞って考えます.このことによって,上記の式は下記のようになります.

f:id:rkoichi2001:20181217022026j:plain
一様分布からの変数変換

いまの前提条件としては,Y = Φ(X) がわかっているという前提で話を進めてきました.なので,上記の変換式から変換後確率変数の分布を計算するのは簡単でした.(エントリ1の例では,y = Φ(x) を y = -log(x) として計算しましたが,x = exp(-y) とできるので,確率密度関数はすんなり求まりました.)

Y = Φ(X) の変換式(X = Φ-1(Y) )がわからないとき

で,実用上興味があるのはたいていこちらのパターンだと思います.つまり,前提条件として,,,,

1.変換前の確率変数Xの確率密度関数がわかっている(=一様分布)

2.変換後の確率変数Yの確率密度関数がわかっている(=乱数を生成したい分布)

で,この二つがわかった状態で,「変換後の確率密度関数に従う乱数を,一様乱数から生成したい」というのが今回の主目的でした.なので,要は Y = Φ(X) が知りたいわけです.で,もう一度下記の数式に戻りたいのですが,

f:id:rkoichi2001:20181217023037j:plain
一様分布からの変数変換

変換後の確率密度関数 p(y) を求めるには dx/dy の微分方程式を解く必要があります.が,ちょっと見方を変えてみると...p(y) は確率密度関数なので,その積分は累積分布関数 F(y) になります.なので,この微分方程式の解は x = F(y) と言えそうです.この時,累積分布関数 x = F(y) の逆関数 y = F-1(x) がわかれば,一様乱数から所望の分布に従う乱数を発生させる方法を手にしたことになります.この逆関数はすんなり求まるものもあれば,そうでないものもあるようです.

y = F-1(x) の幾何学的解釈

x が y の累積分布関数 F(y) という見地から見ると,一様乱数 x を入力として, y = F-1(x) を求めるということは,-∞からyまでの確率積分がxと等しくなるようなyを返しているということになります.ちょっとこんがらがってきたので,下記イラストです.

f:id:rkoichi2001:20181217024649j:plain
積分布関数の逆関数の解釈

ということで,一様分布から所望の分布に従う乱数を生成する方法でした....としたいところですが,具体的なところにまだ全く入ってないんですよね(笑)

Box Muller 法

上で,一様分布に従う確率変数 X から所望の分布に従う確率変数 Y を生成する方法を書きましたが,要は次の条件が満たされてないとだめでした.

1.所望の分布に従う確率変数Yを発生させるためには,Y=Φ(X) の変換式が必要.

2.所望の分布の累積分布関数(x = F(y))の逆関数(y = F-1(x))がわかれば,そいつが変換式になる.

で,今回乱数のケースでは,所望の分布が正規分布です.で,正規分布の累積分布関数が解析的に求まるかというと....求まりません(笑)なので,結局微分方程式にうまく"ハマる"解を見つけてあげなければならず,これを見つけたのが BoxさんとMulllerさんなんです.結果だけを言うと,下記の変換をすれば,y1 と y2 はそれぞれ平均0,分散1の正規分布に従う乱数になります.

f:id:rkoichi2001:20181217031510j:plain
Box Muller 法

上記の考察は,,,あきらめます.代わりに,Pythonで実験です(゚Д゚)ノ

Pythonを使って,Box-Muller 法の実験

Box-Muller法って,一様乱数が2つ必要な代わりに,生成する正規乱数も2つできるんですね.Numerical Recipes in Cには,三角関数を使わない高速バージョンが書かれてたんですが,夜も更けてきたので,一番簡単なものだけ実装しておきます.

import matplotlib.pyplot as plt
import numpy as np
import random
import math

if __name__ == '__main__':
  print(__file__)

  gen_num = 10001
  bin_num = 1001
  draw_range = (-5, 5)
  bin_width = float(((draw_range[1] - draw_range[0]) / float(bin_num)))

  x1 = np.array([random.random() for i in range(gen_num)])  
  x2 = np.array([random.random() for i in range(gen_num)])

  y1 = np.array([math.sqrt(-2*math.log(x1[i])) * math.cos(2 * math.pi * x2[i]) for i in range(gen_num)])
  y2 = np.array([math.sqrt(-2*math.log(x1[i])) * math.sin(2 * math.pi * x2[i]) for i in range(gen_num)])

  hist1 = np.histogram(y1, bins=bin_num, range=draw_range)
  histelem1 = np.array([hist1[0][i] / bin_width for i in range(bin_num)])

  hist2 = np.histogram(y2, bins=bin_num, range=draw_range)
  histelem2 = np.array([hist2[0][i] / bin_width for i in range(bin_num)])

  exp_x = np.array([bin_width * i for i in range(-bin_num/2, bin_num/2)])
  exp_vars = 1 / math.sqrt(2 * math.pi) * np.exp(-np.power(exp_x, 2) / 2.0)

  print(exp_x)
  #plt.plot(histelem)
  plt.plot(exp_x, histelem1)
  plt.plot(exp_x, histelem2)
  plt.plot(exp_x, exp_vars * gen_num)
  plt.show()

実行結果

f:id:rkoichi2001:20181217034615p:plain
Box-Muller法を使った正規乱数の生成

参考文献

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ