ということで,下記のエントリのGMappingの解剖実験進めているんですが,Motion Model の計算のところで正規分布に従う変数を生成する必要があったので簡単にまとめておくことにしました.
確率分布と分布に従う乱数の生成
平均値μ,分散σの正規分布が下記のように書けることは確率統計の本を書けば書いてあると思います.
で上式なんですが,「考えている確率変数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()
一様乱数をつかった他の乱数の生成
で,問題はここからで,生成したい乱数は「正規分布に従う乱数」でした.どうやるかという話なんですが,一様乱数(=一様分布に従う乱数)から他の分布に従う乱数を生成します.ちょっとここから数学の復習ですが...
確率変数Xから変換した確率変数Yの分布を求める.
分布がわかっているある確率変数をXとします.この確率変数Xにある変換をかけてY=Φ(X)を計算するとします.そうすると,このYも確率変数とみなせますが,分布が既にわかっているXの分布を使ってYの分布を計算することはできるでしょうか?XとYの確率密度関数(=分布)をそれぞれ f(x), g(y) と表すとします.この時,,,
ということで,確率分布が既知の確率変数Xを変換してできた確率変数Yの分布を求める方法がわかりました.
変数変換の例(指数分布)
ひとつ変数変換の実際の例を考えてみたいと思います.y = Φ(x) とした変換を具体的に y = -log(x) とすると,y の分布 g(y) は先ほど求めた関係式から下記のように計算されます.
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 生成個数の値がいい感じで一致してます.
このエントリでわかったことは,
確率変数Xの確率密度関数 f(x) と確率変数XからYへの変換式 Y = Φ(X) がわかっている時に,どうやって確率変数Yの確率密度関数を求めるか?
でした.最終的にやりたい・知りたいことは,下記の状況で
どうやって g(y) に従う乱数(=正規乱数)を生成するか?(=一様分布からの変換式 Y = Φ(X) を求めるか)でした.
1.確率変数Xの確率密度関数 f(x) がわかっている.(一様分布)
2. f(x) に従う乱数を生成する方法をもっている.(PCによる一様乱数生成)
参考文献
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ
- 作者: William H. Press,William T. Vetterling,Saul A. Teukolsky,Brian P. Flannery,丹慶勝市,佐藤俊郎,奥村晴彦,小林誠
- 出版社/メーカー: 技術評論社
- 発売日: 1993/06/01
- メディア: 単行本
- 購入: 2人 クリック: 102回
- この商品を含むブログ (33件) を見る