前回のコラムでは，３次元立方格子・４次元超立方格子状におけるランダムウォークの数え上げプログラムを利用して，再帰確率を数値的に計算してみました．その際，３次元では最初の２００項まで加算して0.32，４次元では１００項までで0.19という結果が得られました．

　

　その後にわかったことですが，４次元ランダムウォークの再帰確率は，

　　0.193201673224984...

ですから，最初の１００項まででもかなり近い数値が得られたと評価できるでしょう．

　

　ちなみに，３次元の再帰確率は0.32と計算されましたが，正確には，

　　0.340537329550999...

ですから，同じ項数で比較すると３次元のほうが収束が遅いことが理解されます．したがって，３次元ランダムウォークの再帰確率の精度を上げるためには，より多くの項数まで計算することが必要です．

　

　これらの計算では，漸近確率

　　ｕ2n　〜　２^(1-d)ｄ^(d/2)（πｎ）^(-d/2)

を求めてみて適当な項数で計算を打ち切りましたが，初等的な数え上げであっても，時間さえかければもっと精度よく計算できるものと思われます．しかし，計算に時間がかかるので，パソコンを用いたとしても第１０００項（ｕ2000）の計算などは相当に厄介です．

　

　ところで，前回のコラムに掲げたプログラムでは，アルゴリズムの把握のしやすさを優先させたため，余分な計算時間のかかるプログラムになってしまいました．そこで，今回のテーマとしては，計算所要時間が短くて済むようなものに改良したいと思います．たまにはプログラム演習を試みるのも如何でしょうか．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【１】配列を利用した３次元・４次元酔歩の数え上げプログラム

　

　項数を増やすことによって，計算所要時間は雪だるま式に急増するため，精度を上げるには１０倍も速い計算方法を考え出さなくてはなりません．そのためには，まず，配列を使った再帰確率の計算プログラムを取り上げます．

　

　前回のコラムに掲げたプログラムは，あまりにも真正面から馬鹿正直に組んでしまったため，対数計算が二度手間，三度手間になり，効率が悪いのが欠点でした．そこで，あらかじめ，階乗の対数値をメモリに取り込んでおけば，対数計算に要する時間を大幅に節約することができます．以下のプログラムは原型となるプログラムをわずか数カ所変更しただけですが，項数を増やすほど計算速度のアドバンテージが得られますから，改良前に比べて１０倍以上も計算速度が速くなっています．

　

（３次元）

1660 '

1670 ' ** 配列格納 **

1680 '

1690 NN=20:DIM G(NN*2)

1700 GGG=0:G(0)=0

1710 FOR I=1 TO NN*2

1720　 GGG=GGG+LOG(I)

1730　 G(I)=GGG

1740 NEXT I

1750 '

1760 SS=1

1770 FOR N=1 TO NN

1780　 GOSUB *U2NA

1790　'PRINT N,S

1800　 SS=SS+S

1810　 IF (N MOD 10)=0 THEN PRINT N,SS

1820 NEXT N

1830 PRINT SS,1-1/SS

1840 END

1850 '

1860 ' ** u2n **

1870 '

1880 *U2NA:

1890 S=0

1900 MAXI=N

1910 MINI=-INT(-N/3)

1920 FOR I=MAXI TO MINI STEP -1

1930　 IF (N-I)<I THEN MAXJ=N-I ELSE MAXJ=I

1940　 MINJ=-INT(-(N-I)/2)

1950　 FOR J=MAXJ TO MINJ STEP -1

1960　　 K=N-I-J

1970　　 GOSUB *TIE

1980　　'PRINT I;J;K;C

1990　　 LP=G(N*2)-N*2*LOG(6)-(G(I)+G(J)+G(K))*2+LOG(C)

2000　　 S=S+EXP(LP)

2010　 NEXT J

2020 NEXT I

2030 RETURN

　

（計算結果）

　　ｎ　　　　　Σｕ2n

　１００　　　　１．４６９９

　２００　　　　１．４８３４５

　３００　　　　１．４８９４８

　４００　　　　１．４９３０８

　５００　　　　１．４９５５４

　

　３次元酔歩の数え上げを５００項まで加算してみましたが，再帰確率は

　　Σｆ2n＝１−１／Σｕ2n＝０．３３

に微増しただけでした．もっとも，このことは漸近確率

　　ｕ2n　〜　２^(1-d)ｄ^(d/2)（πｎ）^(-d/2)

より予測されてはいたのですが，・・・

　

（４次元）

1760 '

1770 ' ** 配列格納 **

1780 '

1790 NN=20:DIM G(NN*2)

1800 GGG=0:G(0)=0

1810 FOR I=1 TO NN*2

1820　 GGG=GGG+LOG(I)

1830　 G(I)=GGG

1840 NEXT I

1850 '

1860 SS=1

1870 FOR N=1 TO NN

1880　 GOSUB *U2N4A

1890　'PRINT N,S

1900　 SS=SS+S

1910　 IF (N MOD 10)=0 THEN PRINT N,SS

1920 NEXT N

1930 PRINT SS,1-1/SS

1940 END

1950 '

1960 ' ** u2n **

1970 '

1980 *U2N4A:

1990 S=0

2000 MAXI=N

2010 MINI=-INT(-N/4)

2020 FOR I=MAXI TO MINI STEP -1

2030　 IF (N-I)<I THEN MAXJ=N-I ELSE MAXJ=I

2040　 MINJ=-INT(-(N-I)/3)

2050　 FOR J=MAXJ TO MINJ STEP -1

2060　　 IF (N-I-J)<J THEN MAXK=N-I-J ELSE MAXK=J

2070　　 MINK=-INT(-(N-I-J)/2)

2080　　 FOR K=MAXK TO MINK STEP -1

2090　　　 L=N-I-J-K

2100　　　 GOSUB *TIE

2110　　　'PRINT I;J;K;L;C

2120　　　 LP=G(N*2)-N*2*LOG(8)-(G(I)+G(J)+G(K)+G(L))*2+LOG(C)

2130　　　 S=S+EXP(LP)

2140　　 NEXT K

2150　 NEXT J

2160 NEXT I

2170 RETURN

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【２】漸化式を利用した１次元・２次元酔歩の数え上げプログラム

　

　１次元酔歩では，

　　ｕ2n＝2nＣn／２^(2n)　〜　（πｎ）^(-1/2)

　

　また，２次元酔歩では，２項係数に関する公式

　　ΣnＣknＣn-k＝Σ（nＣk）^2＝2nＣn

が成り立つので，

　　ｕ2n＝１／４^(2n)（2nＣn）^2　〜　（πｎ）^(-1)

と表されます．

　

　３次元酔歩ではｕ2nとΣｕ2nを求めるためには３重のfor･･･nextループ，４次元酔歩では４重のfor･･･nextループが必要になりましたが，２次元酔歩では２重のfor･･･nextループは必要ありません．そのため，３次元・４次元に比べて，１次元・２次元のランダムウォークの数え上げプログラムは格段に易しい問題となります．以下に，配列に用いて正攻法で組んだ１次元・２次元酔歩の数え上げプログラムを掲げます．

　

（１次元）

1000 '

1010 ' *** １次元酔歩 ***

1020 '

1030 ' save "u2n1

1040 '

1050 NN=5000:DIM G(NN*2)

1060 GGG=0:G(0)=0

1070 FOR I=1 TO NN*2

1080　 GGG=GGG+LOG(I)

1090　 G(I)=GGG

1100 NEXT I

1110 '

1120 SS=1

1130 FOR N=1 TO NN

1140　 LP=G(N*2)-N*2*LOG(2)-G(N)*2

1150　 S=EXP(LP)

1160　'PRINT N,S

1170　 SS=SS+S

1180　 IF (N MOD 1000)=0 THEN PRINT N,SS

1190 NEXT N

1200 PRINT SS,1-1/SS

1210 END

　

（２次元）

1000 '

1010 ' *** ２次元酔歩 ***

1020 '

1030 ' save "u2n2

1040 '

1050 NN=5000:DIM G(NN*2)

1060 GGG=0:G(0)=0

1070 FOR I=1 TO NN*2

1080　 GGG=GGG+LOG(I)

1090　 G(I)=GGG

1100 NEXT I

1110 '

1120 SS=1

1130 FOR N=1 TO NN

1140　 LP=G(N*2)*2-N*2*LOG(4)-G(N)*4

1150　 S=EXP(LP)

1160　'PRINT N,S

1170　 SS=SS+S

1180　 IF (N MOD 1000)=0 THEN PRINT N,SS

1190 NEXT N

1200 PRINT SS,1-1/SS

1210 END

　

　ところで，MS-DOS版のN88BASICでは，配列の大きさはいくらでも大きくできるわけではありません．６４ＫＢの限界値が規定され，そのため，単精度実数型配列の添字の最大値は１６３８３が上限となっています．配列を用いる限り，第２００００項（ｕ40000）の計算は不可能なのです．

　

　そこで，計算の工夫が必要になるわけですが，漸化式

　　ｕ2n＝ｕ2(n-1)・２ｎ・（２ｎ−１）／ｎ／ｎ／４

を利用することにします．すなわち，ｕ2(n-1)が計算済みであれば，ｕ0＝１から始めて，ｕ2，ｕ4，・・・と次々と計算していけばよいことになります．

　

　　u0=1,u2=1/2,u4=3/8,u6=5/16,u8=35/128,u10=63/256,u12=231/1024,

　　u14=429/2048,u16=6435/32768,u18=12155/65536,u20=46189/262144

　

（１次元）

1220 '

1230 ' ** 漸化式 **

1240 '

1250 NN=10000

1260 SS=1:S=1

1270 FOR N=1 TO NN

1280　 S=S*(N*2)*(N*2-1)/N/N/4

1290　'PRINT N,S

1300　 SS=SS+S

1310　 IF (N MOD 1000)=0 THEN PRINT N,SS

1320 NEXT N

1330 PRINT SS,1-1/SS

1340 END

　

　２次元酔歩の場合でも，１次元酔歩と同様の考え方により，ｕ0＝１を初めに求めておいて，後は次々にｕ2，ｕ4，・・・と求めていけばよいのです．

　

（２次元）

1220 '

1230 ' ** 漸化式 **

1240 '

1250 NN=10000

1260 SS=1:S=1

1270 FOR N=1 TO NN

1280　 S=S*(N*2)*(N*2)*(N*2-1)*(N*2-1)/N/N/N/N/16

1290　'PRINT N,S

1300　 SS=SS+S

1310　 IF (N MOD 1000)=0 THEN PRINT N,SS

1320 NEXT N

1330 PRINT SS,1-1/SS

1340 END

　

　これらの計算は，漸化関係を使ってデータが入力されるごとに更新しながら求めるものです．精度がよく，データの走査も１回だけですから，計算機向きの数値計算法といえます．以下に計算結果を示しますが，第２００００項まででもあっという間に計算が完了します．なお，この計算は収束せず，いずれ∞に発散してしまいます．

　

（計算結果）

　　　ｎ　　　　　Σｕ2n（ｄ＝１）　　　Σｕ2n（ｄ＝２）

　　２０００　　　　５０．４７２１　　　　３．４８５８４

　　４０００　　　　７１．３７１７　　　　３．７０６４１

　　６０００　　　　８７．４０９４　　　　３．８３５４６

　　８０００　　　１００．９３　　　　　　３．９２７０２

　１００００　　　１１２．８４２　　　　　３．９９８０５

　１２０００　　　１２３．６１２．　　　　４．０５６０８

　１４０００　　　１３３．５１５　　　　　４．１０５１４

　１６０００　　　１４２．７３３　　　　　４．１４７６４

　１８０００　　　１５１．３４１　　　　　４．１８５１２

　２００００　　　１５９．５８　　　　　　４．２１８６８

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【３】漸化式を利用した３次元・４次元酔歩の数え上げプログラム

　

　２次元酔歩では，２項係数に関して

　　ΣnＣknＣn-k＝Σ（nＣk）^2＝2nＣn

が成り立つので，計算が非常に簡単になりましたが，

　　Σ（nＣk）^3＝？，

　　Σ（nＣk）^4＝？

などには単純な公式のないことが証明されています．

　

　そのため，３次元酔歩では，６項分布において，ｉ＋ｊ＋ｋ＝ｎとして

　　ｕ2n＝（２ｎ）！／６^(2n)Σ１／（ｉ！ｊ！ｋ！）^2

を地道に計算していくしかありません．

　

　その係数部分を

　　ｖ2n＝（２ｎ）！／６^(2n)

とおくと，漸化式

　　ｖ2n＝ｖ2(n-1)・２ｎ・（２ｎ−１）／６／６

が成り立ちますので，この部分は芋ヅル式に計算できますが，問題は

　　Σ１／（ｉ！ｊ！ｋ！）^2

の部分です．

　

　では，Σ１／（ｉ！ｊ！ｋ！）^2はどうしたらよいのでしょうか？　この部分は芋ヅル式のようにスマートにはいかず，タコ足式のダサイ書式になってしまうのでしょうか．多項分布の周辺度数を固定して計算することになるのですが，この部分も配列を使わないで，漸化式を使うことにこだわって書いてみました．

　

（３次元）

2040 '

2050 ' ** 漸化式 **

2060 '

2070 NN=20

2080 SS=1:GG=0:GN=0

2090 FOR N=1 TO NN

2100　 GG=GG+LOG(N*2)+LOG(N*2-1)-2*LOG(6)

2110　 GN=GN+LOG(N)

2120　 GOSUB *U2NR

2130　'PRINT N,S

2140　 SS=SS+S

2150　 IF (N MOD 10)=0 THEN PRINT N,SS

2160 NEXT N

2170 PRINT SS,1-1/SS

2180 END

2190 '

2200 ' ** u2n **

2210 '

2220 *U2NR:

2230 I0=N:GI=GN

2240 J0=0:GJ=0

2250 K0=0:GK=0

2260 '

2270 S=0

2280 MAXI=N

2290 MINI=-INT(-N/3)

2300 FOR I=MAXI TO MINI STEP -1

2310　 IF (N-I)<I THEN MAXJ=N-I ELSE MAXJ=I

2320　 MINJ=-INT(-(N-I)/2)

2330　 FOR J=MAXJ TO MINJ STEP -1

2340　　 K=N-I-J

2350　　 GOSUB *TIE

2360　　'PRINT I;J;K;C

2370　　 GOSUB *FACTORIAL.X

2380　　 LP=GG-(GI+GJ+GK)*2+LOG(C)

2390　　 S=S+EXP(LP)

2400　 NEXT J

2410 NEXT I

2420 RETURN

2430 '

2440 ' ** 階乗計算 **

2450 '

2460 *FACTORIAL.X:

2470 F0=I0:F=I:GOSUB *FACTORIAL.SUB:GI=GI+GGG

2480 F0=J0:F=J:GOSUB *FACTORIAL.SUB:GJ=GJ+GGG

2490 F0=K0:F=K:GOSUB *FACTORIAL.SUB:GK=GK+GGG

2500 I0=I:J0=J:K0=K

2510 RETURN

2520 '

2530 ' ** 階乗計算 **

2540 '

2550 *FACTORIAL.SUB:

2560 IF F0=F THEN GGG=0:RETURN

2570 IF F0>F THEN LARGE=F0:SMALL=F ELSE LARGE=F:SMALL=F0

2580 GGG=0

2590 FOR H=SMALL+1 TO LARGE

2600　 GGG=GGG+LOG(H)

2610 NEXT H

2620 GGG=GGG*SGN(F-F0)

2630 RETURN

　

（４次元）

2180 '

2190 ' ** 漸化式 **

2200 '

2210 NN=20

2220 SS=1:GG=0:GN=0

2230 FOR N=1 TO NN

2240　 GG=GG+LOG(N*2)+LOG(N*2-1)-2*LOG(8)

2250　 GN=GN+LOG(N)

2260　 GOSUB *U2N4R

2270　'PRINT N,S

2280　 SS=SS+S

2290　 IF (N MOD 10)=0 THEN PRINT N,SS

2300 NEXT N

2310 PRINT SS,1-1/SS

2320 END

2330 '

2340 ' ** u2n **

2350 '

2360 *U2N4R:

2370 I0=N:GI=GN

2380 J0=0:GJ=0

2390 K0=0:GK=0

2400 L0=0:GL=0

2410 '

2420 S=0

2430 MAXI=N

2440 MINI=-INT(-N/4)

2450 FOR I=MAXI TO MINI STEP -1

2460　 IF (N-I)<I THEN MAXJ=N-I ELSE MAXJ=I

2470　 MINJ=-INT(-(N-I)/3)

2480　 FOR J=MAXJ TO MINJ STEP -1

2490　　 IF (N-I-J)<J THEN MAXK=N-I-J ELSE MAXK=J

2500　　 MINK=-INT(-(N-I-J)/2)

2510　　 FOR K=MAXK TO MINK STEP -1

2520　　　 L=N-I-J-K

2530　　　 GOSUB *TIE

2540　　　'PRINT I;J;K;L;C

2550　　　 GOSUB *FACTORIAL.X

2560　　　 LP=GG-(GI+GJ+GK+GL)*2+LOG(C)

2570　　　 S=S+EXP(LP)

2580　　 NEXT K

2590　 NEXT J

2600 NEXT I

2610 RETURN

2620 '

2630 ' ** 階乗計算 **

2640 '

2650 *FACTORIAL.X:

2660 F0=I0:F=I:GOSUB *FACTORIAL.SUB:GI=GI+GGG

2670 F0=J0:F=J:GOSUB *FACTORIAL.SUB:GJ=GJ+GGG

2680 F0=K0:F=K:GOSUB *FACTORIAL.SUB:GK=GK+GGG

2690 F0=L0:F=L:GOSUB *FACTORIAL.SUB:GL=GL+GGG

2700 I0=I:J0=J:K0=K:L0=L

2710 RETURN

2720 '

2730 ' ** 階乗計算 **

2740 '

2750 *FACTORIAL.SUB:

2760 IF F0=F THEN GGG=0:RETURN

2770 IF F0>F THEN LARGE=F0:SMALL=F ELSE LARGE=F:SMALL=F0

2780 GGG=0

2790 FOR H=SMALL+1 TO LARGE

2800　 GGG=GGG+LOG(H)

2810 NEXT H

2820 GGG=GGG*SGN(F-F0)

2830 RETURN

　

　漸化式を使ってコーディングする際，いくつかの書式・書法が考えられました．実はもっと技法的に凝って速度優先（記述の冗長さにこだわらないで，少々まどろっこしくても計算速度が速いもの）としたものも考えられたのですが，ここには簡潔であることをモットーにして書いたものを掲げました．漸化式を使ったものとしては，多分これ以上簡潔には書けないものと思われます．

　

　また，計算速度に関しては配列を利用したプログラムのほうが数倍速いと思われますが，漸化式を利用したプログラムではメモリ制限がないため，より多くの項数まで計算するのに向いています．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】ｄ次元酔歩の挙動

　

　１歩あたりεだけ移動する１次元酔歩では，ｎステップにおける位置ｘnの平均と分散は，２項分布より，

　　Ｅ（ｘn）＝ｎε（ｐ−ｑ），Ｖ（ｘn）＝４ｎε2ｐｑ

ε（ｐ−ｑ）は１歩あたりの平均移動量と考えることができます．

　

　対称単純ランダムウォーク，すなわち，ｐ＝ｑ＝１／２の場合は

　　Ｅ（ｘn）＝０，Ｖ（ｘn）＝ｎ

で与えられますが，ド・モアブル=ラプラスの定理から，ｎ回のステップののち，その人がｘ＝ｋのところにいる確率は，ｎを十分大にすると平均０，分散σ2=ｎの正規分布Ｎ（０，ｎ）に近づくことを示しています．したがって，原点からの距離が√ｎの領域内に納まる確率は６８．３％，２√ｎの領域内に納まる確率は９５．４％と計算されます．

　

　１次元酔歩が√ｎの領域に入っていることが示されましたが，これは拡散現象では一般的にいえることです．ｄ次元酔歩が互いに独立な１次元酔歩ｄ個から構成されると考えると，ｎステップののち，

　　Ｅ（Ｗn）＝（０，０，，・・・０）

　　Ｅ‖Ｗn‖^2＝ΣＥ｜Ｗk,n｜^2＝ｄ×（√ｎ／ｄ）^2＝ｎ

したがって，

　　Ｅ‖Ｗn‖≦√Ｅ‖Ｗn‖^2＝√ｎ

ですから，原点（０，０，，・・・，０）を中心として拡散し，原点からの距離が√ｎの領域内に納まると考えられます．

　

　格子上のランダムウォークは，ブラウン運動などの拡散モデルとしてよく知られていますが，格子のモデルはブラウン運動の離散化とみなすことができるので，局所的にみると離散的過程であっても，大域的にみると連続空間に分布した連続的なガウス分布とみることができます．

　

　ｄ次元正規分布において，ｄ個の確率変数Ｘ1，・・・，Ｘdの同時確率は

　　ｆ（Ｘ，μ，Σ）＝（２π）^(d/2)Σ^(-1/2)ｅｘｐ｛−（Ｘ−μ）’｜Σ｜^(-1)（Ｘ−μ）／２｝

によって与えられます．

　

　ここで，Ｘはデータベクトル，μは平均ベクトル，Σはｄ×ｄ次の分散共分散行列を指しますが，具体的に書くと，周辺分布がともに平均０，分散σ2の正規分布となる２次元正規分布は

　　p(x,y)=1/2πσ2･exp(-(x2+y2)/2σ2)

３次元正規分布は

　　p(x,y,z)=1/(2π)^(3/2)σ3exp{-(x2+y2+z2)/2σ2)

となります．

　

　多重ガウス積分は不完全ガンマ関数とになることから，原点からの距離がλ＝√ｎの超球内に納まる確率は自由度ｄのχ2分布で与えられることが計算されます．（この詳細についてはコラム「ｎ次元楕円の陰と影」を参照願います．）

　

　それによると，２次元酔歩の場合，自由度２のχ2分布はすなわち指数分布ですから，全体のｐ％がその内側に入るような超球の半径ｃは

　　（ｃ／λ）^2＝−２ｌｎ（１−ｐ）

によって与えられることがわかります．

　

　例えば，ｃ＝λのときｐ＝０．３９４，ｃ＝２λのときｐ＝０．８６５となり，半径が２λの円を描くとその内側には８６．５％が含まれることになります．これらは，それぞれ１辺が２λ，４λの正方形の内側に入る確率，

　　ｐ＝（０．６８３）^2＝０．４６６＞０．３９４

　　ｐ＝（０．９５４）^2＝０．９１０＞０．８６５

よりもさらに小さくなります．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

---