逆関数法による乱数生成

(1)

樋口さぶろお

龍谷大学理工学部数理情報学科

計算科学☆実習B L10(2016-06-27 Mon)

最終更新: Time-stamp: ”2016-06-27 Mon 18:36 JST hig”

今日の目標

逆関数法で,与えられた確率密度関数fR(r) ^にしたがう確率変数Rの擬似乱数が生成できる.

http://hig3.net

樋口さぶろお (数理情報学科) L10逆関数法による乱数生成計算科学☆実習B(2016) 1 / 23

(2)

連続型確率変数とその乱数

L09-Q1

Quiz解答:[a, b) 一様分布

1 C= _b₋¹_a.

2 E[X] = _b₋¹_a¹₂(b²−a^a) = ¹₂(a+b).

3 P(X ≥ a+ 2b

3 ) =

∫ _b

(a+2b)/3

1

b−adx= 1 3.

4 E[X²] = _b−a¹ ¹₃(b³−a³).

5 V[X] = E[X²]−(E[X])² = ₁₂¹(b−a)². L09-Q2

Quiz^解答:[a, b) ^{一様乱数の生成}

(3)

1 d o u b l e g e t r a n d o m (d o u b l e y ){

2 d o u b l e s ;

3 s =5.0∗y−3 . 0 ;

4 r e t u r n s ;

5 }

L09-Q3

Quiz解答:確率変数の変換 f_Q(q) = {

1 (0≤q <1)

0 (他) ^なので,

1

P(R < ¹₂a+b) =P(Q < ¹₂) =

∫ ¹

2 0

fQ(q) dq= ¹₂

(4)

連続型確率変数とその乱数 2 fR(r)dr =fQ(q)dq ^より,

f_R(r) = 1

dr dq

f_Q(q) = 1

a ×f_Q(q) = {1

a (b≤r < a+b) 0 (他)

q ^と r ^{の範囲に注意}.

q · · · 0 · · · 1 · · ·

r a a+b

なお,^この f_R(r) ^{を先に求めた場合は}, 1は次のように計算できる.

P(R < 1

2a+b) =

∫ ¹

2a+b

−∞ f_R(r) dr= 0 +

∫ ¹

2a+b b

1

a dr = 1 2. この f_R(r)は, L09-Q2 のdouble getrandom(double y)^{で生成される} 乱数の確率密度関数. ちゃんと一様になってるでしょ.

(5)

復習

:

正方形クッキーの質量から一辺の長さへの確率変数の変換

I

(6)

復習

:

正方形クッキーの質量から一辺の長さへの確率変数の変換

II

Quiz(確率変数の変換)

あるクッキー製造器が,ボールからとってくるひとかたまりの生地の質量はQg である. 確率変数Q は次の確率密度関数にしたがう.

f_Q(q) = {1

36 (64≤q <100) 0 (^他)

クッキー製造器では,Qg ^{の生地から},^一辺 R=g(Q) =√

Qcm ^の正方形のクッキーが焼ける(本当は比例定数あるけど省略).

1 Qの母平均値と母分散を求めよう.

2 確率 P(Q >82)^{を求めよう}.

3 クッキーの一辺 R cm のしたがう確率密度関数f_R(r) を求めよう.

4 R の母平均値と母分散を求めよう(2つの方法で).

5 確率 P(R >9)^{を求めよう}(2^{つの方法で}).

(7)

1 E[Q] =∫_+∞

−∞ q f_Q(q) dq = ⁶⁴⁺¹⁰⁰₂ V[Q] =· · ·= ₁₂¹(100−64)² = 108..

2 P(Q >82) = E[1_[q>82](Q)] =∫₁₀₀

82 1

36 dq = ¹₂.^これはQ^{が一様分布} だから.

3

fR(r) = 1

1 2

√1q

fQ(q) = 2rfQ(q(r)) = {r

18 (8≤r <10) 0 (^他)

1 y

1 2 s

1 y

1 2 s

(8)

(9)

母期待値

E[ϕ(Q)]=

母平均値

E[R]

の推定

確率変数の変換R=ϕ(Q) =g(Q)を考える.

座標の標本X⁽¹⁾(T),· · ·, X^(N)(T). 標本サイズ(例)N= 1000 標本からの推定

座標Qの母平均値E[Q]の推定値標本平均値 X= 1

N[X⁽¹⁾(T) +· · ·+X^(N)(T)]

.

座標の関数R=ϕ(Q)の母期待値E[R]の推定値標本期待値 ϕ(Q) =R= 1

N[R⁽¹⁾+· · ·+R^(N)], R⁽ⁱ⁾=ϕ(X⁽ⁱ⁾(T)) Rの母分散V[R]のの推定値

不偏標本分散 S_R² = 1

N−1[(R⁽ⁱ⁾−R)²+· · ·(R^(N)−R)²] = N

N−1[R²−(R)²] 母期待値E[ϕ(Q)]の区間推定は,母平均値E[R]の区間推定と同じこと.

確率空間(大学院) 樋口さぶろお (数理情報学科) L10逆関数法による乱数生成計算科学☆実習B(2016) 9 / 23

(10)

(

母分散未知の場合の

)

区間推定

R=ϕ(Q) の標本平均値としてm,不偏標本分散としてs² が得られたとき,^母期待値µ= E[R]^の,

信頼係数1−α ^{の信頼区間は} m−t_α/2(N−1)×√

s²/N < µ < m+t_α/2(N −1)×√ s²/N 異なるシードで

このポリシーで何回も推定した

とき,µ^{がこの不等式を満たす}(=^{信頼区間に含まれる})^確率は1−α.

標本サイズN が大きくなると信頼区間は

小さくなる

信頼係数1−α が大きくなると信頼区間は

大きくなる

計算科学☆実習B(2016)L03 確率統計☆演習I(2015)L11 確率統計☆演習I(2015)L12

(11)

ここまで来たよ

3 連続型確率変数とその乱数

4 逆関数法による乱数生成逆関数法

連続座標のランダムウォーク

(12)

逆関数法による乱数生成逆関数法

逆関数法

Q が[0,1)一様分布にしたがう(f_Qが区分的に定数)ときでも,R =g(Q) の確率密度関数 fR(r) ^{は定数ではない}. ^以下 Q=Y(^一様分布)

逆に,^{与えられた}fR(r) にしたがう擬似乱数を作りたい→ ^うまい g ^で R=g(Y)として作ろう!

変数変換 r=g(y)

[0,1)一様にしたがう y 0 → 1 f_R(r)にしたがう r r_min → r_max rmin, rmax は,f(r)>0 ^となるr ^の下限,^上限. r =g(y) ^{はどんな関数}?

(13)

原理fR(r) dr=fY(y) dy 微分方程式f_R(r)dr

dy =f_Y(y) 両辺を y^′= 0 ^から y ^まで積分

∫ _y

0

f_R(r)dr

dy(y^′)dy^′=

∫ _y

0

f_Y(y^′) dy^′

∫ _r=g(y)

r_min

fR(r^′)dr^′=y

P(R < g(y)) =P(Y < y) 左辺を次のように呼ぶ.

累積分布関数

F(r) =

∫ _r

rmin

f_R(r^′) dr^′(=

∫ _r

−∞f_R(r^′) dr^′)

(14)

累積分布関数

F(r)

の意味

F (r) = P (R < r)

r −∞ rmin rmax +∞

fR(r) =F^′(r) 0← 0 ≥0 0 →0 F(r)

0 ← 0 ↗ 1 → 1

0 1 2 3 4 5

0.0 0.2 0.4 0.6 0.8 1.0

(15)

動画 https://www.youtube.com/watch?v=cbgpdRW_5kQ

(16)

逆関数法の手続き

書き直すと

F(g(y)) =y

乱数を作りたい確率変数 Rの累積分布関数F(r) の逆関数がg(y).

逆関数法(逆変換法)

fR(r)^{に従う乱数を},r =g(y) ^で[0,1)^一様乱数Y ^{から作るには},g(y) ^を次の様に決めればいい.

1 R の累積分布関数 F(r) =

∫ _r

−∞f_R(r^′) dr^′ を計算する.

2 範囲 0≤y <1, rmin≤r < rmax で,y=F(r) ^を解いて,^逆関数 r=F⁻¹(y) =g(y) を求める.

1 d o u b l e g e t r a n d o m (d o u b l e y ){ /∗yは[ 0 , 1 )一様乱数∗/

2 r e t u r n g(y); /∗ 上で求めた式を書く ∗/

3 }

(17)

L10-Q1

Quiz(逆変換法) 確率密度関数

fR(r) = {1

2r (0≤r <2) 0 (^他)

に従う乱数 R を,[0,1)一様乱数 y から r=g(y) で作りたい. g(y) を求めよう.

(18)

L10-Q2

Quiz(逆関数法) 確率密度関数

fR(r) = {3√

2 8

√r (0≤r <2)

0 (他)

に従う乱数 R を,[0,1)一様乱数 y から r=g(y) で作りたい. g(y) を求めよう.

(19)

L10-Q3

Quiz([a, b) 一様乱数の生成)

次の確率密度関数を持つ確率変数Rを考える.

f(r) = {

1/5 (−3≤r <2) 0 (他)

Rに対応する疑似乱数を返す関数double getrandom(double y) ^を書こう. ^ただし,y ^としては[0,1)^{一様乱数を代入する}.

(20)

L10-Q4

Quiz(連続的な値をとる疑似乱数)

次の確率密度関数を持つ確率変数Rを考える.

f(r) =







4/3 (1/4≤r <1/2) 8/3 (1/2≤r <3/4) 0 (それ以外)

Rに対応する疑似乱数を返す関数double getrandom(double y) ^を書こう.

(21)

ここまで来たよ

3 連続型確率変数とその乱数

4 逆関数法による乱数生成逆関数法

連続座標のランダムウォーク

(22)

逆関数法による乱数生成連続座標のランダムウォーク

連続座標のランダムウォーク

t∈Z: 時刻時間離散

X(t+ 1) =X(t) +R(t+ 1)

これまで空間離散(離散座標),時間離散ランダムウォーク

R(t) 整数値をとる離散型確率変数 ⇝ X(t)整数値をとる離散型確率変数離散型確率変数 R(t) ^{は確率関数}P(R=k) =p_k ^{で指定される}.

これから空間連続(連続座標),時間離散ランダムウォーク R(t) ^{連続型確率変数} ⇝ X(t)^{連続型確率変数}

連続型確率変数 R(t) ^{は確率密度関数}fR(r) ^{で指定される}. オイラー表現によるマルコフ連鎖の数値計算はできない…

状態数が無限大だから

確率シミュレーションで.

(23)

お知らせ

任意レポートは 2016-06-29水まで→ 2016-07-06水まで月昼樋口オフィスアワー(1-502)

チューター/Math^{ラウンジ月火水木昼} 1-614

https://manaba.ryukoku.ac.jp

マイページの下の方に manaba^{出席カード} 提出