逆関数法と棄却法、プログラム練習

統計的データ解析２０１４
2015.1.13
林田清
（大阪大学大学院理学研究科）
乱数発生関数の使用例
#include <stdlib.h>
#include <stdio.h>
注: RAND_MAXは
/usr/include/stdlib.hで整数の最大値として定義されている
/* The largest number rand will return (same as INT_MAX). */
#define RAND_MAX
2147483647
x=(double) (random()/(RAND_MAX+1));
とすると
hayasida@oskmule 347 > cc -o testran3 testran3.c
testran3.c: In function ‘main’:
testran3.c:14: 警告: 式の整数がオーバーフローしました
という警告がでて、xは全て0になってしまう。
int main()
{
int seed,i;
double x;
scanf("%d",&seed);
srandom(seed);
for(i=0;i<10;i++) {
x=(double) (random()/(RAND_MAX+1.0));
printf("%d %lf¥n",i, x);
}
}
一様乱数の応用例（簡単なモンテカ
ルロ計算）


一様擬似乱数を与える
円の面積(円周率の計算）




（擬似）一様乱数の組(x,y)をN個生成
半径１の円の中（x2+y2<1）に含まれる点の個数ｍ
(m/N) x4が pの近似値になるはず
面積の計算、従って積分にも応用できる（モンテ
カルロ積分）。
プログラム例

ここではseedは任意の数（で
きれば素数）を入力する。



#include <stdio.h>
#include <stdlib.h>
int main()
{
int seed, n, i, m;
double x,y,pi1,pi2,pi3;
srandom(time(NULL));と
して時間情報を種とするこ
ともできる。この場合は
time.hをincludeする必要
あり。
printf("input seed number¥n");
scanf("%d",&seed);
printf("input number of points ¥n");
scanf("%d",&n);
右の例でpi1は整数の割り算
が最初に行われるので値が
ゼロになる。 pi2, pi3のよう
な記法にせよ。
RAND_MAXは
/usr/include/stdlib.hの中で
定義(#define)されている。
srandom(seed);
m=0;
for(i=0;i<n;i++) {
x=((double)random())/(RAND_MAX+1.0);
y=((double)random())/(RAND_MAX+1.0);
if((x*x+y*y)<1.0) m++;
printf("%d x=%lf y=%lf¥n",i,x,y);
}
pi1 = m/n * 4.0;
pi2 = (m*4.0)/n;
pi3 = ((double) m)/((double) n )*4.0;
printf("pi=%lf %lf %lf¥n",pi1,pi2,pi3);
}
プリプロセッサ

#include <fname.h>


コンパイルされる際に、ソースコードのこの位置に
fname.h(通常は/usr/includeの下にある；具体的には
stdio.h,stdlib.h,math.h等）の内容が挿入される。
#define ABC abc


ソースコードのこの位置以降にあるABCという文字列
があれば、コンパイル時に自動的にabcにおきかえら
れて解釈される。
RAND_MAXは stdlib.hの中でdefineされている。
与えられた分布関数に従う乱数１:逆変換法
分布関数f ( x)に従う乱数xを発生させたい
累積分布関数F ( x)
F ( x)  
x
xmin
f ( x ')dx '
は区間(0,1)に一様にランダムに分布する。
もしFの逆関数F 1をみつけることができればx  F 1 (r )により
擬似一様乱数rから分布関数f ( x)に従う乱数xを求めることができる
例：f ( x)   exp( x)
F ( x)   f ( x ')dx '    exp( x ') 0  1  exp( x)
x
x
0
x
1

ln(1  F ( x))  
1

ln(r )
逆変換法の考え方
1
F ( x)  
一様乱数
xmin
（入力）
rが0<r<1で定義される一様乱数のとき
xの確率分布px ( x)はpx ( x)dx  pr (r )drから
算出できる。
dr
px ( x) 
( x(r  0)  x  x(r  1))
dx
f ( x ')dx '
r
確率分布は pｒ(r )dr  dr (0  r  1)
このrがある関数x(r )でxに変換されるとき
s
f ( x)
0
x
xmin
xmax
x  F -1 (r )と変換された乱数x
（出力）を得る
F -1 (r )が容易に計算できる
逆変換法ではr  F ( x)従ってx  F -1 (r )と
とる。すると
dF ( x)
px ( x) 
 f ( x) ( xmin  x  xmax )
dx
xの確率分布はf ( x)になる。
ことが必要
与えられた分布関数に従う乱数２:棄却法
発生が容易な分布関数（例えば一様分布）g ( x)の定数倍Cg ( x)に
よってf ( x)を完全に包み込む。
1.g ( x)の確率分布に従う乱数x0を発生させる
２．点x0でf ( x0 )とCg ( x0 )を計算する。　(0,1)の乱数rを発生させ、
(a)f ( x0 )  rCg ( x0 )ならば、x0を棄ててステップ１からやり直す
(b)f ( x0 )  rCg ( x0 )ならば、x0が求める乱数
P
Cg ( x0 )
Cg ( x)
x0を棄却
f ( x0 )
x0を採択
0
xmin
x0
f ( x)
x
xmax
正規乱数の発生

積分は誤差関数、解析的に記述できない。x   r  6
12

Box-Muller法か中心極限定理
i 1
ここでriは(0,1)の一様乱数
一般に確率密度関数q( x1 , x2 )に従う変数x1 , x2の関数y1 , y2の確率密度関数は
p( y1 , y2 )dy1 dy2  q( x1 , x2 )
 ( x1 , x2 )
dy1 dy2
 ( y1 , y2 )

 1
 1
exp( y2 2 / 2)  dy1 dy2
exp( y12 / 2)  
p( y1 , y2 )dy1 dy2  

  2p
 2p
となるような変数y1 , y2を生成したい  Box-Muller法
y1  2 ln x1 cos 2p x2 , y2  2 ln x1 sin 2p x2
x1 , x2が区間(0,1)の独立な一様乱数であればq( x1 , x2 )dx1 dx2  dx1 dx2
y
1

 1
arctan 2
x1  exp   ( y12  y2 2 )  , x2 
y1
2p

 2
 ( x1 , x2 )

 1
 1
exp( y2 2 / 2) 
exp( y12 / 2)  
 
 ( y1 , y2 )

  2p
 2p
i
ポアッソン分布に従う乱数
P( x)   x e-  / x !
ただしxは離散値、棄却法はそのままでは使えない
Q( y )dy   [ y ] e-  /[ y ]! という連続分布関数を定義する
ここで[y]はyを超えない整数
Q( y )に従う乱数yを棄却法により発生させ、その整数部を
とってxとする

Download Report