■格子上の確率論(その7)

 2年前に上梓した「雑然か整然か」は,配置が完全に規則的である結晶と,もう一方の極みにあるランダムな配置を扱ったコラムですが,その際に得た教訓は「結晶と完全にランダムな点の配置は対極的ではあるが,どちらも数学的な扱いが容易な構造である.」ということでした.数学的にみると完全な「ランダムネス」は意外に扱いやすいのです.
 
 シリーズ「格子上の確率論」では,ランダムな酔歩問題を正方格子・立方格子上の酔歩としてモデル化して扱いやすくしました.このような理想化したモデルはポリア・ウォークと呼ばれるものですが,今回のコラムではポリア・ウォークから規則的な格子の枠を外して「任意のピッチでδi進む.そこで任意の向きにθiターンする.」の繰り返しによる2次元酔歩(ピアソン・ウォーク)を取り上げたいと思います.
 
 今回のコラムで行う考察は,酔歩に関する論文:
  Montroll EW (1956): random walks in multidimensional spaces, especially on periodic lattices, J. Soc. Indust. Appl. Math. 4, 241-260
の後半部分に相当しますが,その論文の中に小生の愛読書である奥井重彦「電子通信工学のための特殊関数とその応用」森北出版の9.3.3節とまったく同じ内容のモデルを発見したことがきっかけになっています.
 
 そのとき,レイリー分布やライス分布のことが小生の脳裡をよぎったのですが,まったく無関係のように思われる酔歩と電気通信の間に,はたしてどのような関連性があるというのでしょうか?
 
===================================
 
【1】ピアソン酔歩の特性関数・確率密度関数
 
 論文の第7節はピアソン・ウォーク「任意のピッチでδi進む.そこで任意の向きにθiターンする.」する平面をガウス平面をみなして良い結果を得ているように思われますが,まず最初に,ある時刻(ステップ)における位置の分布について調べてみましょう.
 
 θの分布を一様分布:
  f(θ)=1/2π
とすると,最初の一歩のx座標:
  x=δcosθ
はさまざまの値をとります.そこで,特性関数→【補】を用いて,xの確率密度関数と累積分布関数を導き出すことにします.
 
 統計学では,通常,特性関数の定義として逆フーリエ変換のタイプが用いられますが,ここではフーリエ変換タイプを用いることにすると,特性関数は
  φ(u)=1/2π∫(0,2π)exp(-iux)dx
     =1/2π∫(0,2π)exp(-iuδsinθ)dθ
     =J0(uδ)
すなわち,ベッセル関数となることがわかります.
 
 確率密度関数は,特性関数を逆フーリエ変換したものですから,
  f(x)=1/2π∫(-∞,∞)exp(ixu)J0(uδ)du
    =1/π(δ^2-x^2)^(-1/2)=1/πδ{(1-(x/δ)^2}^(-1/2)
その累積分布関数は
  F(x)=1/2+1/2arcsin(x/δ)
のように求められます.
 
 次に,それぞれ向きが互いに独立で一様分布に従う2つのピアソン酔歩の和
  x=δ1cosθ1+δ2cosθ2 
の場合,和の分布の特性関数はそれぞれの分布の特性関数の積
  φx+y(t)=φx(t)*φy(t)
になりますから,
  φ(u)=J0(uδ1)J0(uδ2)
 
 最も簡単な場合として,ピッチは一定で,向きのみランダムに変化するランダムウォークから構成されている場合(δ1=δ2=δ)には
  f(x)=1/2π∫(-∞,∞)exp(ixu){J0(uδ)}^2du
    =1/2π^2δK[{(1-(x/2δ)^2}^(1/2)]
すなわち,第2種完全楕円積分となります.(δ1≠δ2の場合も第2種完全楕円積分になりますが,ここでは省略します.)
 
 ところで,前回のコラムでは1次元ランダムウォーク(ポリア酔歩)の母関数:
  (1-t^2)^(-1/2)
2次元ランダムウォークの母関数:
  2/πK(t)
とベッセル関数との間に,積分公式
  Σu2n=∫(0,∞)exp(-x){I0(x/d)}^ddx
が成り立つことを指摘しておきました.
 
 これは奇しくも本節と同じ状況ですが,単なる偶然とは思われません.それもそのはず,
  J0(ix)=I0(x),J0(x)=I0(-ix)
が成り立ち,変形ベッセル関数I0はベッセル関数J0の虚数部分を表すと考えてよいからです.
 
===================================
 
 一般に,n個のピアソン酔歩をステップ数を固定してアンサンブル(集合)的に観測したとき,ステップを合計したもののx座標の分布の特性関数は,ベッセル関数の積:
  φ(u)=ΠJ0(uδi)
したがって,確率密度関数は
  f(x)=1/2π∫(-∞,∞)exp(ixu)ΠJ0(uδi)du
と書けますが,この積分を直接行うことは簡単ではありません.
 
 しかし,ステップ数nが十分大きいときには
  φ(u)=ΠJ0(uδi) 〜 exp(-u^2σ^2/2)
のように近似できます.ただし,σ^2=1/2Σδi^2です.(δ1=δ2=・・・=δn=δのときσ^2=n/2δ^2)
 
 この場合に確率密度関数は
  f(x)=1/2π∫(-∞,∞)exp(ixu)exp(-u^2σ^2/2)du
    =1/√(2π)σexp(-x^2/2σ^2)
すなわち,平均値0,分散σ^2のガウス分布となります.
 
 このように,ステップ数nが十分大きいとき,位置の分布はガウス分布に従うことがわかりましたが,この問題では格子の形をランダムな構造に変えたところで,結果は格子上の酔歩の場合と同じはずです.そうでないと連続的なモデル(ブラウン運動)と矛盾するからです.
 
===================================
 
【2】ピアソン酔歩における移動距離の分布
 
 引き続いて,n個のステップの和の分布:
  f(x)=1/2π∫(-∞,∞)exp(ixu)ΠJ0(uδi)du
において,原点からの移動距離がrとr+drの間に落ちる確率を導きだしてみましょう.
 
 距離の分布のように非負の値をとるランダム変数の確率密度関数f(r),累積分布関数F(r)を導くには,ハンケル変換後の特性関数を用いるのが便利です.ハンケル変換後の特性関数をφ(λ)とすると
  f(r)=r∫(0,∞)λJ0(λr)φ(λ)dλ
  φ(λ)=∫(0,∞)J0(λr)p(r)dr
  F(r)=r∫(0,∞)J1(λr)φ(λ)dλ
の関係があります.ここで積分範囲が0〜∞となっていることに注意して下さい.
 
 ここでもまた,特性関数はベッセル関数の積
  φ(λ)=ΠJ0(λri)
と表されますが,前節と同様にして,nが大きい場合には,近似的に
  ΠJ0(λri) 〜 exp(-λ^2σ^2/2)   (σ^2=1/2Σri^2)
ですから,その確率密度関数は次のように求められます.
  f(r)=r∫(0,∞)λJ0(λr)exp(-λ^2σ^2/2)dλ
    =r/σ^2exp(-r^2/2σ^2)
すなわち,レイリー分布に帰着します.
 
 以上より,レイリー分布は同程度のピッチのピアソン酔歩の場合の原点からの距離の分布と考えることができるのですが,もしそこに歩幅がr0の一歩が加わったらどうなるのでしょうか? もちろん,その歩幅r0はほぼ一定の歩幅riの10倍もあるジャンプのことも,あるいは,フライトと呼んでも構わないような大ジャンプのことさえあるわけです.
 
 その場合,ハンケル変換後の特性関数は
  φ(λ)=J0(λr0)ΠJ0(λri)
したがって,nが大きい場合,
  φ(λ) 〜 J0(λr0)exp(-λ^2σ^2/2)
  f(r)=r∫(0,∞)λJ0(λr)J0(λr0)exp(-λ^2σ^2/2)dλ
    =r/σ^2exp{-(r^2+r0^2)/2σ^2}I0(rr0/2σ^2)
 
 この確率密度関数はライス分布と呼ばれるものですが,r0=0ならばI0(0)=1よりレイリー分布に帰着,また,r0》σならば,
  I0(x) 〜 exp(x)/√(2πx)
より,ライス分布は平均値r0,分散σ^2のガウス分布に近づきます.このことは直感的にも理解できるでしょう.
 
===================================
 
 以上のことは2次元平面で考えたものですが,3次元空間におけるピアソン酔歩も同様に考えることができます.
 
 小生の認識の誤りかもしれませんが,ベッセル関数は膜の振動だとか,円孔スリットによる回折にでてくる関数なので,平面的な関数というイメージがありました.そのため,この論文のピアソン酔歩を初めて読んだとき,上で求められた式は2次元でしか成り立たない,3次元以上ではベッセル関数はもはや出現しないのでは?という疑問を抱きました.
 
 しかし,幾度か読んでみるうちに,これらは3次元以上でも同じであって,特性関数はベッセル関数の積で表されるし,同様に,確率密度関数は逆正弦分布,楕円積分分布(と呼ぶことにする),ガウス分布になると考えるようになりました.それを確認するためには,単位円上のフーリエ解析ではなくて,いわば球面上におけるフーリエ解析のようなものが必要になると思われますが,このように考えないと連続的なモデル(ブラウン運動)の場合と矛盾するからです.
 
 さて,前にも述べましたが,ピアソン酔歩では格子の形をランダムな構造に変えたところで,結果はポリア酔歩すなわち格子上の酔歩の場合と同じはずです.そこで,任意のd次元における原点からの距離の分布についても述べておくことにします.
 
 d次元正規分布において,原点からのユークリッド距離の確率分布は,χ^2分布の平方根分布であるχ分布になります.自由度nのχ2分布の確率密度関数
  f(x)=1/{2^(n/2)Γ(n/2)}・(x)^(n/2-1)・exp(-x/2)   0≦x<∞
において,x=y2と変数変換すると,dx=2ydyより,χ分布の確率密度関数
  f(x)=1/{2^(n/2-1)Γ(n/2)}・(x)^(n-1)・exp(-x^2/2)   0≦x<∞
が得られます.
 
 コラム「標的問題の解とχ分布」では,周辺分布がともに平均0,分散σ^2の正規分布となるd次元正規分布を極座標変換して,rとr+drの間に落ちる確率の導き方を紹介していますので,是非ご参照下さい.
 
 とくに,自由度1のχ分布は,半正規分布(half normal distribution):
  f(x)=1/σ√(2/π)exp(-x^2/2σ^2)
であり,この分布は期待値が0の正規分布:
  f(x)=1/σ√(2π)exp(-x^2/2σ^2)
をy軸で折り返した分布になっています.
 
 また,自由度2のχ分布はレイリー分布:
  f(x)=x/σ^2exp(-x^2/2σ^2)
自由度3のχ分布はマクスウェル分布:
  f(x)=2^(3/2)/σ^3x^2exp(-x^2/2σ^2)
と命名されています.マクスウェル分布は3次元酔歩における原点からの遊走距離の分布になっています.
 
 一方,χ分布を拡張する方向としては,ひとつには自由度を増すこと,もうひとつは非心パラメータをつけることが考えられます.非心χ2分布の平方根分布が非心χ分布で,その確率密度関数は
  f(x)=1/σ2exp(-μ2/σ2)1/μ^(n/2-1)(x)^(n/2)exp(-x2/σ2)I(n/2-1,√μx/σ2)
で表されます.ここでI(ν,x)は第1種の変形ベッセル関数です.
 
 χ分布の自由度を増すと,半正規分布→レイリー分布→マクスウェル分布→・・・となりますが,非心χ分布の自由度を増すと,折り重ね正規分布→ライス分布→・・・が得られます.すなわち,非心χ分布において,自由度1の場合が折り重ねられた正規分布(folded normal distribution)で,期待値が0でない正規分布をy軸(x=0)で折り返して得られた分布になっています.また,自由度2の非心χ分布が上述したライス分布なのです.
 
===================================
 
【3】ピアソン酔歩の応用
 
 1920年代に,ウィーナーは離散時間の酔歩をブラウン運動の数学モデルとして定式化しましたが,その過程で,
  1)ブラウン粒子が出発する点が一度確定すると,その後の位置の分布はガウス分布に従うこと
  2)ブラウン運動の軌跡は連続的であるにも関わらず,どの点でも滑らかでないこと
を証明しています.
 
 さらに,ウィーナーはラジオの電波の干渉にブラウン運動のモデルを応用しました.実は今回のコラムに掲げた議論は,ピッチを振幅,ターンを位相,酔歩を電波に変えてもまったく同様に成り立ちます.電波が,振幅および位相がランダムに変化する正弦波から構成されている状況を考えてみることにしましょう.
 
 通信では電波が互いに干渉しあって,うねりのような強弱がついて非常に聴取しづらい現象が起こります.受信点では様々な伝搬経路を通ってきた多数の電波が合成されるからであり,この干渉をフェージング(fading)といいます.地表波と空間波の干渉を近距離フェージング,また,経路の異なる空間波どうしの干渉によるものを遠距離フェージングといいます.
 
 フェージングはその強度が等しいとき最も激しく起こりますが,合成された電波の位相ならびに振幅が不規則に変化するので,そのデータはランダムな性格を有し,明確な数式で記述されるよりは確率的記述と統計的平均値で表されなければなりません.このように,ラジオの電波の伝達を撹乱し放送局と局の間にダイヤルを合わせると聞こえる雑音には非常にランダムな性質がありますから,ウィーナーはバックグラウンド・ノイズからシグナルを分離する電子フィルターを設計する段階で,ブラウン運動の数学モデルを利用しています.その成果は第2次世界大戦中のレーダーの開発にも活かされ,長い間,軍事機密扱いにされていました.
 
 また,フェージングを受けた受信波包絡線の基本的な確率分布として,レイリー分布を一般化したχ分布が知られています.レイリー分布は,受信波が同程度の振幅の多重波からなるフェージングの包絡線変動を表すのに用いられますが,近年,χ分布は一般化されたレイリー分布(generalized Rayleigh distribution)として論文にも引用されることが多くなっています.χ分布はとくに電気通信分野で広い応用範囲を有して,その分野ではm分布とも呼ばれています.
 
 さらに定常波に多重波が加わる場合もありますが,その際に必要とされるのがライス分布です.ライス分布はマーカムのQ関数と関連していて,電気通信分野ではn分布の名で通っています.ライス分布において,μ=0ならばレイリー分布に一致します.また,μ→∞のときの極限は正規分布になります.
 
 ライス分布は第2次大戦中,米国ベル研のライスと本邦の仲上稔教授によって独立に研究されたものですが,とくに,仲上氏のフェージングの研究は戦時中のおそらく実験設備も十分でない環境にもかかわらず,先進国米国の研究に先立つオリジナルな成果として発表されています.しかし,originally made in Japan であるにもかかわらず,本邦ではその重要性が認められませんでした.ライス分布と名づけられていますが,この機会に仲上・ライス分布として紹介したいと考えています.
 
===================================
 
【補】積率母関数と特性関数
 
 関数f(x)としてsinxやxの無限積分を考えると前者は不定,後者は発散してしまいます.一方,これらの関数にexp(-x)を掛け合わせたf(x)exp(-x)を無限積分すると収束します.exp(-x)はxのどんなベキよりも速くx→∞で0になるからです.
 
 x>0のとき,ガンマ関数
  Γ(s)=∫(0-∞)x^(s-1)exp(-x)dx
の存在が知られているルーツにもこのような理由があるのです.すなわち,exp(-x)は無限積分において不定や発散する関数を収束させる働きをもっていることが理解されます.
 
 このことより,exp(-x)の代わりにもうひとつの変数sを含んだexp(-sx)を考え,
  F(s)=∫(-∞,∞)f(x)exp(-sx)dx
とおくと,無限積分の後,xの関数はsの関数に変換されます.この操作をラプラス変換と呼びます.
 
 ラプラス変換において変数sは複素変数であり,フーリエ変換はラプラス変換におけるパラメータsの実部が0である場合に相当します.フーリエ変換とは,普通の空間から周波数空間に観点を変えることに相当します.
 
 たとえば,周期が2πの三角関数をフーリエ変換してやると,その周期のところにピークがくるグラフになりますが,このことは波という観点からすべてをみることだと考えてもよいわけですし,また,波を粒子に変換しているという解釈も許されるのです.
 
[1]積率母関数
 
 h(x)=exp(tx)とおいた場合の期待値をtの関数M(t)とみて,積率母関数といいます.
  M(t)=E[exp(tx)]=∫(-∞,∞)f(t)exp(tx)dt
ラプラス変換のカーネルはh(x)=exp(-tx)ですから,積率母関数は1種のラプラス変換(+ラプラス変換)と考えることができます.
 
 また,
  M(t)=E[etx]=E[1+tx+(tx)^2/2!+(tx)^3/3+・・・]
  =1+tE[x]+t^2/2!E[x^2]+t^3/3!E[x^3]+・・・
より,
  M(t)=Σμ'k(t)^k/k!   k=0,1,2,・・・
を得ます.すなわち,積率母関数は原点まわりの積率の指数型母関数になっていることが理解されます.積率母関数を微分してt=0とおくことによって原点まわりの積率μ'kが順次求まります.
  d(k)M/dt(k)|t=0=μ'k=E[x^k]
 
(例題)正規分布の積率母関数は,M(t)=exp(μt+σ2t2/2)より,μ1,μ2を求めよ.
 
 M'(t)=(μ+σ2t)exp(μt+σ2t2/2)より,E[x]=M'(0)=μ
 M"(t)=(σ2+(μ+σ2t)^2)exp(μt+σ2t2/2)より,E[x^2]=M"(0)=σ2+μ2
を得ることができます.したがって,μ1=μ,μ2=σ2
 
 積率母関数には,「和の分布の積率母関数は積率母関数の積で表される.」という重要な性質があります.すなわち,x1,x2,...,xnが独立で,それぞれの積率母関数をMx1(t),Mx2(t),・・・,Mxn(t)とするとy=x1+x2+・・・+xnの積率母関数My(t)は
  My(t)=ΠMxi(t)
で表されるというものです.とくに,x1,x2,・・・,xnの積率母関数が同じを積率母関数Mx(t)をもつとき,My(t)=[Mx(t)]^nとなります.
 
 正規分布の和の分布について考えてみましょう.xがN(μx,σx2)に,YがN(μy,σy2)にしたがい,両者が独立であれば,x+yの積率母関数は
  Mx+y(t)=Mx(t)*My(t)=exp(μxt+σx2t2/2)exp(μyt+σy2t2/2)
  =exp((μx+μy)t+(σx2+σy2)t2/2)
 
 これはN(μx+μy,σx2+σy2)の積率母関数にほかなりません.したがって,正規分布の和の分布はまた正規分布となります.これを正規分布の再生性といいます.なお,ポアソン分布や負の2項分布,コーシー分布やガンマ分布も再生性を有しています.
 
 一方,差の分布の積率母関数は,
  Mx-y(t)=Mx(t)*My(-t)
で表されます.例題と同様に,正規分布の差の分布は
  Mx-y(t)=Mx(t)*My(-t)=exp(μxt+σx2t2/2)exp(-μyt+σy2t2/2)
  =exp((μx-μy)t+(σx2+σy2)t2/2)
すなわち,N(μx-μy,σx2+σy2)の正規分布になることを示すことができます.ところが,ポアソン分布の差の分布はポアソン分布にはならず,第1種変形ベッセル関数を用いて表される形になります.(ちなみに正規分布の積の分布は第2種変形ベッセル関数になります.)
 
[2]特性関数
 
 積率母関数の欠点は,積率をもたないコーシー分布やブラウンノイズ関数などに対しては積率母関数が定義されないということです.そこで,複素関数を導入して
  h(x)=exp(itx)
としたものが特性関数です.
  φ(t)=E[exp(itx)]=∫(-∞,∞)f(t)exp(itx)dt
 
 フーリエ変換のカーネルはh(x)=exp(-itx)ですから,特性関数は1種のフーリエ変換(+iフーリエ変換)と考えることができます.また,上式において,exp(itx)にオイラーの公式
  exp(itx)=cos(tx)+isin(tx)
を適用すると,特性関数は次のように表現できます.
  φ(t)=∫(-∞,∞)f(t)cos(tx)dt+i∫(-∞,∞)f(t)sin(tx)dt
 
 正規分布N(μ,σ2)の積率母関数は,M(t)=exp(μt+σ2t2/2)ですから,その特性関数はexp(iμt-σ2t2/2)となります.また,特性関数はすべての確率分布に対して存在し,コーシー分布の特性関数はexp(iμt-|t|σと表されます.
 
 また,畳み込みのフーリエ変換はフーリエ変換の単なる積になりますから,畳込みの特性関数はそれぞれの分布の特性関数の積
  φx+y(t)=φx(t)*φy(t)
で表されます.また,差の分布の特性関数は,
  φx-y(t)=φx(t)*φy(-t)
で表されます.
 
 特性関数を逆変換することによって,確率密度関数が得られます.たとえば,ガウス分布とコーシー分布の畳み込みを逆変換すること,すなわち,
  1/2π∫(-∞,∞)φx+y(t)exp(-itx)dt
によって,
  f(x)=1/π∫(0,∞)exp(-σ^2t^2/2-βt)cos{(γ-x)t}dt
が得られますが,この積分記号を含む関数はフォークト(Voigt)関数と呼ばれ,分光学の分野では非常に重要な関数となっています.
 
 積率に関しても,積率母関数と同様なことは特性関数でもいえて,
  φ(t)=Σμ'r(it)^r/r!
  μ'r=(-i)^rd(r)φ/dt(r)|t=0
が示されます.
 
(例題)特性関数も積率母関数同様に和の分布を求めるときなど利用されていますが,ここでは,区間(0,1)の一様分布の特性関数が
  φ(t)=exp(it/2)sin(t/2)/(t/2)
となることを利用して,一様乱数riをn個合計したものの分布が,n→∞の極限で正規分布になることを示してみましょう.
 
 一様乱数をn個の合計のしたものの分布の特性関数は
  [φ(t)]^n=exp(int/2){sin(t/2)/(t/2)}^n
 
 一方,シンク関数
  sinx/x=1-1/3!x2+5!x4-・・・(ベキ級数表示)
  =Σ(-1)^mx^2m/(2m+1)!
の解が±π,±2π,±3π,±4π,・・・となることを利用して,無限積表示すると
  sinx/x=(1-x2/π2)(1-x2/4π2)(1-x2/9π2)(1-x2/16π2)・・・
  =Π(1-x2/k2π2)  (無限積表示)
 
 kが大きいとき,
  (1-x2/k2π2)〜exp(-x2/k2π2)
  Π(1-x2/k2π2)〜exp(-x2/π2Σ1/k2)
ここで,Σ1/k2=π2/6=ζ(2)より,
  Π(1-x2/k2π2)〜exp(-x2/6)
これより,{sin(t/2)/(t/2)}^n→exp(-nt2/24)ですから,
  [φ(t)]^n→exp(int/2-nt2/24)
 
 正規分布N(μ,σ2)の特性関数はexp(iμt-σ2t2/2)ですから,この結果はn個の独立した一様乱数の和の分布は平均値n/2,分散n/12の正規分布に近づくことを示していて,これは,一様分布の場合について中心極限定理を示したものとなっています.
 
===================================