区間 [0,1) の一様 (疑似) 乱数を生成する関数 double uniform(void)

計科 II-5.1

＜注＞

区間 [0,1) の一様 (疑似) 乱数を生成する関数 double uniform(void)

が用意されているとし，計科 II.5 の [問] で if( uniform() < Pr ){

S= 1;

} else if ( uniform() < Pr+Ps ){

S= 0;

} else { S= -1;

}

とすると，確率変数 S は

S =

 



 

1 ; 確率 P r

0 ; 確率 (1 − P r) × (P r + P s)

− 1 ; 確率 (1 − P r) × (1 − P r − P s)

(5.1.1)

に従って生成されることになり，意図と違う結果になることに注意。

計科 II.13 ( 参考 ) サンプルの不偏分散の式で 1/N − 1 をかける理由

不偏分散 S

²

の期待値を考えてみよう。

X

1

− X ¯

N

= (

1 − 1 N

)

X

1

− 1 N

∑

N

k=2

X

_k

(13.1)

なので

( X

₁

− X ¯

_N

)

2

= (

1 − 1 N

)

₂

(X

₁

)

²

− 2 (

1 − 1 N

) X

₁

1

N

∑

N

k=2

X

_k

+ 1 N

²

∑

N

k=2

(X

_k

)

²

+ 1 N

²

∑

k,l,k6=l

X

_k

X

_l

(13.2)

となる。この式の期待値 E( (

X

₁

− X ¯

_N

)

2

) は { X

_k

, k = 1, · · · , N } ^{がそれぞれ独立なので} (

1 − 1 N

)

₂

E(X

²

) − 2 (

1 − 1 N

) N − 1

N E(X)

²

+ N − 1

N

²

E(X

²

) + (N − 1)(N − 2)

N

²

E(X)

²

= N − 1 N

(

E(X

²

) − E(X)

²

)

= N − 1

N σ

²

(13.3)

となる。 S

²

の期待値は上式の N

N − 1 倍なので不偏分散の期待値が母分散に等しいことがわかる：

E(S

²

) = σ

²

(13.4)

(例) ランダムウォークの母平均と母分散の時間発展 一般の時刻 t でのランダムウォークの母平均

µ(t) =

∑

∞ x=−∞

xP (x, t) (13.5)

を計算してみよう。いきなり µ(t) を求めるのは難しそうなので，時間が 1 ステップ増えると µ(t) がどれだけ変化するかを考える。式 (8.1) より

µ(t + 1) =

∑

∞ x=−∞

xP (x, t + 1) =

∑

∞ x=−∞

x(pP (x − 1, t) + qP (x + 1, t)) (13.6)

となる。上式右辺第 1 項は p

∑

∞ x=−∞

xP (x − 1, t) = p

∑

∞ y=−∞

(y + 1)P (y, t) = p

∑

∞ y=−∞

yP (y, t) + p

∑

∞ y=−∞

P (y, t)

= p µ(t) + p (13.7)

となる。y = x − 1 であり，

∑

∞ y=−∞

P (y, t) = 1 を用いた。同様に

q

∑

∞ x=−∞

xP (x + 1, t) = q

∑

∞ y⁰=−∞

(y

⁰

− 1)P (y

⁰

, t) = q

∑

∞ y⁰=−∞

y

⁰

P (y

⁰

, t) − q

∑

∞ y⁰=−∞

P (y

⁰

, t)

= q µ(t) − q (13.8)

となるので

µ(t + 1) = µ(t) + p − q (13.9)

が得られる。µ(0) = 0 なので

µ(t) = (p − q) t (13.10)

となる。時間とともに p > q なら母平均は増加し，p < q なら減少する。

次に母分散を計算してみよう。まず，時刻 t + 1 でのランダムウォーカーの位置 X(t + 1) の 2 乗の期待値は

E (

X(t + 1)

²

)

=

∑

∞ x=−∞

x

²

P (x, t + 1) = p

∑

∞ x=−∞

x

²

P (x − 1, t) + q

∑

∞ x=−∞

x

²

P (x + 1, t)

= p

∑

∞ y=−∞

y

²

P (y, t) + 2p

∑

∞ x=−∞

y P (y, t) + p

∑

∞ x=−∞

P (y, t)

+ q

∑

∞ y⁰=−∞

y

⁰²

P (y

⁰

, t) − 2q

∑

∞ y⁰=−∞

y

⁰

P (y

⁰

, t) + q

∑

∞ y⁰=−∞

P (y

⁰

, t)

= (p + q)E(X(t)

²

) + 2(p − q)E(X(t)) + p + q

= E(X(t)

²

) + 2(p − q)µ(t) + 1 (14.11)

のように時刻 t の量と関係がつく。これより σ

²

(t + 1) = E

(

X(t + 1)

²

) − µ(t + 1)

²

= E(X(t)

²

) + 2(p − q)µ(t) + 1 − (

µ(t) + p − q )

2

= σ

²

(t) + 4pq (14.12) が得られ，母分散は t が 1 増えるごとに 4pq だけ増加することがわかる。。σ

²

(0) = 0 なので

σ

²

(t) = 4 p q t (14.13)

となる。

【問】

ランダムウォーカーが確率 p

R

で 右へ，確率 p

L

で左へ，確率 p

S

でその位置にとどまる場合の 時刻 t の母平均と母分散を求めなさい。

【答】

上と同様に時刻 t + 1 と t の量の関係をつけると

µ(t + 1) = µ(t) + p

_R

− p

_L

, σ

²

(t + 1) = σ

²

(t) + p

_R

+ p

_L

− (

p

_R

− p

_L

)

2

(14.14) が得られる。初期条件 µ(0) = 0 , σ

²

(0) = 0 より

µ(t) = (p

_R

− p

_L

) t , σ

²

(t) = {

p

_R

+ p

_L

− (

p

_R

− p

_L

)

2

}

t (14.15)

となる。尚，p

R

+ p

_L

+ p

_S

= 1 なので σ

²

(t) は次のようにも表せる。

σ

²

(t) = (

4p

R

p

L

+ p

S

(1 − p

S

) )

t (14.16)

計科 II-14.1

【注】生成関数

Z(s, t) =

∑

∞ x=−∞

s

^x

P (x, t) = (

ps + qs

⁻¹

)

t

(14.1.1) から µ(t) と σ

²

(t) を求める方法

µ(t) = ∑

x

xP (x, t) = dZ(s, t) ds

¯¯ ¯¯

s=1

= d ds

(

ps + qs

⁻¹

)

t

¯¯

¯¯

s=1

= [

t (

ps + qs

⁻¹

)

t−1

(

p − qs

⁻²

)]

s=1

= t(p + q)

^t−1

(p − q)

^p+q=1

= t(p − q) (14.1.2)

となり，確かに (13.10) と一致する。

E (

X(t)

²

)

= ∑

x

x

²

P (x, t) = [ d

ds s d

ds Z (s, t) ]

s=1

= t d ds

((

ps + qs

⁻¹

)

t−1

(

ps − qs

⁻¹

))¯¯

¯¯

s=1

= t [

(t − 1) (

ps + qs

⁻¹

)

t−2

(

p − qs

⁻²

)(

ps − qs

⁻¹

)

+ (

ps + qs

⁻¹

)

t−1

(

p + qs

⁻²

)]

s=1

= t [

(t − 1)(p + q)

^t−2

(p − q)

²

+ (p + q)

^t−1

(p + q) ]

p+q=1

= t (

1 + (t − 1)(p − q)

²

)

= t (

1 − (p − q)

²

)

+ t

²

(p − q)

²

(14.1.3)

σ

²

(t) = E (

X(t)

²

) − µ(t)

²

= E (

X(t)

²

) − t

²

(p − q)

²

= t (

1 − (p − q)

²

)

_(p+q)2=1

= 4pqt (14.1.4)

となり (14.12) と確かに一致する。

・確率変数のとる値が連続 (実数) の場合

確率変数 X が離散的な値 { x

₁

, x

₂

, · · · } をとる場合は X = x

_i

となる確率 P (x

_i

) を考えたが，

X のとる値が連続な場合は，X がある 1 点の値をとる確率，例えば (X = 1.0 となる確率) は 0 となる。(点の長さは 0 であるため。)

X のとる値がある区間 [a , b] に入る確率，P (a ≤ X ≤ b) を考える必要がある。

P (a ≤ X ≤ b) =

∫

_b

a

p(x) dx (15.1) である場合 p(x) ( ≥ 0) は 確率密度関数 と呼ば れる。 ∫

_∞

−∞

p(x) dx = 1 (15.2) となっている。また，X のとる値が x 以下である 確率 F (x) = P ( −∞ < X ≤ x) を x の関数と考え て 累積分布関数 と呼ぶ。

F (x) =

∫

x

−∞

p(x

⁰

) dx

⁰

(15.3)

-2 -1 0.5 1 2

0.1 0.2 0.3 0.4 0.5 0.6

p(x) = 1

√ π e

^−x2

＜注＞ dF (x)/dx = p(x) の関係がある。すなわち，F (x) は p(x) の原始関数の一つ。

・母平均 (離散的な場合の式 (10.1) に対応) µ = E(X) =

∫

_∞

−∞

x p(x) dx (15.4)

・関数 f(X) の期待値 (離散的な場合の式 (11.2) に対応) E(f (X)) =

∫

_∞

−∞

f(x) p(x) dx (15.5)

・母分散 (離散的な場合の式 (12.1) に対応) σ

²

= E

(

(X − µ)

²

)

=

∫

_∞

−∞

(x − µ)

²

p(x) dx =

∫

_∞

−∞

x

²

p(x) dx − µ

²

(15.6) (例) 区間 [0 , 1) の一様乱数

p(x) = {

1 ; 0 ≤ x < 1

0 ; x < 0 || 1 ≤ x (15.7)

double uniform(void){

double r;

r = ((double) rand())/((double)RAND_MAX + 1.0);

return(r);

}

他の確率分布に従う確率変数は uniform() から作る。

計科 II.16

・累積分布関数

F (x) =

 



 

0 ; x < 0 x ; 0 ≤ x < 1 1 ; 1 ≤ x

(16.1)

・母平均

µ =

∫

1 0

x dx = x

²

2

¯¯ ¯¯

^x=1

x=0

= 1

2 (16.2)

・母分散

σ

²

=

∫

₁

0

x

²

dx − µ

²

= x

³

3

¯¯ ¯¯

^x=1

x=0

− 1 4 = 1

12 (16.3)

(例) 区間 [a , b) の一様乱数

p(x) = {

1/(b − a) ; a ≤ x < b

0 ; x < a || b ≤ x (16.4)

double get_myrandom(void){

double x;

double r;

r=uniform();

x= r*(b -a) + a ; return x;

}

F(x)

x

a = 1，b = 3 の場合

・累積分布関数

F (x) =

 



 

0 ; x < a (x − a)/(b − a) ; a ≤ x < b

1 ; b ≤ x

(16.5)

・母平均

µ =

∫

b a

x

b − a dx = x

²

2(b − a)

¯¯ ¯¯

^x=b

x=a

= a + b

2 (16.6)

・母分散

σ

²

=

∫

b a

x

²

b − a dx − µ

²

= x

³

3(b − a)

¯¯ ¯¯

^x=b

x=a

− (a + b)

²

4 = (b − a)

²

12 (16.7)

(例)

p(x) =

 



 

1/2 ; 1 ≤ x < 2 1/4 ; 2 ≤ x < 4

0 ; x < 1 || 4 ≤ x

(16.8)

double get_myrandom(void){

double x;

double r;

r=uniform();

if( r < 0.5 ){

x= 2.0*r + 1.0;

} else {

x= 4.0*r ; }

return x;

}

F(x)

x

x p(x)

・累積分布関数

F (x) =

 

 



 

 

0 ; x < 1 (x − 1)/2 ; 1 ≤ x < 2 (x − 2)/4 + 1/2 ; 2 ≤ x < 4

1 ; 4 ≤ x

(17.1)

・母平均

µ =

∫

2 1

x 2 dx +

∫

4 2

x

4 dx = x

²

4

¯¯ ¯¯

^x=2

x=1

+ x

²

8

¯¯ ¯¯

^x=4

x=2

= 9

4 (17.2)

・母分散

∫

_∞

−∞

x

²

p(x) dx =

∫

2 1

x

²

2 dx +

∫

4 2

x

²

4 dx = x

³

6

¯¯ ¯¯

^x=2

x=1

+ x

³

12

¯¯ ¯¯

^x=4

x=2

= 35 6 σ

²

= 35

6 − ( 9

4 )

2

= 37

48 (17.3)

・確率密度 p(x) に従う確率変数 X を区間 [0, 1) の一様乱数 r から作る方法，逆関数法 確率変数 X が区間 [0, 1) の一様乱数 r から X = f(r) と作れたとする。このとき

∫

b a

p(x) dx = P (a ≤ X < b) = P (f

⁻¹

(a) ≤ r < f

⁻¹

(b)) = f

⁻¹

(b) − f

⁻¹

(a) (17.4) という関係がある。ただし f

⁻¹

は f の逆関数を表す。式 (15.3) と比較すると，f

⁻¹

を X の累 積分布関数 F (x) とすればよいことがわかる。従って

X = F

⁻¹

(r) , F : 累積分布関数 (17.5)

によって確率変数 X が区間 [0, 1) の一様乱数 r から作れる。

計科 II.18 (例) コーシー (Cauchy) 分布

p(x) = a π

1

a

²

+ x

²

(18.1) F (x) = 1

π arctan ( x

a )

+ 1

2 (18.2)

X = a tan (

(r − 1 2 )π

)

(18.3)

F(x)

p(x)

x

x

a = 1 の場合の図 (例) 指数分布；稀にしか起こらない現象が次に起きるまでの時間の分布

p(x) = a e

^−ax

, x ≥ 0 (18.4)

F (x) = 1 − e

^−ax

(18.5)

X = − 1

a log(1 − r) (18.6)

x

x p(x)

F(x)

a = 1 の場合の図

【問】

確率密度 p(x)

p(x) = {

1 − x/2 ; 0 ≤ x < 2

0 ; x < 0 || 2 ≤ x (18.7)

に従う確率変数 X を一様乱数 uniform() から生成する関数を表す C のプログラム (の一部) を 書きなさい。

【答】

累積分布関数は

F (x) = x − x

²

4 (18.8)

となる。一様乱数を r で表すと，X は r = X − X

²

4 を X について解いて X = 2 − 2 √

1 − r (18.9)

と表される。以下は C で書かれたプログラムの例を示す；

double Finverse(double r){

double x;

x= 2.0 - 2.0 * sqrt( 1.0 - r);

return x;

}

double get_myrandom(void){

double x;

double r;

r=uniform();

x=Finverse(r); /* Finverse は累積分布関数の逆関数で, 上で定義される */

return x;

}

・ 2 つの確率変数 X と Y の間に関数関係 Y = f (X) がある場合の確率密度の関係 A = f (a) , B = f (b) であるとき，

P (a ≤ X ≤ b) = P (A ≤ Y ≤ B) =

∫

_B

A

p

_Y

(y) dy (19.1)

となる。ところが，関係 y = f(x) を用いて，積分変数を y から x に変えると

∫

B A

p

Y

(y) dy =

∫

b a

p

Y

(f(x)) df (x)

dx dx (19.2)

となるので結局

P (a ≤ X ≤ b) =

∫

_b

a

p

_Y

(f (x)) df (x)

dx dx (19.3)

という等式が得られる。この式と (15.1) を比べると X の確率密度 p

X

(x) と Y の確率密度 p

Y

(y) の間に次の関係があることがわかる；

p

_X

(x) = p

_Y

(f(x)) df (x)

dx (19.4)

【問】

摩擦のない水平面 (x-y 平面) 上の点 (0 , − a) から物体をいろいろな向きに滑らせる。物体をど の向きにも一様にランダムに滑らすとする。物体が x 軸を横切るときの x 座標を確率変数 X とする。X に対する確率密度 p(x) を求めなさい。

【答】

物体を滑らす向きと y 軸との角度を θ とする。物体が x 軸を横切る場合のみ考えると， θ は

− π/2 < θ < π/2 の範囲にあるので θ に対する確率密度 p

_θ

(θ) は p

_θ

(θ) = 1

π , − π

2 < θ < π

2 (19.5)

計科 II.20 となる。X = a tan(θ) の関係があるので式 (19.4) より

p(x) = a π

1

a

²

+ x

²

(20.1)

となる。式 (18.1) と比べると X はコーシー分布に従うことがわかる。

・2 つの確率変数に対する確率密度

2 つの確率変数 X と Y を考える。点 (X , Y ) が x-y 平面の領域 R に属する確率 P (

(X , Y ) ∈ R )

は確率密度 p(x, y) を用いて P

(

(X , Y ) ∈ R )

=

∫ ∫

R

p(x, y) dxdy (20.2)

と表される。ここで，右辺は 2 変数の関数 p(x, y) の領域 R での 2 重積分を表す。

【問】

1 回にジャンプする変位 S = X(t + 1) − X(t) が連続な値をとるランダムウォークを考える。S が確率密度 w(s) に従う場合を考える。時刻 t のウォーカーの位置 X(t) が a ≤ X(t) < b とな る確率は確率密度 p(x, t) を用いて

P (a ≤ X(t) < b) =

∫

_b

a

p(x, t) dx (20.3)

と表される。時刻 t + 1 のウォーカーの位置 X(t + 1) に対する確率密度 p(x, t + 1) と p(x, t) の 間の関係を表す式を求めなさい。

【答】

時刻 t + 1 のウォーカーの位置 X(t + 1) が a ≤ X(t + 1) < b となる確率 P (a ≤ X(t + 1) < b) は a ≤ X(t) + S < b となる確率に等しい。この確率は 2 つの確率変数 S と X(t) に対する確率 密度 w(s) p(x, t) を領域 a ≤ X(t) + S < b で積分すれば得られる；

P (a ≤ X(t + 1) < b) =

∫ ∫

a≤x+s<b

w(s) p(x, t) dxds (20.4) 積分変数 x の代わりに y = x + s を用いると

P (a ≤ X(t + 1) < b) =

∫

b a

dy

∫

_∞

−∞

ds w(s) p(y − s, t) (20.5) となる。一方，X(t + 1) に対する確率密度 p(x, t + 1) を用いると

P (a ≤ X(t + 1) < b) =

∫

_b

a

p(x, t + 1) dx (20.6)

であるので，式 (20.5) と (20.6) を比べて次の関係式が得られる；

p(x , t + 1) =

∫

_∞

−∞

w(s) p(x − s , t) ds (20.7)

【問】

1 回にジャンプする変位 S = X(t + 1) − X(t) が確率密度 w(s) に従うランダムウォークを考え る。時刻 t = 0 でウォーカーが位置 X(0) = X

0

から出発する場合の，任意の時刻 t = 0, 1, 2, · · · での母平均と母分散を求めなさい。

【答】

一般の時刻 t でのランダムウォークの母平均 µ(t) =

∫

_∞

−∞

xp(x, t) dx (21.1)

をプリント.13 と同様に計算してみよう。まず時間が 1 ステップ増えると µ(t) がどれだけ変化 するかを考える。式 (20.7) より

µ(t + 1) =

∫

_∞

−∞

xp(x, t + 1) dx =

∫

_∞

−∞

dx x

∫

_∞

−∞

ds w(s) p(x − s , t) ds (21.2) となる。積分変数 x の代わりに y = x − s を用いると

µ(t + 1) =

∫

_∞

−∞

dy

∫

_∞

−∞

ds (y + s) w(s) p(y , t)

=

∫

−∞

dy y p(y , t)

∫

_∞

−∞

ds w(s) +

∫

−∞

dy p(y , t)

∫

_∞

−∞

ds s w(s)

= µ(t) + µ

S

(21.3)

が得られる。ここで

µ

_S

= E(S) =

∫

_∞

−∞

ds s w(s) (21.4)

は 1 回にジャンプする変位 S の母平均を表す。µ(0) = X

0

なので

µ(t) = µ

_S

t + X

₀

(21.5)

となる。

次に母分散を計算してみよう。時刻 t + 1 でのランダムウォーカーの位置 X(t + 1) の 2 乗の期 待値は

E (

X(t + 1)

²

)

=

∫

_∞

−∞

x

²

p(x, t + 1) =

∫

_∞

−∞

dx x

²

∫

_∞

−∞

ds w(s) p(x − s , t) ds (21.6) となる。上と同様に積分変数 x の代わりに y = x − s を用いると次が得られる：

E (

X(t + 1)

²

)

=

∫

_∞

−∞

dy

∫

_∞

−∞

ds (y + s)

²

w(s) p(y , t)

=

∫

_∞

−∞

dy y

²

p(y , t)

∫

_∞

−∞

ds w(s) + 2

∫

_∞

−∞

dy y p(y , t)

∫

_∞

−∞

ds s w(s) +

∫

_∞

−∞

dy p(y , t)

∫

_∞

−∞

ds s

²

w(s)

= E

( X(t)

²

)

+ 2µ(t) µ

_S

+ E (

S

²

)

(21.7)

計科 II.22 ここで

E (

S

²

)

=

∫

_∞

−∞

ds s

²

w(s) (22.1)

である。式 (21.3) と式 (21.7) より σ

²

(t + 1) = E

(

X(t + 1)

²

) − µ(t + 1)

²

= E(X(t)

²

) + 2µ(t) µ

_S

+ E (

S

²

) − (

µ(t) + µ

_S

)

2

= E(X(t)

²

) − µ(t)

²

+ E (

S

²

) − µ

²_S

= σ

²

(t) + σ

²_S

(22.2)

が得られる。

σ

_S²

= E (

S

²

) − µ

²_S

(22.3)

は 1 回にジャンプする変位 S の母分散を表す。ランダムウォーカーの母分散は t が 1 増えるご とに σ

²_S

だけ増加することがわかる。。σ

²

(0) = 0 なので

σ

²

(t) = σ

²_S

t (22.4)

となる。

【問】

一回にジャンプする変位 S(t) = X(t + 1) − X(t) が時刻 t によって異なる確率密度 w(s, t) に 従うランダムウォークを考える。この場合には任意の時刻 t = 0, 1, 2, · · · での母平均と母分散 はどんな式で表されるだろうか？

【答】

上と同様に考えると

µ(t + 1) = µ(t) + µ

_S(t)

(22.5)

σ

²

(t + 1) = σ

²

(t) + σ

_S(t)²

(22.6) となる。ここで

µ

_S(t)

= E (

S(t) )

=

∫

_∞

−∞

s w(s, t) ds (22.7)

は S(t) の母平均，

σ

_S(t)²

= E (

(S(t) − µ

_S(t)

)

²

) )

=

∫

_∞

−∞

s

²

w(s, t) ds − (µ

_S(t)

)

²

(22.8) は S(t) の母分散を表す。これより

µ(t) =

t−1

∑

k=0

µ

_S(k)

+ µ(0) (22.9)

σ(t)

²

=

t−1

∑

k=0

σ

_S(k)²

+ σ(0)

²

(22.10)

と表せる。

・サンプル平均の母平均と母分散

1 回にジャンプする変位 S = X(t + 1) − X(t) が同じ確率密度 w(s) に従うランダムウォークを 考える。

X(t + 1) = X(t) + S(t) (23.1)

なので，ウォーカーが X(0) = 0 から出発する場合を考えると X(t) = S(0) + S(1) + · · · + S(t − 1) =

t−1

∑

k=0

S(k) (23.2)

となる。従って (21.5) と (22.4) より t 個の独立な同じ確率密度に従う確率変数の和

S(0) + S(1) + · · · + S(t − 1) (23.3) の母平均と母分散は，和の中のどれか 1 個の確率変数 { S(k) , k = 0, · · · , t − 1 } の母平均や母分 散の t 倍になることがわかる。

今，母平均が µ

₁

，母分散が σ

₁²

である確率密度に従う確率変数を考える。この確率変数 N 個 から作られるサンプル平均

X ¯

_N

= 1

N (X

₁

+ X

₂

+ · · · + X

_N

) (23.4) の母平均 µ

_N

と母分散 σ

_N²

は

µ

_N

= µ

₁

, σ

_N²

= 1

N σ

²₁

(23.5)

となる。分散の平方根は確率変数の分布のばらつきの大きさの程度を表す量なので，

N 回の試行から得られたサンプル平均のばらつきの大きさは 1

√ N に比例して減少し ていく。

確率変数 X が何かの測定値である場合を考えよう。いろいろな測定誤差のために得られた測定 値 X は真の値 µ

₁

からずれる。しかし，何回も測定を繰り返して平均値を求めることで真の値 µ

_N

= µ

₁

からのずれが小さくなる確率が高まる。また，サンプルの不偏分散から母分散を推定 することにより得られたサンプル平均 X ¯

_N

に含まれる誤差の大きさを見積もることができる。

＜注＞ 中心極限定理

サンプル平均 X ¯

N

の従う確率密度は N → ∞ で平均 µ

N

，分散 σ

_N²

の 正規分布 p(x) = 1

√ 2π σ

_N

exp (

− (x − µ

_N

)

²

2σ

_N²

)

(23.6)

に近づくことが知られている。

計科 II.24

・中心極限定理を利用して平均 0，分散 1 の正規分布を作る方法

r

₁

, r

₂

, · · · r

_N

をそれぞれ区間 [0, 1) の一様乱数とする。，式 (16.2)，(16.3) より r

_k

の平均と分散 は 1/2 と 1/12 なので

z =

( r

₁

+ r

₂

+ · · · + r

_N

N − 1

2

)/ √ 1

12N (24.1)

は平均 0，分散 1 の確率変数となる。中心極限定理より N が大きくなると z の従う確率密度 は平均 0，分散 1 の正規分布に近づく。実際上は N = 6 ぐらいで z は正規分布に近くなる。

＜注＞ 2 つの確率変数 X と Y が 独立 であるということ

2 つの確率変数 X と Y が独立であるとは任意の関数 f (X) と g(Y ) に対して E

(

f(X) g(Y ) )

= E (

f(X) )

E (

g(Y ) )

(24.2) が成り立つことを意味する。

X と Y が離散的な値をとる確率変数の場合，これは積事象の確率がそれぞれの事象の確率の 積となること

P (X = x && Y = y) = P

_X

(X = x) P

_Y

(Y = y) (24.3) を意味する。また，X と Y が連続的な値をとる確率変数の場合，これは X，Y の確率密度が

p(x , y) = p

_X

(x) p

_Y

(y) (24.4)

のように 2 つの関数の積となることを意味する。

今まで，考えたきたランダムウォークはウォーカーのいる位置 X(t) と次にジャンプする変位 S(t) が独立な確率変数の場合なので，式 (14.11) や式 (21.7) の結果は次のようにしても得られる；

E (

X(t + 1)

²

)

= E

(

(X(t) + S(t))

²

)

= E (

X(t)

²

+ 2X(t) S(t) + S(t)

²

)

= E

( X(t)

²

) + E

(

2X(t) S(t) )

+ E (

S(t)

²

)

= E

( X(t)

²

) + 2 E

( X(t)

) E

( S(t)

) + E

( S(t)

²

)

= E

( X(t)

²

)

+ 2 µ(t) µ

S

+ E (

S(t)

²

)

(24.5)

・独立でない確率変数の例 (例 1) つながった２つのコイン

2 つのコイン A,B を投げ A が表の場合 X = 1，裏の場合 X = − 1 とする。また B が表の場合 Y = 1，裏の場合 Y = − 1 とする。2 つのコイン A,B をつなげて A が表の場合，必ず B が表に なるようにする。この場合

P

_X

(X = 1) = P

_X

(X = − 1) = 1

2 , P

_Y

(Y = 1) = P

_Y

(Y = − 1) = 1

2 (25.1)

であるが

P (X = 1 && Y = 1) = 1

2 , P (X = 0 && Y = 0) = 1 2

P (X = 1 && Y = 0) = P (X = 0 && Y = 1) = 0 (25.2) となり，X と Y は独立ではない。

( 参考 )

下の図はイジング (Ising) モデルと呼ばれる強磁性体 ( 磁石になる物質 ) の確率モデルの計算例を示す。

(http://www.apph.tohoku.ac.jp/matsubara-lab/simIsing/Isng.html を参照 )

強磁性体は小さな分子磁石の集まりとみなせる。下の図では各セルの色で分子磁石が上を向いているか 下を向いているかを表している。左側の図は温度が高い場合で各セルの色はほぼ独立な確率変数とみな せる。右側の図は温度が低い場合で同じ色のセルが集まる傾向があり，各セルの色は周りのセルの色と 独立ではなくなる。同じ向きを向いた分子磁石が集まることによりこの物質は磁石になる。

温度が高い場合，各セルの色はほぼ独立な確率変数とみなせる。 温度が低い場合，同じ色のセルが集まる傾向があり，各セルの色は独立ではなくなる。

(例 2) 1 回にジャンプする変位 S(t) = X(t + 1) − X(t) の従う確率密度が位置 X(t) によって異 なるランダムウォーク

ウォーカーの位置が x であるとき S(t) の従う確率密度を w(s, x) とすると式 (20.7) は p(x , t + 1) =

∫

_∞

−∞

w(s, x − s) p(x − s , t) ds (25.3)

となる。

計科 II.26

・OpenGL による作画

[参考文献]

[1] エドワード・エンジェル，“OpenGL 入門”, ピアソン・エデュケーション (2002) [2] 床井浩平，“GLUT による「手抜き」OpenGL 入門”，

http://www.wakayama-u.ac.jp/˜tokoi/opengl/libglut.html

例 1：2 次元の静止画 /**

OpenGL 2 次元グラフィックスによる静止画 Time-stamp: "2005/11/25 金 21:04 hig"

を少し変更

*/

#ifdef __unix__ /* Unix の場合 */

#include <unistd.h>/* これは usleep を使うためだけ */

#else /* Windows の場合 */

#include <windows.h>/* これは Sleep を使うためだけ */

#define _USE_MATH_DEFINES /* VSN で M_PI などを使うにはこれがいる */

#endif

#include <stdio.h>/* scanf/printf を使う */

#include <stdlib.h>/* 乱数を使う */

#include <math.h>/* sqrt */

#include <GL/glut.h>/* OpenGL/glut を使う */

/** ウィンドウのデフォルトの縦横, 位置 */

int wdefault=400; /* 幅 width */

int hdefault=400; /* 高さ height */

int wposition=100; /* 左上隅の横方向の位置 */

int hposition=80; /* 左上隅の縦方向の位置 */

/** ウィンドウ内に設定される座標 */

double xmin=-50.0; /* 左端 */

double xmax=+50.0; /* 右端 */

double ymin=-50.0; /* 下端 */

double ymax=+50.0; /* 上端 */

int t; /* display() が呼ばれた回数 */

/** 四角形の頂点のデータ */

double xrect[]={-30.0,-30.0,+30.0,+30.0};

double yrect[]={+10.0,-10.0,-10.0,+10.0};

/* 関数プロトタイプ宣言 これらは自分で作った関数 */

void init(void);

void display(void);

int main(int argc, char **argv){

/* 決り文句 glut の初期化 */

glutInit(&argc,argv);

/* RGB カラー/ ダブルバッファリングを使う */

glutInitDisplayMode( GLUT_RGB | GLUT_DOUBLE );

/* つくるウィンドウの大きさ */

glutInitWindowSize(wdefault,hdefault);

/* つくるウィンドウの位置 (画面左上から) */

glutInitWindowPosition(wposition,hposition);

/* ウィンドウをつくる. ウィンドウ名を指定 */

glutCreateWindow("OpenGL Example 1");

/* 初期設定 */

init();

/* 決り文句. 画面を描く関数は display という名前 */

glutDisplayFunc(display);

/* 決り文句. イベント待ち. 今の場合, while(1) display(); みたいなもの */

glutMainLoop();

return 0;

}

/** 画面を描く関数 */

void display(void){

int i;

glClear( GL_COLOR_BUFFER_BIT ); /* 画面クリア */

/* 2 次元の場合の決り文句 */

glMatrixMode(GL_PROJECTION);

glLoadIdentity();

計科 II.28 /* 画面に描かれる範囲. 左端の座標, 右端の座標, 下端の座標, 上端の座標

残り 2 つは後で説明します */

glOrtho(xmin,xmax,ymin,ymax,-1.0,1.0);

/* 長方形を描く */

glBegin(GL_LINES);

glColor3d(1.0,0.0,0.0); /* RGB 値. 赤で描く */

glVertex3d( xrect[0], yrect[0], 0.0 );

glVertex3d( xrect[1], yrect[1], 0.0 );

glColor3d(0.0,1.0,0.0); /* RGB 値. 緑で描く */

glVertex3d( xrect[1], yrect[1], 0.0 );

glVertex3d( xrect[2], yrect[2], 0.0 );

glColor3d(0.0,0.0,1.0); /* RGB 値. 青で描く */

glVertex3d( xrect[2], yrect[2], 0.0 );

glVertex3d( xrect[3], yrect[3], 0.0 );

glColor3d(0.0,0.0,0.0); /* RGB 値. 黒で描く */

glVertex3d( xrect[3], yrect[3], 0.0 );

glVertex3d( xrect[0], yrect[0], 0.0 );

glEnd();

glutSwapBuffers();

printf("display() is called %d times.\n",t);

t++;

}

/** 初期化 */

void init(void){

glClearColor(1.0,1.0,1.0,0.0); /* 画面クリアに使う色の RGBA 値, 背景は白 */

t=1;

}

・図形のタイプ glBegin() の引数 mode に指定できる図形のタイプ

GL POINTS 点を打つ。

GL LINES ２点を対にして、その間を直線で結ぶ。

GL LINE STRIP 折れ線を描く。

GL LINE LOOP 折れ線を描く。始点と終点の間も結ぶ。

GL TRIANGLES / GL QUADS ３／４点を組にして、三角形／四角形を描く。

GL TRIANGLE STRIP / GL QUAD STRIP 一辺を共有しながら帯状に三角形／

四角形を描く。

GL TRIANGLE FAN 一辺を共有しながら扇状に三角形を描く。

GL POLYGON 凸多角形を描く。

図形のタイプ

例 1 の作画例 例 2 の作画例

計科 II.30 例 2：2 次元のアニメーション

/**

OpenGL 2 次元グラフィックスによるアニメーション

Time-stamp: "2005/11/22 火 21:04 hig"

を少し変更

*/

#ifdef __unix__ /* Unix の場合 */

#include <unistd.h>/* これは usleep を使うためだけ */

#else /* Windows の場合 */

#include <windows.h>/* これは Sleep を使うためだけ */

#define _USE_MATH_DEFINES /* VSN で M_PI などを使うにはこれがいる */

#endif

#include <stdio.h>/* scanf/printf を使う */

#include <stdlib.h>/* 乱数を使う */

#include <math.h>/* sqrt */

#include <GL/glut.h>/* OpenGL/glut を使う */

#define NV 4 /* 頂点の個数 */

/** ウィンドウのデフォルトの縦横, 位置 */

int wdefault=400; /* 幅 width */

int hdefault=400; /* 高さ height */

int wposition=100; /* 左上隅の横方向の位置 */

int hposition=80; /* 左上隅の縦方向の位置 */

/** ウィンドウ内に設定される座標 */

double xmin=-50.0; /* 左端 */

double xmax=+50.0; /* 右端 */

double ymin=-50.0; /* 下端 */

double ymax=+50.0; /* 上端 */

int t; /* 時刻 */

int tmax=100; /* 半周するのにかかる時間 */

/** 四角形の頂点のデータ */

double xrect[]={-30.0,-30.0,+30.0,+30.0};

double yrect[]={+10.0,-10.0,-10.0,+10.0};

/* 関数プロトタイプ宣言 これらは自分で作った関数 */

void init(void);

void idle(void);

void display(void);

int main(int argc, char **argv){

/* 決り文句 glut の初期化 */

glutInit(&argc,argv);

/* RGB カラー/ ダブルバッファリングを使う */

glutInitDisplayMode( GLUT_RGB | GLUT_DOUBLE );

/* つくるウィンドウの大きさ */

glutInitWindowSize(wdefault,hdefault);

/* つくるウィンドウの位置 (画面左上から) */

glutInitWindowPosition(wposition,hposition);

/* ウィンドウをつくる. ウィンドウ名を指定 */

glutCreateWindow("OpenGL Example 2");

/* 初期設定 */

init();

/* 決り文句. 画面を描く関数は display という名前 */

glutDisplayFunc(display);

/* 決り文句. することがない (idle) とき呼ぶ関数は idle という名前 */

glutIdleFunc(idle);

/* 決り文句. イベント待ち. 今の場合, while(1) display(); みたいなもの */

glutMainLoop();

return 0;

}

/** 画面を描く関数 */

void display(void){

int i;

double cosine, sine;

glClear( GL_COLOR_BUFFER_BIT ); /* 画面クリア */

/* 2 次元の場合の決り文句 */

glMatrixMode(GL_PROJECTION);

glLoadIdentity();

/* 画面に描かれる範囲. 左端の座標, 右端の座標, 下端の座標, 上端の座標

残り 2 つは後で説明します */

計科 II.32 glOrtho(xmin,xmax,ymin,ymax,-1.0,1.0);

cosine=cos(M_PI*t/tmax);

sine=sin(M_PI*t/tmax);

/* 青長方形を描く */

/* RGB 値. 赤 緑 青 つまり青で描く */

glColor3d(0.0,0.0,1.0);

glBegin(GL_POLYGON);

for(i=0; i<NV; i++){

glVertex3d(cosine* xrect[i] - sine* yrect[i], sine* xrect[i] + cosine* yrect[i],

0.0 );

}

glEnd();

#ifdef __unix__ /* Unix の場合 */

usleep(100000);/* usleep(n) は n マイクロ秒だけ待つマイクロ秒=0.000001 秒 */

#else /* Windows の場合 */

Sleep(100); /* Sleep(n) は n ミリ秒=だけ待つ. ミリ秒=0.001 秒*/

#endif

glutSwapBuffers(); /* 決り文句. 画面を更新 */

/** 時刻を更新 */

t++;

printf("This is in display(). t is %d\n",t);

}

/** 初期化 */

void init(void){

glClearColor(1.0,1.0,1.0,0.0); /* 画面クリアに使う色の RGBA 値, 背景は白 */

t=0;

}

/* 何もすることがないときに呼ばれる関数 */

void idle(void){

glutPostRedisplay(); /* display 関数を呼んで描き直す */

}

図形の変換

プログラム中に変換行列 M

₁

, M

₂

と図形の点を表すベクトル ~ r が M

₁

M

₂

~ r

の順序で出てきた場合，ベクトル ~ r は M

₁

M

₂

~ r と変換される。

・回転

glRotated(GLdouble angle, GLdouble x, GLdouble y, GLdouble z)

変換行列に回転の行列を乗じる。引数はいずれも GLdouble 型 (double と等価)。１つ目の引数

angle は回転角 (単位は度)，残りの３つの引数 x, y, z は回転軸の方向ベクトル。引数が float

型なら glRotatef() を使う。

例 3：例 2 の下記の部分を変更

'

&

$

%

cosine=cos(M_PI*t/tmax);

sine=sin(M_PI*t/tmax);

/* 青長方形を描く */

/* RGB 値. 赤 緑 青 つまり青で描く */

glColor3d(0.0,0.0,1.0);

glBegin(GL_POLYGON);

for(i=0; i<NV; i++){

glVertex3d(cosine* xrect[i] - sine* yrect[i], sine* xrect[i] + cosine* yrect[i],

0.0 );

}

glEnd();

'

⇓

&

$

%

glMatrixMode(GL_MODELVIEW);

glLoadIdentity(); /* 変換行列を単位行列にする */

theta=180.0*t/tmax;

glRotated(theta,0.0,0.0,1.0); /* z 軸のまわりの theta 度 の回転 */

/* 青長方形を描く */

/* RGB 値. 赤 緑 青 つまり青で描く */

glColor3d(0.0,0.0,1.0);

glBegin(GL_POLYGON);

for(i=0; i<NV; i++){

glVertex3d(xrect[i],yrect[i],0.0 );

}

glEnd();

計科 II.34 以下のように書き換えても良い

例 4

'

&

$

%

glMatrixMode(GL_MODELVIEW);

theta=180.0/tmax;

glRotated(theta,0.0,0.0,1.0); /* z 軸のまわりに theta 度 の回転 */

/* 青長方形を描く */

/* RGB 値. 赤 緑 青 つまり青で描く */

glColor3d(0.0,0.0,1.0);

glBegin(GL_POLYGON);

for(i=0; i<NV; i++){

glVertex3d(xrect[i],yrect[i],0.0 );

}

glEnd();

この例では関数 display が呼ばれるたびに，theta=180.0/tmax だけの回転行列が変換行列に乗 ぜられる。

・平行移動

glTranslated(GLdouble x, GLdouble y, GLdouble z)

引数はいずれも GLdouble 型 (double と等価)。３つの引数 x, y, z は現在の位置からの相対的 な移動量。引数が float 型なら glTranslatef() を使う。

º

¹

·

¸

/** 四角形の頂点のデータ */

double xrect[]={-3.0,-3.0,+3.0,+3.0};

double yrect[]={+1.0,-1.0,-1.0,+1.0};

例 5-1 の作画例 例 5-2 の作画例