確率 P ( x;t ) によるランダムウォークの表現

(1)

.

... 確率 P(x, t) によるランダムウォークの表現

樋口さぶろお

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

計算科学☆演習II L05(2013-05-08 Wed)

今日の目標 .

..

1 P(x, t) の定義と意味が説明できる .

2.. オイラー表現とラグランジュ表現が説明できる .

..

3 ラグランジュ表現をオイラー表現に書き替えら

れる http://hig3.net

樋口さぶろお (数理情報学科) L05確率P(x, t)によるランダムウォークの表現計算科学☆演習II(2013) 1 / 24

(2)

ランダムウォークの座標の平均値と分散 Quiz解説

ここまで来たよ

1... ランダムウォークの座標の平均値と分散

Quiz解説

2... 確率 P(x, t) によるランダムウォークの表現

確率 P(x, t)

ラグランジュ表示とオイラー表示 P(x, t) ^の漸化式

P(x, t) の初期条件

(3)

Quiz解答:確率シミュレーション

ソースコード 1: もどってくる確率

1 /∗ i n c l u d e など前略∗/

2

3 int g e t r a n d o m ( d o u b l e y );

4 d o u b l e g e t u n i f o r m ();

5

6 int m a i n (){

7 int d ; /∗ ^{シード} ∗/

8

9 int x ; /∗ ウォーカーの座標 ∗/

10 int x s t a r t =0; /∗ スタート時刻の座標 ∗/

11 int x t a r g e t =0; /∗ ゴールの座標 ∗/

12

13 int n ; /∗ 試行のカウンタ ∗/

14 int n m a x = 1 0 0 0 0 ; /∗ ^{試行の総数}=サンプルサイズ ∗/

15 int c o u n t ; /∗ 条件が成立した試行の個数 ∗/

16

17 int t ; /∗ ^{時刻} ∗/

(4)

18 int t m a x = 2 0 ; /∗ ^{最終時刻} ∗/

19

20 s c a n f ( " % d " ,& d );

21 s r a n d ( d );

22 c o u n t =0;

23 for ( n =0; n < n m a x ; n + + ) {

24 x = x s t a r t ;

25 for ( t =0; t < t m a x ; t + + ) {

26 x = x + g e t r a n d o m ( g e t u n i f o r m ( ) ) ;

27 }

28 if ( x == x t a r g e t ){

29 c o u n t ++;

30 }

31 }

32 /∗ 確率の推定値 ∗/

33 p r i n t f ( " % lf \ n " ,( d o u b l e ) c o u n t / n m a x );

34 r e t u r n 0;

35 }

36

37 /∗^{関数定義など後略} ∗/

(5)

int^{で済む変数を}double^{にしている}

intのほうが速いしメモリも少ししか使わないし正確なので,int^で済むものはintで.

x,t,サンプル内通し番号n,条件を満たす試行の個数count^は整数値しかとらないのでint^で済む. srand^の引数は unsigned int^と決まってる.

double a; ^に対してa++;などとすると意図しない結果になる. forloop ^{のカウンタに}double^{を使うのは超悪趣味}. ^遅いし,^誤差が蓄積する.

タイプキャスト(型変換)してない

int a=1; int b=3; ^のとき,a/b^は0.

((double)a)/((double)b)^ならdouble^{同士の演算になって} 0.33333となる. (double)a/b やa/(double)b でもOK.

(6)

割り算でタイプキャストするのが面倒だから最初からdoubleを使っておく,のは超悪趣味. できるところまでint^{で計算して},必要なところでタイプキャストする.

for^{ループの範囲を意識}

時刻 t= 0,1,2,3, . . . ,20 =T =tmax. 問題文に20 と書いてある. サンプル内通し番号=^試行番号n= 0,1,2, . . . , N−1 =nmax-1. ^サンプルサイズ N は問題文では指定してないが,^{何か自分で決めない} と,最後に count/(double)nmax^{を計算できないのでは}.

nとN,tとT を区別.

(7)

確率P(x, t)によるランダムウォークの表現確率P(x, t)

1... ランダムウォークの座標の平均値と分散 Quiz解説

確率 P(x, t)

(8)

ランダムウォーカーの座標の母分布 例: ^{ランダムウォーク}

X_t+1 =X_t+R_t+1 例えば,初期条件X₀= 0.

例えば,

R 確率

... 0

−1 0

0 ¹₃ 1 ²₃

2 0

3 0

... ...

⇝

X₀ 確率

... 0

−1 0

0 1

1 0

2 0

3 0

... ...

X₁ 確率

... 0

−1 0

0 ¹₃ 1 ²₃

2 0

3 0

... ...

X₂ 確率

... 0

−1 0

0 ¹₉ 1 ⁴₉ 2 ⁴₉

3 0

... ...

· · ·

時間 t →

(9)

一般に, 1名のウォーカーが歩いているとき,P(x, t)を次のように定義.

Xt 確率

... ...

−1 P(−1, t) 0 P(0, t) 1 P(1, t) 2 P(2, t) 3 P(3, t)

... ... .確率 P(x, t) ..

...

x: 座標(整数), t: 時刻(整数)

定義 P(x, t) =^時刻 t^に,^{ウォーカーが}x^{にいる確率}.

性質

+∞

∑

x=−∞

P (x, t) = 1. (t : ^任意 )

(10)

別の表

t\x · · · −1 0 1 2 · · · x · · ·

0 0 0 1 0 0 · · · 0 · · ·

1 0 0 ¹₃ ²₃ 0 · · · 0 · · ·

2 0 0 ¹₉ ⁴₉ ⁴₉ · · · 0 · · ·

...

t P(x, t)

...

(11)

確率P(x, t)によるランダムウォークの表現ラグランジュ表示とオイラー表示

確率 P(x, t)

(12)

ラグランジュ表現とオイラー表現

確率は忘れて,^{ウォーカーが大勢}(^下では6^人)^{いる状況を考えよう}.

ラグランジュ表現数式的

xm,t: ^{ウォーカー番号}m^番の,^時刻 t^の座標.

上の状況なら x_0,t= 1, x_1,t= 2, x_2,t= 2, x_3,t= 3, x_4,t = 1, x_5,t = 2.

C的

x[m] ^{ウォーカー番号}m^番の座標(^時刻 t^とともに,^{この変数を更新}) int x[6]; /*配列の宣言*/

または,

int x[]={-1,0,0,1,-1,0};/*^{配列の宣言兼代入}*/

(13)

オイラー表現

数式的 U(x, t): ^時刻 t^に,^座標 x にいるウォーカーの人数. 上の状況なら

U(0, t) = 0, U(1, t) = 2, U(2, t) = 3, U(3, t) = 1, U(他, t) = 0.

C^的 U[x]^座標 x にいるウォーカーの人数(^時刻 t ^{とともに更新}) int U[100]; /*^{配列の宣言}. 100=x^{座標の上限}*/

または

int U[]={0,2,3,1,0,0,...};/*^{配列の宣言兼代入}*/

(14)

ラグランジュ表現とオイラー表現の比較 ラグランジュ表現オイラー表現

座標の値 int^でもdouble^でも int^限定(配列の添字) ウォーカー

の区別

ありなし

得意な問

彼はどこ ? ^{そこに何人} ?

シューティング

自機,

ボスキャラ雑魚キャラ

ブロック崩し

弾ブロック

テトリス落下前落下後

ランダムウォーク (例

X

_t

P (x, t)

(15)

.Quiz(P(x,t)の意味) ..

...

次のうち,確率P(x, t)について正しいものをいくつでも答えよう. .

..

1 P(x, t) ^をx ^{について加えると}1^になる. .

2.. P(x, t) ^をt^{について加えると}1^になる. .

..

3 P(x, t) をx,tの両方について加えると1になる. .

..

4 P(x, t) ^{の値は座標で単位は}m^である. .

5.. P(x,1)の値は時刻で単位は秒である. .

..

6 P(x,1)の値は確率で単位はない.

(16)

U とP の対応

xm,t, U(x, t) ^で,ランダムウォークのサイズN のサンプルを表現することもできる.

x_1,t, x_2,t, . . . , x_N,t. 度数 U(x, t). ^性質

∑∞ x=0

U(x, t) =N.

相対度数 u(x, t) = _N¹U(x, t). ^性質

∑∞ x=0

u(x, t) = 1.

サンプルサイズ N →+∞ ^で, 相対度数 u(x, t) → ^母分布P(x, t)

と期待されるので,P(x, t) ^{のことを考えるのに},u, U ^{でイメージしても} いい.

(17)

確率P(x, t)によるランダムウォークの表現 P(x, t)の漸化式

確率 P(x, t)

(18)

P(x, t) の漸化式

ランダムウォークで,X_t+1 =X_t+R_t+1, R=

R ^確率

−1 q = 1−p

+1 p

とする.

P(x, t) ^{で書きたい}

時刻 tにx=x₀ にいる U(x₀, t)≃N×P(x₀, t)人のうち,時刻 t+ 1には,平均的には

U(x0, t)×p ^人が x0+ 1^に U(x₀, t)×q 人が x₀−1 に去るはず.

(19)

逆に考えると,時刻 t+ 1に, x0−1^から

U (x

₀

− 1, t) × p

人が x0 に x0+ 1^から

U (x

₀

+ 1, t) × q

人が x0 にやってくるはず.

P (x

₀

, t + 1) = pP (x

₀

− 1, t) + qP (x

₀

+ 1, t)

(20)

考え方1

t\x · · · −1 0 1 2 · · · x · · ·

0 · · · · · ·

1 · · · · · ·

2 · · · · · ·

...

t P(x, t)

...

P(x, t) の漸化式を適用

t\x · · · −1 0 1 2 · · · x · · ·

0 · · · · · ·

1 · · · · · ·

2 · · · · · ·

...

t P(x, t)

(21)

確率P(x, t)によるランダムウォークの表現 P(x, t)の初期条件

確率 P(x, t)

(22)

P(x, t) の初期条件

例1

X₀ 確率

... 0

−1 0

0 1

+1 0

... 0

→ P(x,0) =

{ 1 (x = 0) 0 (x ̸ = 0)

例2

X0 確率

... 0 0 ¹₂ ... 0

1

→P(x,0) =

 

 

 



1

2

(x = 0)

1

2

(x = 10)

0 ( ^他 )

(23)

.Quiz(ラグランジュ表現とオイラー表現)

..

...

(座標が整数値のみをとる離散型の)ランダムウォークを考える.

6^{羽のペンギンが},^座標x= 0,1,2, . . . ,9の範囲をランダムウォークする. ある時刻tに,x= 1 に2羽,x= 3 に3羽,x= 8 に1羽いるとする.

. ..

1 ラグランジュ表現を用いたとき,配列x[] のサイズはどれだけ必要か. ^また,^時刻tに配列の各要素はどのような値をとるか.

.

2.. オイラー表現を用いたとき,^配列 u[]のサイズはどれだけ必要か. また,時刻tに配列の各要素はどのような値をとるか.

.Quiz(P(x, t)の漸化式) ..

...

各時間ステップ tで確率pで−2,確率qで−1,確率rで+2,確率 1−p−q−rで0(その場にとどまる) だけ移動するランダムウォークを考える.

P(x, t) =u[x]からP(x, t+ 1) =unext[x]を定める式を書こう.

(24)

予習復習問題これからは毎週

金11:05^{締切の予習復習問題は} RaMMoodle

http://el.math.ryukoku.ac.jp/moodle→^計算科学II(講義) でやってね.

水13:30^{締切の予習復習問題は} C-learning http://asp.c-learning.jp/s/^{でやってね}.

連休プチテスト対応. 実際に締切があるのは, 2013-05-08水 (C-learning), 2013-05-15^水(C-learning), ^{そのあとふつう}.

演習の春のプチテスト2013-05-10^金2. ^案内参照. ^出題計画(^{って範囲が} それくらいしかない) hist1, rw1, rand2. 数値だけ変えた,よりは変化は大きい予定.

自宅で演習の課題をやろうVisual Studio^には Express Edition^という‘^無料版’があります. 数理情報学科の学生はDreamSpark 経由で Visual

Studio 製品を自宅で使えます.

確率 P ( x;t ) によるランダムウォークの表現

時間 t →

∑

P (x, t) = 1. (t : 任意 )

彼はどこ ? そこに何人 ?

ボスキャラ 雑魚キャラ

弾 ブロック

X

P (x, t)

U (x

− 1, t) × p

U (x

+ 1, t) × q

P (x

, t + 1) = pP (x

− 1, t) + qP (x

+ 1, t)

{ 1 (x = 0) 0 (x ̸ = 0)

 

 

 



(x = 0)

(x = 10)

0 ( 他 )

P (x, t) = 1. (t : ^任意 )

彼はどこ ? ^{そこに何人} ?

ボスキャラ雑魚キャラ

弾ブロック

0 ( ^他 )