前回のコラムでは，立方格子上の３次元ランダムウォークの再帰確率（0.34）が宿題として残されました．今回のコラムでは，まず，横浜国大・今野紀雄先生よりご教示された解答を提示することにしますが，特性関数を用いる再帰確率の計算法が，小生の行った初等的数え上げ法とはあまりにもかけ離れた方法だったのには大変驚かされました．

　

　その後，３次元・４次元格子上のランダムウォークの再帰確率を初等的な方法で計算してみますが，その前に，今野先生とのメールのやりとりから紹介したいと思います．

　

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

　

　　今野紀雄先生

　

　はじめまして．そして，新年あけましておめでとうございます．佐藤郁郎＠宮城県立がんセンター・研究所・病理と申します．

　

　小生は数学を専攻している者ではありませんが，多少なりとも科学の一端を担いたいという気持ちから，先生のご著書をいつも楽しく読まさせていただいております．

　

　ぶしつけではありますが，本日は３次元ランダムウォークの再帰確率について，つねづね疑問に思っていたことを質問させていただきます．

　

　ナツメ社からのご著書「確率モデル」によりますと，再帰確率0.34とあり，また，先生が監訳されたＩ．ピーターソンの「カオスと偶然の数学」にも0.34とありますから，この数値は数学者の間では正しいものと認知されている数値なのでしょう．しかし，私にとっては0.34が何の説明もなしに唐突にでてくるので，この漸近確率がどのようにして求められたものなのか，初めてご著書を購読したときからの疑問でありました．

　

　「マルコフ連鎖から格子確率モデルへ」（今野紀雄訳）およびその参考文献にもある「ルベーグ積分から確率論」（志賀徳造）を参考にして，

　　（中略）→前回のコラム参照

と計算してみましたが，0.34には一致しませんでした．0.34はどのようにして計算された値なのでしょうか？　よろしくご回答お願い申し上げます．

　

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

　

　　佐藤 郁郎　様

　

　拙著をお読み頂き有り難うございます．

　

　執筆するときに参考にした書籍の関係する部分のコピーを送りますのでご覧下さい．但し，この本は英語なので多少読みづらいかもしれませんが，佐藤様はかなり周辺のことを詳しくお知りのようですので，問題ないかと思います．

　

　取り急ぎ以上，用件のみにて失礼致します．　今野紀雄

　

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

　

【１】特性関数を用いた再帰確率の計算

　

　今野先生にはご多忙のおり，いろいろご教示下さり誠にありがとうございました．今後ともご指導よろしくお願い申し上げます．

　

　さて，送られてきた参考書籍：

　　Durrett: Probability: theories and examples

にある結論だけ示しますと，ｄ次元超立方格子上のランダムウォークにおいては，

　　Σｕ2n＝（２π）^(-d)∫(-π,π)Re(1-φ(t))^(-1)dt

　　　　φ(t)：特性関数φ(t)

が成り立つというものです．

　

　とくに，３次元の場合は，

　　Σｕ2n＝（２π）^(-3)∫(-π,π)（１−1/3Σｃｏｓt）^(-1)dt

　　　　　＝（√6/32π^3）Γ(1/24)Γ(5/24)Γ(7/24)Γ(11/24)

　　　　　＝１．５１・・・

と評価され，したがって，

　　Σｆ2n＝（Σｕ2n−１）／Σｕ2n＝１−１／Σｕ2n

　　　　　＝０．３４・・・

と計算できるとあります．

　

　私は３次元ランダムウォークの母関数さえも思いつきませんでしたから，特性関数を用いた計算をみて，眼からウロコ状態になりました．

　

　また，その文献には，志賀徳造「ルベーグ積分から確率論」共立出版にあるのと同じ漸近確率

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

が掲げられていますから，今後の計算ではこの評価式を使用することにします．

　

　なお，紹介された文献によると，４次元〜９次元超立方格子上のランダムウォークの再帰確率については，

　　Kondo and Hara (1987)

に掲げられているとのことでした．

　

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

　

【２】３次元酔歩の再帰確率の数え上げプログラム

　

　ところで，前回のコラムでは初等的な数え上げプログラムによる数値計算を紹介しました．

　　Σｆ2n＝（Σｕ2n−１）／Σｕ2n＝１−１／Σｕ2n

したがって，Σｕ2nがわかればΣｆ2nが求まるという方針自体は正しかったようですが，それならば，なぜ計算は不調に終わったのでしょうか？

　

　Σｕ2nを求めるプログラムに間違いがあったのかも・・・等々，さんざん悩みましたが，プログラム自体は９９％正解を計算しておりました．

　

1000 ' save "u2n

1010 NN=100

1020 SS=0

1030 FOR N=1 TO NN

1040　 FFF=N*2:GOSUB *FACTORIAL:GN=GGG

1050　 GOSUB *U2N

1060　'PRINT S

1070　 SS=SS+S

1080　 IF (N MOD 20)=0 THEN PRINT N,SS

1090 NEXT N

1100 PRINT SS,1-1/SS

1110 END

1120 '

1130 ' ** u2n **

1140 '

1150 *U2N:

1160 S=0

1170 MAXI=N

1180 MINI=-INT(-N/3)

1190 FOR I=MAXI TO MINI STEP -1

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

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

1220　 FOR J=MAXJ TO MINJ STEP -1

1230　　 K=N-I-J

1240　　 GOSUB *TIE

1250　　'PRINT I;J;K;C

1260　　 FFF=I:GOSUB *FACTORIAL2:GI=GGG

1270　　 LP=GN-N*2*LOG(6)-(GI+GJ+GK)*2+LOG(C)

1280　　 S=S+EXP(LP)

1290　 NEXT J

1300 NEXT I

1310 RETURN

1320 '

1330 ' ** tie **

1340 '

1350 *TIE:

1360 IF I=J AND J=K THEN C=1:GOTO 1390

1370 IF I=J OR J=K THEN C=3 :GOTO 1390

1380 C=6

1390 RETURN

1400 '

1410 ' ** 階乗計算 **

1420 '

1430 *FACTORIAL:

1440 IF FFF=0 THEN GGG=0:RETURN

1450 GGG=0

1460 FOR H=1 TO FFF

1470　 GGG=GGG+LOG(H)

1480 NEXT H

1490 RETURN

1500 '

1510 ' ** 階乗計算 **

1520 '

1530 *FACTORIAL2:

1540 IF FFF=0 THEN GGG=0:RETURN

1550 GGG=0

1560 FOR H=1 TO FFF

1570　 GGG=GGG+LOG(H)

1580　 IF H=K THEN GK=GGG

1590　 IF H=J THEN GJ=GGG

1600 NEXT H

1610 RETURN

　

　３次元対称単純ランダムウォークでは，６項分布において，ｉ＋ｊ＋ｋ＝ｎとして

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

と計算されます．簡単なプログラムなので解説するまでもないと思いますが，階乗計算がオーバーフローしないように対数に直して計算し，その後,１２８０行で元に戻しています．

　

　このプログラムで最も工夫した点は１１９０行，１２２０行のループです．ｉ＋ｊ＋ｋ＝ｎとなる（ｉ，ｊ，ｋ）の組を順次探索しながら累積確率を求めているのですが，探索に時間がかかるため，対称性を利用して探索範囲を絞り込み，ｎ≧ｉ≧ｊ≧ｋ≧０となるような最小範囲の探索で済むようにしています．

　

　これによって，計算量を約１／６に減らすことができますから，その分，あとで因子Ｃを掛けることにします．同順位がない場合，Ｃにはその組合せ数

　　Ｃ＝３！／１！１！１！＝６

を割り当てます．

　一方，同順位の数があった場合，たとえば，ワンペアでは

　　Ｃ＝３！／２！１！＝３

スリーカードでは，

　　Ｃ＝３！／３！＝１

となり，Ｃにはそれぞれ３と１が割り付けられます．

　

　計算結果は，次のようになりました．

　　ｎ　　　　　Σｕ2n

　　２０　　　　．４１３８３

　　４０　　　　．４４３１６

　　６０　　　　．４５６３８７

　　８０　　　　．４６４３２６

　１００　　　　．４６９７６３

　

　これが１に満たないようでは到底，再帰確率＝０．３４を示すことはできませんが，ここで単純なミスに気がつきました．ｕ0＝１を加えるのをすっかり見落としていたのです．したがって，１０２０行を

　　1020 SS=1

と変更します．そうすれば，

　　ｎ　　　　　Σｕ2n

　　２０　　　　１．４１３８３

　　４０　　　　１．４４３１６

　　６０　　　　１．４５６３８７

　　８０　　　　１．４６４３２６

　１００　　　　１．４６９７６３

　

　第１００項までで

　　Σｕ2n＝１．４７

すなわち，

　　Σｆ2n＝０．３２

であって，ほぼ正解の値が得られます．要するにお粗末な数え上げ計算だったというわけですが，わずか２０項まででも，

　　Σｆ2n＝０．３０

と評価されている点は驚きです．

　

　さて，ここでは最初の１００項まで加算しましたが，それでは２００項まで計算したらどうなるでしょうか？　そこで，ｕ2nの漸近確率を調べてみることにしましょう．

　

　ｕ2nの漸近確率は

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

で表されます．すなわち，ｄ＝３では

　　ｕ200　〜　2.33291E-4

　　ｕ400　〜　8.24808E-5

したがって，ここで計算を打ち切ったとしても

　　|Σｕ400−Σｕ200|≦2.33291E-2

となるはずで，第１００項まででもかなりよい近似値が得られているものと考えられます．実際に計算してみることにしましょう．

　

　　ｎ　　　　　Σｕ2n

　１２０　　　　１．４７３７８５

　１４０　　　　１．４７６９１６

　１６０　　　　１．４７９４４４

　１８０　　　　１．４８１５４

　２００　　　　１．４８３３１

　

　予想通り，第２００項まで計算したとしても

　　Σｆ2n＝０．３２

はほとんど変わらないという結果が得られました．ただし，ｎの増加に伴って計算時間が雪だるま式に急増するため，数え上げ法ではこのあたりが限界かと思われました．

　

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

　

【３】４次元酔歩の再帰確率の数え上げプログラム

　

　前節に掲げた数え上げプログラムは，４次元ランダムウォークの再帰確率の計算にも簡単に拡張できます．

　

1000 ' save "u2n4

1010 NN=100

1020 SS=1

1030 FOR N=1 TO NN

1040　 FFF=N*2:GOSUB *FACTORIAL:GN=GGG

1050　 GOSUB *U2N4

1060　'PRINT S

1070　 SS=SS+S

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

1090 NEXT N

1100 PRINT SS,1-1/SS

1110 END

1120 '

1130 ' ** u2n4 **

1140 '

1150 *U2N4:

1160 S=0

1170 MAXI=N

1180 MINI=-INT(-N/4)

1190 FOR I=MAXI TO MINI STEP -1

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

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

1220　 FOR J=MAXJ TO MINJ STEP -1

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

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

1250　　 FOR K=MAXK TO MINK STEP -1

1260　　　 L=N-I-J-K

1270　　　 GOSUB *TIE

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

1290　　　 FFF=I:GOSUB *FACTORIAL2:GI=GGG

1300　　　 LP=GN-N*2*LOG(8)-(GI+GJ+GK+GL)*2+LOG(C)

1310　　　 S=S+EXP(LP)

1320　　 NEXT K

1330　 NEXT J

1340 NEXT I

1350 RETURN

1360 '

1370 ' ** tie **

1380 '

1390 *TIE:

1400 IF I=J AND J=K AND K=L THEN C=1:GOTO 1480

1410 IF I=J AND J=K THEN C=4　　　　:GOTO 1480

1420 IF J=K AND K=L THEN C=4　　　　:GOTO 1480

1430 IF I=J AND K=L THEN C=6　　　　:GOTO 1480

1440 IF I=J THEN C=12　　　　　　　 :GOTO 1480

1450 IF J=K THEN C=12　　　　　　　 :GOTO 1480

1460 IF K=L THEN C=12　　　　　　　 :GOTO 1480

1470 C=24

1480 RETURN

1490 '

1500 ' ** 階乗計算 **

1510 '

1520 *FACTORIAL:

1530 IF FFF=0 THEN GGG=0:RETURN

1540 GGG=0

1550 FOR H=1 TO FFF

1560　 GGG=GGG+LOG(H)

1570 NEXT H

1580 RETURN

1590 '

1600 ' ** 階乗計算 **

1610 '

1620 *FACTORIAL2:

1630 IF FFF=0 THEN GGG=0:RETURN

1640 GGG=0

1650 FOR H=1 TO FFF

1660　 GGG=GGG+LOG(H)

1670　 IF H=L THEN GL=GGG

1680　 IF H=K THEN GK=GGG

1690　 IF H=J THEN GJ=GGG

1700 NEXT H

1710 RETURN

　

　３次元に比べてネスティングが深くなっているので，計算に時間がかかりますし，漸近確率

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

において，ｄ＝４では

　　ｕ200　〜　2.02643E-5

　　ｕ400　〜　5.06607E-6

ですから，最初の１００項まで計算した結果を以下に示します．

　

　　ｎ　　　　　Σｕ2n

　　２０　　　　１．２２９６７８

　　４０　　　　１．２３４４６９

　　６０　　　　１．２３６１０６

　　８０　　　　１．２３６９３２

　１００　　　　１．２３７４３

　

　これより，

　　Σｆ2n＝０．１９

という結果が得られました．すなわち，４次元では３次元よりもさらに進める方向が多くなるため，偶然に出発点に戻ってくるのはますます簡単にはいかないことになります．この結果を文献（kondo and Hara）にある数値と比較してみたいものです．

　

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