標準正規分布

　　φ(x)=1/√(2π)exp(-x^2/2)

において，累積分布関数を，

　　Φ(x)=∫(-x,x)f(t)dt

と定義する．もちろん，一般的な定義は

　　Φ(x)=∫(-∞,x)f(t)dt

なのだが，以下で取り上げる近似式との関係を考慮すると，前者の方がエレガントになるからである．

　

　Φ(x)の近似式に，

　　F(x)={1-exp(-2x^2/π)}^(1/2)

があり，簡単な割にはよい近似が得られる．

　

　ところで，昔，大学の物理の先生が言ったことであるが，近似式の「発明」はかなり泥臭い分野であるそうだ．要するに結果オーライで，発見に至った手法は問わず，実験結果と一致していればよいというわけである．

　

　数学の近似式でも，そのような例はあるはずだ．むしろ，経験的あるいは力づくで得られた近似式の方が多いのではないだろうかと思われる．たとえば，Mathematica関連書籍に載っている記事で，

　　Sin[x]+ArcSin[x] = x

というのがある．この式はむろん正しくないのだが，グラフを描画しているときに，たまたま上記の式が成り立つようなグラフが得られたそうだ．

　

　近似式と似たような問題に，電子回路の等価回路の問題がある．トランジスタなどは非線形回路なので数学的に扱いにくいので，回路シミュレーションなどでは，線形回路を組み合わせた等価回路に置き換えられる．しかし，この置き換えの法則たるや全くないようである．

　

　これらの指摘の如く，確かに近似式には泥臭いものが多い．しかし，

　　F(x)={1-exp(-2x^2/π)}^(1/2)

のような美しい近似式にはめったにお目にかかれない．まずないと言ってもよいものであろう．小生の偏見かもしれないが，このように美しい式が経験的に得られたものであるはずはなく，何らかの理論に基づいているに違いないのである．

　

　妙に気にかかる式なので，この近似式の導出過程を考えてみることにした．正規分布は熱方程式（拡散方程式）の解であるので，まず熱方程式との関係が頭に浮かび，そこから近似式

　　F(x)={1-exp(-2x^2/π)}^(1/2)

が求められたのではないかと考えた．熱方程式と確率論とは密接な関わりがあるので，あながち荒唐無稽ともいえないと思われたからである．ところが，熱方程式云々というアイディアでは結局うまくいかなかった．

　

　そこで，論文を取り寄せて調べてみたところ，この近似式は熱方程式とは一切関係がなく，非常に簡単な原理に基づいていて簡潔な導出が可能であることが判明したのである．

　

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

　

【１】Williamsの近似式

　

　統計数値表（日本規格協会）には，近似式

　　Φ(x)≒ {1-exp(-2x^2/π)}^(1/2)　　　（ｘ≧０）

　　　　 ≒-{1-exp(-2x^2/π)}^(1/2)　　　（ｘ≦０）

は，Williamsの近似式として紹介されている．

　

　一方，Continuous univariate distributions (John & Wiley)では，Polya（1945）によって得られたと記されていて，また，Williams=山内の近似式：

　　Φ(x)≒[1-exp(-2x^2/π){1+2(π-3)/3π^2･x^4}]^(1/2)

の基になったと思われる

　　Φ(x)≒[1-exp(-2x^2/π-2(π-3)/3π^2･x^4)]^(1/2)

は，Cadwellの近似式（1951）とある．

　

　これらの原典を取り寄せてみると，Williamsの論文が最も簡潔に記されていたので，以下は，Williamsに従って近似式の導出過程について解説したい．

　

　ガウス積分

　　Φ(x)=1/√(2π)∫(-x,x)exp(-t^2/2)dt

を２次元に拡張し，２重積分

　(2π)^(-1)∫(-x,x)∫(-x,x)exp{-(t1^2+t2^2)/2}dt1dt2

を考えると，これより直ちに

　　{Φ(x)}^2

を得ることができる．この式の積分領域は１辺の長さが２ｘの正方形であり，その面積は(2x)^2である．

　

　一方，周辺分布がともに平均０，分散１の正規分布となる２次元正規分布

　　p(x,y)dxdy=1/2π･exp(-(x^2+y^2)/2)dxdy

において，x=rcosθ，y=rsinθと極座標変換すると，ヤコビアンは

　　D(x,y)/D(r,θ)=r

であるから

　　p(x,y)dxdy=rexp(-r2/2)dr*1/2πdθ

　

　この式は，ｒとｒ＋ｄｒの間に落ちる確率が

　　rexp(-r2/2)dr

すなわち，振幅ｒの確率分布はレイリー分布となること，一方，位相θの分布は

　　p(θ)=1/2π

すなわち一様分布となることを意味している．

　

　ここで，１辺の長さが２ｘの正方形の積分領域を円で近似することを考える．

　　πr^2＝(2x)^2

すなわち，

　　r=2x/√π

にｒを定めれば，当該の正方形と面積の等しい円が得られることになる．

　

　このとき，

　　{Φ(x)}^2=(2π)^(-1)∫(-x,x)∫(-x,x)exp{-(t1^2+t2^2)/2}dt1dt2

　　　　　　 <∫(0,r)rexp(-r^2/2)dr∫(0,2π)1/2πdθ

　　　　　　 =∫(0,r)rexp(-r^2/2)dr=1-exp(-2x^2/π)

より，不等式

　　Φ(x)<{1-exp(-2x^2/π)}^(1/2)

が得られるのであるが，実際に計算してみると相対誤差は１％を越えず，

　　F(x)={1-exp(-2x^2/π)}^(1/2)

　　　　={1-exp(-.63662x^2)}^(1/2)

はΦ(x)のよい近似式になっているというわけである．

　

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

　

【２】近似式の精緻化

　

　F(x)を微分して密度関数を求めてみると，左右対称な関数

　　f(x)=|x|exp(-2x^2/π)/π{1-exp(-2x^2/π)}^(1/2)

が得られる．その分散は

　　σ^2=∫(-∞,∞)x^2f(x)dx

　　　　=2∫(0,∞)x^3exp(-2x^2/π)/π{1-exp(-2x^2/π)}^(1/2)dx

　　　　=-π/4∫(0,1)log(y)/(1-y)^(1/2)dy　　

　　　　=π(1-log2)=0.964

であるから，ｘ→ｘσと置き換えて，分散が１になるように規格化すると

　　F(x)={1-exp(-2x^2σ^2/π)}^(1/2)

　　　　={1-(e/2)^(-2x^2)}^(1/2)

　　　　={1-exp(-.613706x^2)}^(1/2)

となるが，これでは相対誤差は２％になるという．

　

　数値計算の結果，

　　F(x)={1-exp(-cx^2)}^(1/2)

という形では

　　F(x)={1-exp(-.6302x^2)}^(1/2)

が最良であり，相対誤差は約０．５％ということであった．

　

　一方，多項式を考えると，近似の精度をあげることができる．Cadwellは変数をリファインすることによって，近似の精度を高め，

　　Φ(x)≒[1-exp(-2x^2/π-2(π-3)/3π^2･x^4)]^(1/2)

　　　　 ≒{1-exp(-.63662*x^2-.0095642*x^4)}^(1/2)

さらに

　　Φ(x)≒{1-exp(-.63662*x^2-.0095642*x^4-.0005*x^6+.00002*x^8)}

を得ている．

　

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

　

【３】不等式の上界と下界

　

　２つの独立な標準正規変数x1,x2の同時分布，すなわち，２変数正規分布の議論から，正方形領域を半径２ｘ／√πの円で近似することを考えてきたが，正方形に内接する円（半径ｘ），外接する円（半径√２ｘ）を考えることによって，Φ(x)の上界・下界が，不等式

　　(1-exp(-x^2/2))^(1/2)≦Φ(x)≦(1-exp(-x^2))^(1/2)

によって与えられることが示される．

　

　さらに，左辺の1-exp(-x^2/2)は

　　1-exp(-x^2/2)+(2/π-1/2)^2exp(-x^2)

右辺の1-exp(-x^2)は

　　1-exp(-x^2)-(1-2/π)^2exp(-x^2)

と置き換えることができるとのことである．

　

　　(1-exp(-x^2/2)+(2/π-1/2)^2exp(-x^2))^(1/2)≦Φ(x)≦(1-exp(-x^2)-(1-2/π)^2exp(-x^2))^(1/2)

　

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

　

【４】不等式の高次元化

　

　以上述べてきたWilliamsの導出方法は，２次元正規分布を使ったものだが，同様のことを３次元で行うと，３次元空間の直角座標（ｘ，ｙ，ｚ）←→球面座標（ｒ，θ，φ）の座標変換は

　　x=rsinθcosφ,y=rsinθsinφ,z=rcosφ

ヤコビアンは

　　D(x,y,z)/D(r,θ,φ)=r^2sinθ

　

　ここで，方向を表すベクトルを球面座標でｓ＝（θ，φ）とおき，

　　ds=sinθdθdφ,dxdydz=r^2drds

のような変換を行えば，３次元正規分布

　　p(x,y,z)dxdydz=√(2/π)r^2exp{-r^2/2)dr*1/4πds

に変換され，ｒの確率分布はマクスウェル分布であること，また，ｓは球面上で確率密度1/4πの一様分布をすることも理解される．

　

　２次元の場合と同様にして，

　　4/3πr^3＝(2x)^3

　　r=2x/(4/3π)^(1/3)

なるｒを用いて，立方体を球で近似すると，不等式

　　{Φ(x)}^3<∫(0,r)√(2/π)r^2exp(-r^2/2)dr∫(0,4π)1/4πds

　　　　　　 =√(2/π)∫(0,r)r^2exp(-r^2/2)dr

が得られる．

　

　ここで，右辺の積分は，部分積分によって，

　　∫(0,r)r^2exp(-r^2/2)dr=-rexp(-r^2/2)+∫(0,r)exp(-r^2/2)dr

　　　　　　　　　　　　　 =-rexp(-r^2/2)+1/2Φ(r)

となり，再び，誤差関数が出現してしまう．これでは，２次元の場合のように簡単には表されないという不満が残ってしまうので，もうひとつだけ次元を上げることにしよう．

　

　４次元の場合，

　　1/2π^2r^4＝(2x)^4

　　r=2x/(1/2π^2)^(1/4)

なるｒを用いると，不等式

　　{Φ(x)}^4<1/2∫(0,r)r^3exp(-r^2/2)dr

が得られる．

　

　また，

　　∫(0,r)r^3exp(-r^2/2)dr=-r^2exp(-r^2/2)+2∫(0,r)rexp(-r^2/2)dr

　　　　　　　　　　　　　 =-r^2exp(-r^2/2)+2(1-exp(-r^2/2))

となって，偶数次元では明示的な公式が得られることがおわかり頂けよう．

　

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

　

　一般のｎ次元についても不等式の導出を試みたい．周辺分布が平均０，分散１の正規分布となるｎ次元正規分布は，

　　p(x1,x2,x3,･･･,xn)=(2π)^(-n/2)exp{-(x1^2+x2^2+･･･+xn^2)/2}

で与えられる．また，ｎ次元ユークリッド空間の点(x1,x2,x3,･･･,xn)は

　　ｒ＞0,0≦θ1,θ2,･･･,θn-2≦π,0≦θn-1≦2π

を満たすｒ,θ1,θ2,･･･,θn-1によって，

　　x1=rcosθ1

　　x2=rsinθ1cosθ2

　　x3=rsinθ1sinθ2cosθ3

　　････････････････････

　　xn-1=rsinθ1sinθ2･･･sinθn-2cosθn-1

　　xn=rsinθ1sinθ2･･･sinθn-2sinθn-1

と表すことができる．ただし，ｎ＝２のときは，周知のとおり，ｘ1＝ｒｃｏｓθ1，ｘ2＝ｒｓｉｎθ1とする．

　

　（ｒ,θ1,θ2,･･･,θn-1）がｎ次元極座標で，そのとき，ヤコビアンD(x1,･･･,xn)/D(r,θ1,･･･,θn-1)は

　　ｒ^(n-1)sin^(n-2)θ1･･･sin^2θn-3sinθn-2

となるので，

　　ds=sin^(n-2)θ1･･･sin^2θn-3sinθn-2dθ1dθ2･･･dθn-1

とおくと，　

　　dx1dx2･･･dxn=ｒ^(n-1)drds

　

　したがって，

　　p(x1,x2,x3,･･･,xn)dx1dx2･･･dxn=(2π)^(-n/2)exp(-r^2/2)r^(n-1)drds

　

　ここで，ｎ次元単位超球の体積を

　　ｖn=π^(n/2)/Γ(n/2+1)

とすると，ｎ次元単位超球の表面積は

　　ｓn-1=ｎｖn=(2*π^(n/2))/Γ(n/2)

であり，

　　p(x1,x2,x3,･･･,xn)dx1dx2･･･dxn

=nVn/(2π)^(n/2)exp(-r^2/2)r^(n-1)dr*1/nVnds

=1/(2^(n/2-1)Γ(n/2))exp(-r^2/2)r^(n-1)dr*Γ(n/2)/(2*π^(n/2))ds

が得られる．

　

　したがって，必要とされる密度関数は，ｎ次のχ分布

　　1/(2^(n/2-1)Γ(n/2))r^(n-1)exp(-r^2/2)

となり，このとき，

　　{Φ(x)}^n＜1/(2^(n/2-1)Γ(n/2))∫(0,r)r^(n-1)exp{-r^2/2}dr

　

　積分領域である１辺の長さが２ｘのｎ次元超立方体の体積は(2x)^nであり，半径ｒのｎ次元超球の体積は，ｎ次元単位超球の体積をｖnとすると

　　ｖnｒ^n

であるから

　　ｖnｒ^n=(2x)^n

　　ｒ=2x/ｖn^(1/n)=2xπ^(-1/2){Γ(n/2+1)}^(1/n)

となるようにｒを定めれば，当該の超立方体と体積の等しい超球が得られることになる．

　

［補］本コラムとは直接関係のない話であるが，似たような式があるので掲げておきたい．ハウスドルフ測度Ｈはルベーグ測度Ｌを一般化したものになっていて，ｎ次元単位球の体積をｖnとすれば，

　　Ｌ＝２^(-n)ｖnＨ＝２^(-n)π^(n/2)/Γ(n/2+1)Ｈ

の関係があることが知られている．

　

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

　

【５】失敗談

　

　お恥ずかしい話であるが，ときには失敗に至った道筋を掲げておくのも意味のあることだろう．近似式と熱方程式との関係をあげておきたい．

　

　母平均がμ，母分散がσ^2の正規分布はN(μ,σ^2)と書き表わされるが，

　　f(x)=1/√(2π)σ･exp{-(x-μ)^2/2σ^2}

の式は正規分布の確率密度関数と呼ばれるものである．

　

　年齢を固定したときの人間の身長の分布など，連続な値をとる測定値に対してヒストグラムを作成したときには単峰性で左右対称な形になることが多い．その場合，ヒストグラムにあてはまる関数として上式がよく想定される．

　

　正規分布の形状はμを位置母数，σを尺度母数として左右対称で左右に長くすそをひく釣り鐘型の分布曲線になる．そのグラフは山の高さがσに反比例して小さくなり，山の裾がσに比例して広がる．すなわち，本質的な意味での形状母数に相当するものはない．

　

　静かな水面にインクを一滴たらすとインクで染められた部分がどんどん拡散していくが，以下に述べるように，この濃度分布は正規分布においてσ^2を時間ｔに置換した式になっていて，この分布はインク分子が一定時間内に移動する距離の確率分布としても用いられる．

　　u(x,t)=1/√(4πt)exp(-x^2/4t)

すなわち，典型的な拡散過程では，時刻ｔには初期値から√ｔのオーダー離れた場所にいることを示しているのである．

　

　最も簡単な一次元熱伝導方程式（拡散方程式）

　　ｕt=ｕxx

は放物型偏微分方程式であり，初期条件をu(x,0)=F(x)とするとその解はフーリエ変換などを用いて解くことができ，正規分布関数と関数F(x)のconvolutionのような形で求められる．

　　u(x,t)=1/√(4πt)∫F(ξ)･exp{-(x-ξ)^2/4t}dξ

　

　ところで，積分変換解析法というと，まずラプラス変換やフーリエ変換，それに近年ではウェーブレット変換などが思い浮かぶが，積分変換解析法：

　 h(s)=∫k(s,x)f(x)dx

において，関数k(s,x)をカーネル（核関数）と呼ぶ．もっと，シンボリックに，

　　（出力）＝∫（入力）×（核）

と書いてもよいだろう．核は一般には複素数値関数であり，ラプラス変換のカーネルはk(s,x)=exp(-sx)，メリン変換ではk(s,x)=x^(s-1)というわけである．

　

　上の式において，

　　k(t,x,ξ)=1/√(4πt)exp{-(x-ξ)^2/4t}

は熱核と呼ばれるものである．熱核は２つの基本的な性質をもっている．１つはもちろん偏微分方程式であるが，もう１つは「近似」である．実は，よく知らないことなので明言することははばかられるが，積分変換には「再生核」と呼ばれるものがあり，これは変換前後でノルムを保存する等距離写像と関連しているらしい．

　

　正規分布のフーリエ変換は

　　∫(-∞,∞)exp(-πx^2)exp(-i2πxs)dx=exp(-πs^2)

より再び正規分布になるのであるが，これを用いたら誤差関数の近似が可能になるのではないか？　また，半直線上の熱方程式の解は

　　u(x,t)=∫(0,t)F(ξ)xexp(-x^2/4(t-ξ))/2√π(t-ξ)^(3/2)dξ

となるらしい．この式は，本コラムで問題とした関数

　　F(x)={1-exp(-2x^2/π)}^(1/2)

を微分した式

　　f(x)=|x|exp(-2x^2/π)/{1-exp(-2x^2/π)}^(1/2)

とよく似ている．これは偶然の一致だろうか？　拡散方程式と確率論（ブラウン運動）とは密接な関わりがあるので，偶然とはいえないのではなかろうか？

などと考えて，

　　斉藤三郎「再生核の理論入門」牧野書店

を読んでいたとき，頼んでいた論文が届き，このアイディアは見かけ倒しということがわかった次第である．

　

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

　

【６】思わぬ拾い物

　

　正弦積分とは，

　　Si(x)=∫(0,t)sin(t)/tdt=x-x^3/3･3!+x^5/5･5!-･･･

として定義される特殊関数（初等関数によって表し得ない関数）で，その特殊値

　　Si(∞)=∫(0,∞)sin(t)/tdt=π/2

は，ディリクレ積分と呼ばれる．

　

　コラム「シンク関数の数学的諸性質」でディリクレ積分の数学的諸性質，とくに（超）立方体の断面積との関係について取り上げた際，ａ，ｂ，ｃを正の定数０＜ａ≦ｂ≦ｃとするとき，

(1)2/π∫(0,∞)sin(ax)sin(bx)/x^2dx=a

(2)ａ＋ｂ≦ｃならば，

　　2/π∫(0,∞)sin(ax)sin(bx)sin(cx)/x^3dx=ab

(3)ａ＋ｂ＞ｃならば，

　　2/π∫(0,∞)sin(ax)sin(bx)sin(cx)/x^3dx=ab-1/4(a+b-c)^2

が成り立つ，などを紹介した．

　

　(1)の場合は積分値がｂ，(2)の場合は積分値がｃの値によらないことに注意していただきたいのだが，丹野修吉先生はこれらを超立方体と超平面の交わり部分の体積として証明していて，ａ＋ｂ≦ｃならば積分値がｃに依存しないことは，平面が立方体の上面に交わらないことに対応するもので，この驚くべき結果も超立方体と超平面の関係を考えると理解できるというものであった．

　

　斉藤三郎「再生核の理論入門」牧野書店を読んでいて，公式(1)：

　　K(x,y)=min(x,y)=2/π∫(0,∞)sin(xt)sin(yt)/t^2dt

と再生核の理論を用いると，関数φ(x)の逆関数{φ(x)}^(-1)が

　　{φ(x)}^(-1)=∫(0,b){min(φ(ξ),x)}'/φ'(ξ)dξ

=∫(0,b)(2/π∫(0,∞)sin(φ(ξ)t)sin(xt)/t^2dt)'/φ'(ξ)dξ

=2/π∫(0,∞)∫(0,∞)cos(φ(ξ)t)sin(xt)/tdtdξ

　

　したがって，例えば，単調関数φ(ξ)=ξ^nのとき，

　　x^(1/n)=2/π∫(0,∞)∫(0,∞)cos(ξ^nt)sin(xt)/tdtdξ

のように，逆関数が積分型で表現されることを知った．当該の近似式とは関係しなかったが，思わぬところでいい拾い物をした気分となった．

　

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