確率分布からの乱数の生成1 確率変数の変換

ということで,下記のエントリのGMappingの解剖実験進めているんですが,Motion Model の計算のところで正規分布に従う変数を生成する必要があったので簡単にまとめておくことにしました.

確率分布と分布に従う乱数の生成

平均値μ,分散σの正規分布が下記のように書けることは確率統計の本を書けば書いてあると思います.

f:id:rkoichi2001:20181216223448j:plain
正規分布確率密度関数と分布の様子

で上式なんですが,「考えている確率変数Xが平均μ,分散σに従うとき,取り出した値がxになる確率はいくらか?」ってことには答えてくれます.ですが,正規分布に従う乱数を発生させる」ことにはそのまま使えません.

ここがこのエントリのテーマなんですが,ロボットのオドメトリとかセンサの計測モデルとかを考えるうえで,適当な誤差をエイやっと載せて計測の不確実性をモデルに取り込みたいときが多々あるわけです.正規分布に従う乱数発生器は Python とかならライブラリとして提供されていますが,C/C++だとどこかから探してこないといけません.

パソコンで発生させられる乱数

で,C/C++だと乱数発生器は全くないのか?という話なんですが,一様乱数を発生させる関数としていくつか関数が標準ライブラリに用意されてます.下記,0から32ビット最大値までの乱数を生成する関数のサンプルプログラムです.

#include <iostream>
#include <random>

int main(int, char**) {

  std::cout << __FILE__ << std::endl;
  std::mt19937 mt((int)time(0));

  for (int i = 0; i < 10; i++) {
    std::cout << mt() << std::endl;
  }

  return 0;
}

Pythonを使って実験&プロット

C++だと可視化するのがめんどくさいので,ここからはPythonを使って実験します.上記の一様乱数をPythonを使って生成し,ヒストグラムを作ってみます.上のC++の関数は32ビットの最大値まで生成しましたが,ここでは 0 ~ 1 までの乱数を生成します.

import matplotlib.pyplot as plt

import numpy as np
from scipy.stats import norm
import random

def sample_normal_distribution(avg, var):
 
  vars = np.array([random.uniform(-var, var) for i in range(12)])
  return 0.5 * np.sum(vars)

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

  n = np.linspace(-10, 10, 201)

  random_nums = np.array([random.random() for i in range(1000)])
  print(random_nums)
 
  plt.xlim((0, 1))
  plt.ylim((0, 20))
  plt.hist(random_nums, bins=101)
  plt.show()

Pythonを使って一様乱数のヒストグラムを生成

実行結果は下記のようになります.0から1までの一様乱数を1000個生成して,ビンが100個のヒストグラムを作ったので,各ビンで10個づつくらい要素があればだいたい一様乱数と言えるかと.

f:id:rkoichi2001:20181216230845p:plain
0から1までの一様乱数を1000個生成し,ヒストグラムを作成

一様乱数をつかった他の乱数の生成

で,問題はここからで,生成したい乱数は正規分布に従う乱数」でした.どうやるかという話なんですが,一様乱数(=一様分布に従う乱数)から他の分布に従う乱数を生成します.ちょっとここから数学の復習ですが...

確率変数Xから変換した確率変数Yの分布を求める.

分布がわかっているある確率変数をXとします.この確率変数Xにある変換をかけてY=Φ(X)を計算するとします.そうすると,このYも確率変数とみなせますが,分布が既にわかっているXの分布を使ってYの分布を計算することはできるでしょうか?XとYの確率密度関数(=分布)をそれぞれ f(x), g(y) と表すとします.この時,,,

f:id:rkoichi2001:20181216234409p:plain
確率変数の変換方法

ということで,確率分布が既知の確率変数Xを変換してできた確率変数Yの分布を求める方法がわかりました.

変数変換の例(指数分布)

ひとつ変数変換の実際の例を考えてみたいと思います.y = Φ(x) とした変換を具体的に y = -log(x) とすると,y の分布 g(y) は先ほど求めた関係式から下記のように計算されます.

f:id:rkoichi2001:20181217000029j:plain
一様分布を指数分布に変換した例

Pythonを使って確認.

ということで,上記の変数変換の結果が本当にあっているかどうか見てみます.

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

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

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

  random_nums = np.array([random.random() for i in range(gen_num)])
  trans_vars = -np.log(random_nums)
  hist = np.histogram(trans_vars, bins=bin_num, range=draw_range)
  histelem = np.array([hist[0][i] / bin_width for i in range(bin_num)])

  exp_x = np.array([bin_width * i for i in range(bin_num)])
  exp_vars = np.exp(-exp_x)

  plt.plot(exp_x, histelem)
  plt.plot(exp_x, exp_vars * gen_num)
  plt.show()

上記コードで生成した結果が,下記のプロットになります.ヒストグラムの結果と確率密度関数 x 生成個数の値がいい感じで一致してます.

f:id:rkoichi2001:20181217003919p:plain
指数乱数の生成

このエントリでわかったことは,

確率変数Xの確率密度関数 f(x) と確率変数XからYへの変換式 Y = Φ(X) がわかっている時に,どうやって確率変数Yの確率密度関数を求めるか?

でした.最終的にやりたい・知りたいことは,下記の状況で

どうやって g(y) に従う乱数(=正規乱数)を生成するか?(=一様分布からの変換式 Y = Φ(X) を求めるか)でした.

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

2. f(x) に従う乱数を生成する方法をもっている.(PCによる一様乱数生成)

3.確率変数Yの確率密度関数 g(y) もわかっている.(正規分布

ということで,次のエントリでやっとこさ乱数の生成に入ります...

参考文献

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

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