積分上限ｘを変数とし，次の積分で定義される関数を誤差関数と呼びます．

　　erf(x)=2/√π∫(0,x)exp(-t^2)dt

関数の定義は書物によって微妙に違い，係数2/√πを除いた形で定義されることもあるので注意してください（例：岩波数学公式）．

　

　また，統計学の分野では，標準正規分布の確率密度関数を

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

として，累積分布関数

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

によって定義された関数がよく使われます（下側確率：ΦL）．誤差関数との間には

　　Φ(x)=1/2[1+erf(x/√2)]

の関係があります．

　

　累積分布関数の一般的な定義は下側確率なのですが，この定義も本によって異なり，

　　Φ(x)=∫(x,∞)φ(t)dt　　　（上側確率：ΦU）

や

　　Φ(x)=∫(-x,x)φ(t)dt　　　（中心確率：ΦC）

で定義される場合もあります．左右対称な確率分布では中心確率も使いやすい定義です．

　

　今回のコラムで取り上げるテーマは「誤差関数の近似式」ですが，中心確率で定義した場合，近似式

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

が，簡単なわりには誤差関数のかなりよい近似になるので，知っておくと便利です．

　

　いずれにしても，どのような流儀で定義されているかは実に様々であり，同じ「誤差関数」という名前でも定義が異なることがあるので，他の書物を参照するときは注意し，使う前によく確かめておくことが大切です．

　

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

　

【１】解析的公式（級数展開→項別積分）

　

　正規分布の累積分布関数の計算では，

　　ａ）解析的な式に基づいて級数を足し上げる方法を採るもの（精度優先）

　　ｂ）近似式を用いているもの（速度優先）

があります．

　

　前者の例としては，正規分布の密度関数φ（ｘ）の指数部分をテイラー展開（マクローリン展開）して積分すると

　　Φ(x)=1/2+1/√2π{x-x^3/2･1!･3+x^5/2^2･2!･5-･･･+(-1)^(2k+1)x^(2k+1)/2^k･k!･(2k+1)+･･･}

が得られます．

　

　あるいは同じことですが，超幾何関数を使うと，以下のように表現されます．

　　Φ（ｘ）＝１／２＋１／√２πｘ1Ｆ1(1/2,3/2,-x^2/2)

また，数値計算の桁落ちを避けるために，この式にクンマー変換を施すと

　　Φ（ｘ）＝１／２＋１／√２πｘexp(-x^2/2)1F1(1,3/2,x^2/2)

が得られます．

　

　収束が遅いことを厭わなければ，累積確率はこれらの展開式から精度よく計算することができます．しかし，これらの展開式は収束が遅く，計算式としてはあまり実用的ではありません．

　

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

　

【２】正規分布の代用として

　

　そこで，正規分布の累積確率の近似計算について考えてみることにしましょう．プリミティブな発想ですが，正規分布によく似た釣り鐘型曲線で，しかも積分値を求めやすい関数を使って代用するというアイディアが浮かびます．

　

［１］ロジスティック分布

　

　f(x)=exp(-(x-α)/β)/{1+exp(-(x-α)/β)}^2

　　　=1/4βsech^2{(x-α)/2β}

　

　F(x)=1/{1+exp(-(x-α)/β)}

　　　=1/2[1+tanh{(x-α)/2β}]

　

　-∞<x<∞

　mean=α　 (mode,medianとも)

　variance=β^2π^2/3

　

　累積分布関数F(x)がロジスティック曲線になる分布がロジスティック分布です．α=0,β=1の場合が標準ロジスティック分布

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

でその平均値０，分散はπ^2／３となります．

　

　ロジスティック分布も裾の重い分布の１つで（−２，２）の外の確率が０．２３８４，（−３，３）の外の確率が０．０９４９です．この分布の裾の重さはほぼ自由度８〜９のｔ分布に相当します．また，尖度は正規分布が３に対してロジスティック分布が４．２です→［補］

　

　正規分布もロジスティック分布もベル型曲線を描きますが，正規分布

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

の粗い近似式として，平均値０，分散１に規格化したロジスティック分布

　f(x)=exp(-πx/√3)/{1+exp(-πx/√3)}^2

　F(x)=1/{1+exp(-πx/√3)}

もしばしば用いられます．

　

　また，確認したわけではありませんが，√３ではなく１．７にすると正規近似の精度がよくなるという話もあります．

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

　F(x)=1/{1+exp(-πx/1.7)}

　

　この式を形式的に二項展開すると，

　F(x)=1/{1+exp(-πx/√3)}

　　　=1-exp(-πx/√3)+exp(-2πx/√3)-exp(-3πx/√3)+exp(-4πx/√3)-･･･

となるわけですが，

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

　　　　　=1-1/2･exp(-2x^2/π)-1/8･exp(-4x^2/π)-1/16･exp(-6x^2/π)+･･･

とは似ても似つきません．

　

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

　

［２］ワイブル分布

　

　f(x)=α/β(x/β)^(α-1)exp{-(x/β)^α)

　F(x)=1-exp{-(x/β)^α)

　

　0≦x＜∞

　mean=βΓ(1/α+1)

　mode=β(l-1/α)^(1/α)

　median=β(log2)^(1/2)

　variance=β^2｛Γ(2/α+1)-[Γ(1/α+1)]^2｝

　

　ワイブル分布は上式の関数形で与えられますが，この式でαは形状母数，βは尺度母数と呼ばれます．ワイブル分布ではシェイプパラメータαの値を変えると形が種々に変化し，α＝１のときこの分布は指数分布となり，α＝２の場合がレイリー分布です．αの値が小さいほどガンマ分布より非対称性の強い分布を与え，α＞３〜４（３．２５）の場合にはほぼ正規分布の代用となるような対称的な形になります．

　

　α＝２，β＝√(π/2)のレイリー分布（左右非対称）では，

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

ですが，

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

の指数：^(1/2)がつきません．α＝３．２５（左右対称）の場合も，

　F(x)=1-exp{-(x/β)^3.25}

となるだけです．

　

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

　

【３】近似公式

　

　答えを正確に求めることが可能ならばそれはすばらしいことですが，ときには近似のほうが適切な場合があり，その答えについてある程度のことが知りたいとき，完全なものにこだわる必要はありません．

　

　しかし，正規分布の代用曲線を用いる方法では近似の精度が粗くなってしまいます．そのため，正規分布の累積確率を求めるための近似式は多数提案されています．なかでも

　　ａ）ラプラスやシェントンによる連分数展開式

　　ｂ）ヘイスティングスの近似式

は収束が速くかつ精度も高い式としてよく知られています．これらについては［補］で取り上げるにとどめます．

　

　今回のコラムで紹介したい近似式はこれらではなく，Williamsの近似式（1946）：

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

です．

　

　小生の偏見かも知れませんが，この近似式は近似式の中でもとりわけ美しく感じられるものです．近似式は純粋な数理統計学者にとってはとても気持ちの悪いことで，とくに小さな自由度に関して信頼できませんから，近似式は受け入れられないと思う人も多いかもしれません．しかし，このような美しい近似式ならば，受け入れてもよいという人は少なからずいるはずです．

　

　ただし，ここで

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

です．中心確率で定義した方がエレガントですし，正規分布は左右対称ですから使いやすい定義でもあります．

　

　下側確率：

　　ΦL(x)=∫(-∞,x)exp(-t^2/2)dt

で定義することにすると，

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

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

が得られます．この関数と正しい値との差はｘ＝１．６の近くで最も大きい（最大絶対誤差：0.003）のですが，それでも相対誤差は１％を越えません．

　

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

　

　　ΦL(x)=∫(-∞,x)exp(-t^2/2)dt

と定義することにすると，

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

が得られますが，近似の程度を一歩進めると，

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

となります．前者同様，ｘ＝２．５付近で誤差は最大となるのですが，最大絶対誤差（0.0007）は前者よりも小さくなります．

　

プログラム例：

　

1000 '

1010 ' *** 正規分布の下側確率 ***

1020 '

1030 *NORMAL.SUB:

1040 UU=SQR(1-EXP(-.63662*Z^2)*(1+.0095642*Z^4))/2

1050 P0=.5+UU*(SGN(Z)):RETURN

1060 RETURN

　

　さらに，

　　ΦL(x)≒1/2+1/2[1-exp(-2x^2/π){1+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6}]^(1/2)

と続きますが，項数を多くとるほど最大誤差を与える点はｘの大きい方へ移動し，ｘの大きいところでは精度が良くなるとは限りません．

　

　また，ここまでくると複雑すぎて，初期モデル

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

のエレガントさは失われてしまった感があります．近似式の目的からすれば，Williams=山内の近似式：

　　ΦL(x)=1/2+1/2[1-exp(-2x^2/π){1+x^4(0.0055+0.0551/(x^2+14.4)}]^(1/2)

あたりが限界といったところでしょう．

　

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

　

　ところで，

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

はどのようにして得られたものなのでしょうか？

　

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

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

になることがわかりますが，これ以上の進展は望めそうにありません．また，

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

　　　　　=1-1/2･exp(-2x^2/π)-1/8･exp(-4x^2/π)-1/16･exp(-6x^2/π)+･･･

などと，二項展開しても埒があきません．

　

　今回のコラムではこの近似式について書いたのですが，どのようにして得られたものなのか，今のところまったくわかりません．しかし，経験的に得られたものであるはずはなく，何らかの理論に基づいていると思われます．宿題として遺すことになりましたが，この近似式は前々から妙に気になっていたので，後日，その由来を調べて（その２）として報告したいと考えています．

　

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

　

［補］裾の重い分布

　

　たとえば，ラプラス分布（両側指数分布）の尖度はβ2＝６ですから，この分布も裾の重い分布の１つで（−２，２）の外の確率が０．１３５３，（−３，３）の外の確率が０．０４９４です．長い裾をもつほど飛び離れた値をもつ確率は高くなるわけですが，ここで正規分布および分散を１に規格化した対称性分布・・・ｔ分布（自由度５，１０），両側指数分布，ロジスティック分布の両側５％点，１％点，０．１％点を比較してみましょう．

　

　　　　　　　　　　 ５％点　　　１％点　　　０．１％点

正規分布　　　　　　 1.96059　　 2.57527　　 3.28897

t分布(df=5)　　　　　1.98992　　 3.10464　　 5.1894

　　 (df=10)　　　　 1.99354　　 2.83264　　 4.09149

両側指数分布　　　　 2.1183　　　3.25635　　 4.88449

ロジスティック分布　 2.01983　　 2.91835　　 4.19025

　

　１％点で比べると，正規分布，自由度１０のｔ分布，ロジスティック分布，自由度５のｔ分布，両側指数分布ときて，規格化できないコーシー分布が最も長い裾をもつことがわかります．

　

［注］この数値は小生が数表を用いずに，近似計算で求めた値であるから信頼率は95%以下と思われる．

　

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

　

［補］ラプラスやシェントンによる連分数展開式

　

1000 '

1010 ' ** NORMAL DISTRIBUTION **

1020 '

1030 *NORMAL.PROB1:

1040 PI=3.14159

1050 AZ=ABS(UUN)

1060 ZZ=AZ*AZ

1070 PZ=EXP(-ZZ/2)/SQR(PI*2)

1080 NZ=7:ID=NZ:SN=1

1090 IF AZ<2 THEN GOSUB *SHENTON ELSE GOSUB *LAPLACE

1100 IF UUN<0 THEN PP=1-PP

1110 RETURN

1120 '

1130 *SHENTON:

1140 PP=0

1150 FOR IZ=1 TO NZ

1160　 PP=ID*ZZ/(2*ID+1+SN*PP)

1170　 SN=-SN:ID=ID-1

1180 NEXT IZ

1190 PP=.5-PZ*AZ/(1-PP)

1200 RETURN

1210 '

1220 *LAPLACE:

1230 PP=0

1240 FOR IZ=1 TO NZ

1250　 PP=ID/(AZ+PP)

1260　 ID=ID-1

1270 NEXT IZ

1280 PP=PZ/(AZ+PP)

1290 RETURN

　

　シェントンの連分数展開式は確率変数uunの小さいところで収束がよく，ラプラスの連分数展開式はuunの大きいところで収束がよいので，両式を使い分ける境界値はuun=2としています（１０９０行）．

　

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

　

［補］ヘイスティングスの近似式

　

　正規分布の片側確率を求めるための近似式は多数提案されていますが，なかでもヘイスティングスの近似式はよく知られています．

　

1300 '

1310 ' ** NORMAL DISTRIBUTION **

1320 '

1330 *NORMAL.PROB:

1340 AA1=.196854

1350 AA2=.115194

1360 AA3=.000344

1370 AA4=.019527

1380 W=ABS(UUN)

1390 P0=1+W*(AA1+W*(AA2+W*(AA3+W*AA4)))

1400 P0=P0^4

1410 P0=1-.5/P0

1420 P0=.5+(P0-.5)*SGN(UUN)

1430 PP=1-P0

1440 RETURN

　

　ヘイスティングスの近似式と呼ばれるものには，もっと精度の良い関数近似もあります．

　

1450 '

1460 ' ** NORMAL DISTRIBUTION **

1470 '

1480 *NORMAL.PROB2:

1490 AA1=.049867347#

1500 AA2=.0211410061#

1510 AA3=.0032776263#

1520 AA4=3.80036E-05

1530 AA5=4.88906E-05

1540 AA6=5.383E-06

1550 W=ABS(UUN)

1560 P0=1+W*(AA1+W*(AA2+W*(AA3+W*(AA4+W*(AA5+W*AA6)))))

1570 P0=P0^16

1580 P0=1-.5/P0

1590 P0=.5+(P0-.5)*SGN(UUN)

1600 PP=1-P0

1610 RETURN

　

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