（その３）では，近似式の高次元化，すなわち，超立方体を同じ体積の超球で置換する方法について掲げた．近似式は，いずれの次元でも

　　Φ(x)≒{1-exp(-2x^2/π･定数)(多項式)}^(1/n)

の形になったが，高次元化することによって式は複雑になり，近似精度もかえって悪化したことから，高次元化はまったく無意味であることがわかった．

　

　今回のコラムでは，近似式

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

の精緻化の問題を取り上げるが，このことから，近似式を精緻化するためには，２次元で素朴に考える方がよい結果をもたらすだろうと推測された．すなわち，

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

を形を利用して近似度を高めるのであるが，その際，話の流れからいっても，近似式は

　　Φ(x)≒{1-exp(-2x^2/π)(多項式)}^(1/2)

の形がふさわしいと考えられる．

　

　今回のコラムでは，このように設定された問題を，変分法を用いて解くことによって，最大絶対誤差が最小となる最良近似式がどのように表されるかを考えてみることにした．

　

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

　

【１】ベキ級数展開

　

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

の指数部分をマクローリン展開して項別積分すると

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

　　　　 =√(2/π)x{1-x^2/6+x^4/40-･･･}

が得られる．あるいは，同じことであるが，超幾何関数を使うと，

　　Φ(x)=√(2/π)x1F1(1/2,3/2,-x^2/2)

　　　　 =√(2/π)xexp(-x^2/2)1F1(1,3/2,x^2/2)

と表現される．

　

　この多項式を収束するまで加算するのは，コンピュータを使ったとしても結構骨の折れる作業である．そこで，近似式が必要とされるわけであるが，このベキ級数展開を有限のところで打ち切った多項式型近似式ではいかにもつまらない．そこで，求めるべき近似式として

　　Φ(x)≒{1-exp(-2x^2/π)(多項式)}^(1/2)

の形を想定するのである．

　

　一方，近似式

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

をベキ級数展開すると

　　exp(-2x^2/π)=1-2x^2/π+1/2･(2x^2/π)^2-1/6･(2x^2/π)^3+･･･

より，

　　F(x)=√(2/π)x{1-1/2･(2x^2/π)+1/6･(2x^2/π)^2-･･･}^(1/2)

　　　　=√(2/π)x{1-x^2/2π+5x^4/24π^2-･･･}

が得られる．

　

　　Φ(x)=√(2/π)x{1-x^2/6+･･･}

と

　　F(x)=√(2/π)x{1-x^2/2π+･･･}

を比較することにより，（ｘが小さいとき）

　Φ(x)>F(x)

になることが理解されるであろう．

　

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

　

　さて，近似式

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

は，２次元正規分布において正方形と面積の等しい円に置換することによって得られたものであるが，まず考えられるのは，面積ではなく，その内側に入る確率が等しくなるように円の半径を定めることによって，近似式を精緻化する方法である．

　

　１辺の長さが2Xの正方形の内側に入る確率は{Φ(x)}^2，半径ｒの円の内側に入る確率は

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

であるから，

　　1-exp(-r^2/2)={Φ(x)}^2

となるような関数r(x)を求めることに帰着されるが，この方法では再びベキ級数に戻り着いてしまうだけであって，それ以上の進展は望めない．

　

　次に考えつく精緻化としては

　　Φ(x)/F(x)=1+(1/2π-1/6)x^2+･･･

より

　　Φ(x)≒{1-exp(-2x^2/π)}^(1/2)･{1+(1/6-1/2π)x^2+･･･}

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

　

　さらに，

　　{1+(1/3-1/π)x^2+･･･}

を指数関数で近似して

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

とすることであるが，これでは実際によい近似になっているかどうかすら疑わしい．

　

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

　

【２】変分法の適用

　

　そこで，変分法を適用して，絶対誤差が最小となる近似式は何か考えてみることにする．もちろん，この条件を満たす関数は無数に存在するわけであるから，関数の候補を，y=f(x)を任意の多項式関数として，

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

に限定する．

　

　また，絶対値記号があると扱いにくいので，目的とする関数の関数，すなわち，汎関数を絶対誤差の２乗

　　I[y]={Φ(x)-F(x)}^2

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

とおくと，最大絶対誤差（の２乗）を最小にするような変分問題に帰着されることになる．

　

　ここで，関数ｙ＝ｆ（ｘ）の近傍における比較関数

　　y+δy=f(x)+αη(x)

を考え，これをI[y]に代入すればαの関数

　　I[y+αη]={Φ(x)}^2-2Φ(x){1-exp(-2x^2/π)･(f(x)+αη(x))}^(1/2)+1-exp(-2x^2/π)･(f(x)+αη(x))

が得られる．

　

　したがって，

　　∂I/∂α=-exp(-2x^2/π)η(x)[Φ(x){1-exp(-2x^2/π)･(f(x)+αη(x))}^(-1/2)-1]

α＝０で極値をとらなければならないという必要条件，

　　∂I/∂α|(α=0)=0

から，解は

　　f(x)=exp(2x^2/π)[1-{Φ(x)}^2]

になる．

　

　しかし，f(x)は多項式であることから，

　　exp(2x^2/π)=1+2/π･x^2+1/2･(2/π)^2･x^4+1/6･(2/π)^3･x^6+･･･

　　{Φ(x)}^2=(2/π)･x^2{1-x^2/6+x^4/40-･･･}^2

　　　　　　 =(2/π)･x^2{1-1/3･x^2+7/45･x^4-･･･}

　　1-{Φ(x)}^2=1-2/π･x^2+2/3π･x^4-14/45π･x^6+･･･

と，ベキ級数展開すると，たったこれだけの項数の計算で

　　f(x)=exp(2x^2/π)[1-{Φ(x)}^2]

　　　　=1+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6+･･･

になることがわかるのである．

　

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

　　f(x)=1+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6+･･･

　

　高木貞治の名言に「微分のことは微分でやりなさい」というのがあるが，ここでは「変分のことは変分でやりなさい」であって，変分は微分の親戚であることがおわかり頂けたであろうか？

　

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

　

　次に，関数の候補を

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

すなわち，多項式を指数関数の引数中に繰り込んだ場合，同様に変分によって

　　-2x^2/π･g(x)=log[1-{Φ(x)}^2]

　

　ここで，対数関数の展開は

　　log[1-{Φ(x)}^2]=-{Φ(x)}^2-1/2･{Φ(x)}^4-1/3･{Φ(x)}^6-･･･

であって，そこへ

　　{Φ(x)}^2=(2/π)x^2{1-x^2/6+x^4/40-･･･}^2

を代入すると，

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

　

　これより，

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

　　-2x^2/π･g(x)=-2/π･x^2+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6+･･･

となる．

　

　この計算も小生が数式処理ソフトを用いずに，手計算で求めた値であるから信頼率は50%以下である．検算しておくことをお勧めする次第であるが，f(x)の場合と同じ式が出現するのは単なる偶然とは考えにくい．多分，計算があっている証だと思われる．

　

　また，多項式を指数関数中に引き込んでも，

　　+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6+･･･

の式は共通している．これにより，Williams=山内の近似式

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

において，指数関数近似

　　(1+2(π-3)/3π^2･x^4)≒exp(2(π-3)/3π^2･x^4)

すると，Cadwellの近似式が得られることがわかる．

　

　なお，g(x)を定数関数として，

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

の形の最良近似式を得ることも考えられるが，ｃの値はｘに依存することになるので，上の議論からそれはできない相談であることもおわかり頂けたかと思う．

　

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

　

【３】統計数値表との比較

　

　ところで，統計数値表（日本規格協会）には，当該の多項式は

　　f(x)=Σ(-1)^k･a2k･x^(2k)

　　a2k＝1/k!2^k･4/π∫(0,π/4)(sec^2θ-π/4)^kdθ

であって，実際に計算すると

　　a0=1,a2=0,a4=2(π-3)/3π^2,a6=(7π^2-60π+120)/45π^3,･･･

すなわち，

　　k=0　→　Φ(x)≒{1-exp(-2x^2/π)}^(1/2)

　　k=1　→　Φ(x)≒{1-exp(-2x^2/π)}^(1/2)

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

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

が得られるとある．

　

　積分

　　　a2k＝1/k!2^k･4/π∫(0,π/4)(sec^2θ-π/4)^kdθ

において，sec^2θやπ/4がでてくるおおよその理由については検討はつくが，変分法の結果と同一のものであるかどうかはまだ検討していない．

　

　変分法では，このような積分を用いなくても同じ多項式を簡単に求めることができたが，奥深いところでこの積分が関係しているものと想像される．一体どのようにして求められた積分なのだろうか？　この解答は次の機会までの宿題として自分自身に課すことにしたい．

　

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

　

　さらに，統計数値表では，下側確率：

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

で定義して，近似式

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

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

　　Φ3(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)

の順に最大絶対誤差が小さくなっている様子を明らかにしている．

　

　それと同時に，項数を多くとるほど最大誤差を与える点はｘの大きい方へ移動し，ｘの大きいところでは精度が良くなるとは限らないことにも触れている．この最大誤差を与える点は，

　　{Φ(x)-F(x)}

あるいは同じことではあるが，

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

が極値となる点である．

　

　微分すると，

　　2Φ(x)Φ'(x)=exp(-2x^2/π){4x/πf(x)ｰf'(x)}

また，

　　Φ'(x)=2/√2πexp(-x^2/2)

より

　　Φ(x)･exp(x^2(2/π-1/2))={√(2/π)xf(x)-1/2･√(π/2)f'(x)}

となる．

　

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

の場合，極値では

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

となるが，この解は左辺の曲線：y=Φ(x)exp(x^2(2/π-1/2))と右辺の直線：y=√(2/π)xの交点からx=1.76と求めることができる．

　

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

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

と進んでも，左辺の曲線は一定で，右辺が原点を通る多項式に変わるだけである．

　

　このことから直ちに最大誤差を与える点がｘの大きい方へシフトする様子を窺い知ることができるだろう．

　

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