超幾何関数を用いた確率分布の計算（その３）

　方程式ｆ（ｘ）＝０が因数分解できる場合は，その根は簡単に求めることができますが，解析的に求まらない場合は数値計算法によらなければなりません．もちろん，問題の解としては可能な限り厳密なものが望まれますが，非線形方程式では，多くの場合，厳密な解析解を得るのは非常にやっかいであり，その点，数値的近似解ならアプローチしやすく手早く求めることができるため，たいていの非線形方程式の求根問題には逐次法（反復法）が用いられます．

　確率分布の計算においても，十分な精度が得られるまで反復して解を求めますが，コラム「超幾何関数を用いた確率分布の計算」その１・その２では，精度と速度のせめぎ合いについて述べてきました．数値計算では精度と速度の２つが要請されますが，一般的にいって，精度を上げれば速度が下がり，速度を上げれば精度が下がります．すなわち，この２つの条件は競合的かつ背反的で，一方を減らすと他方が増えるというトレードオフの関係にあり，同時に向上させることは難しいとされます．これまで掲げてきた確率分布の計算プログラムではもっぱら精度を追求してきましたが，数値計算の精度・速度を論ずる以前の問題として，安定的に解が得られることが必要になります．

　逐次法では，まずおおざっぱな近似解を何らかの方法で見積もり，次に近似解を逐次改良していき望みの精度が得られるまで続けます．ところが，初期値によっては真の根より遠のく場合があります．そして，このようなことが連鎖的に起こると計算は収束せずに発散してしまいます．解決策としては，初期値ｘ0を決めるときに，ｆ（ｘ0）とｆ”（ｘ0）とが同符号，すなわち，ｆ（ｘ0）・ｆ”（ｘ0）＞０になるようなｘ0を選べば収束が安定化します．

　このように，収束の安定化のためには２次までの微分を計算する必要がありますが，今回のコラムでは，ニュートン法のサブルーチンを改良し，収束を安定化を図ってみます．収束の安定化は，発散を防止し，最終収束までの所要時間を短縮しますから，結果的に収束を速くするためにも寄与してくれます．

　ニュートン法は２次収束，ベイリー法は３次収束なので，ベイリー法はニュートン法より速く収束しますが，反復回数は少ないかわりに２次導関数を計算する必要があるため，計算量が増すので結局大差ないことになります．一般的にいって，ニュートン法のほうがかえって効率的であることが多いようです．

としてΔｘを求めています．この式は根号の中が負にならないときのみ適用可能で，複号はｆ’≧０のとき－を，ｆ’＜０のとき＋を選択します．また，根号の中が負のときはΔｘ＝－ｆ／ｆ’（ニュートン法）とします．

　改良後のプログラムでは，フラグ変数ｎｅｗｔｏｎを０にすると１次の微分を用いたニュートン法，フラグを立てる（０以外にする）と２次までの微分を用いたニュートン法で計算してくれます．

2370 P20=P10*ZZZ+((DFX*(DFX-1)*Z^(DFX-2)*(1-Z)^DFY-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)+Z^DFX*DFY*(DFY-1)*(1-Z)^(DFY-2))*G1#+(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G5#+Z^DFX*(1-Z)^DFY*G6#)*ZZ*ZZ

2650 P20=P10*ZZZ+((DFX*(DFX-1)*Z^(DFX-2)*(1-Z)^DFY-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)+Z^DFX*DFY*(DFY-1)*(1-Z)^(DFY-2))*G1#+(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G5#+Z^DFX*(1-Z)^DFY*G6#)*ZZ*ZZ

　改良前のプログラムでは，χ2分布の下側５％点を計算しようとすると途中で計算が止まってしまいました．原因を調べてみると，χ2分布のパーセント点の近似値を求める段階にバグがあり，負であるはずの変数が非負になったためのエラートラップであることが判明しました．今回の修正ではその点をバグフィックスし，記述の見苦しい箇所を若干手直ししてあります．

　確率分布の計算において，これまで１次の微分を用いたニュートン法では発散するが，２次までの微分を用いたニュートン法では見事に収束するといったシーンに遭遇したことはありません．

　２次のニュートン法では，２次導関数を計算するぶん計算量が増し，どうしても余分に所要時間がかかってしまうし，１次のニュートン法でも十分安定に収束しますから，確率分布の計算にとって，２次までの微分を用いたニュートン法は必要ないのかもしれません．しかし，もしも計算時間をほとんど要さずに，２次の微分値を計算することができたならば，話は別です．

　前述したように，超幾何関数には微分・積分してもふたたび超幾何関数になるという特性があり，さらに，微分値を超幾何関数の漸化式を用いて計算すると，時間のかかる計算を回避することができます．２次までの微分を用いたニュートン法は超幾何関数向きの計算法になっているといえましょう．

1800 P20=P10*ZZZ+((DFX*(DFX-1)*Z^(DFX-2)*(1-Z)^DFY-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)+Z^DFX*DFY*(DFY-1)*(1-Z)^(DFY-2))*G1#+(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G5#+Z^DFX*(1-Z)^DFY*G6#)*ZZ*ZZ

2070 P20=P10*ZZZ+((DFX*(DFX-1)*Z^(DFX-2)*(1-Z)^DFY-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)-DFX*Z^(DFX-1)*DFY*(1-Z)^(DFY-1)+Z^DFX*DFY*(DFY-1)*(1-Z)^(DFY-2))*G1#+(DFX*Z^(DFX-1)*(1-Z)^DFY-Z^DFX*DFY*(1-Z)^(DFY-1))*G5#+Z^DFX*(1-Z)^DFY*G6#)*ZZ*ZZ

　プログラムの全文を載せると長くなるので，改良部分だけを掲載しました．ここで用いた漸化式については難しいので掲載しませんが，漸化式の計算では，「電子通信工学のための特殊関数とその応用」（奥井重彦著：森北出版）を参考にしました．

　確率分布には様々の特殊関数が用いられていますが，なかでもガンマ関数・ベータ関数などは必須項目であり，さらに非心分布ではベッセル関数・超幾何関数などの理解も必要になります．しかし，これらの特殊関数はプロの統計研究者にとっても馴染みの薄いものであり，ましてや一般のひとにとって（無論，私にとっても）とっつきにくいものです．「電子通信工学のための特殊関数とその応用」は，一松信「特殊関数入門」森北出版よりもプラクティカルな面に注力していて，現在，私の仕事に最も役立っている本となっています．特殊関数についていろいろな話題を交えながらわかりやすく解説しており，「電子通信工学のための・・・」と銘打っていますが，数理統計学における方法的基礎の副読本としてもぜひともお勧めしたい書籍です．

　ニュートン法やベイリー法では元の方程式の導関数の計算が必要になりますが，方程式の形が複雑な場合には導関数から解析的に勾配を求めることが大変になってきます．その場合には，ｘに対して微小の幅で数値差分をとり，微分式の代用とします．

　１次微分係数ｆ’（ｘ）の符号はもとの関数ｆ（ｘ）の増減と関連し，２次微分係数ｆ”（ｘ）の符号はもとの関数の凹凸（曲率）と関連していることはご存じと思いますが，数値差分（数値微分）というのは，例えば，ｙ＝ｘ^2を微分してｙ’＝２ｘのように解析的に求めるのではなく，ｘのある点における関数ｙ＝ｆ（ｘ）の傾きを差分近似を用いて数値的に求める数値解法をいいます．

で差分近似して計算するほうが取り扱いが容易になります．いずれも導関数を差分商で置き換えるのですが，（１）を前進差分近似，（２）を後退差分近似，（３）を中心差分近似と呼びます．

　幅ｈが０に近づくとき割線は接線に近づきますから，差分を使うとき，ｈは原理的には小さいほど理想的ですが，あまりにも小さいと差分の中に占める数値計算誤差の比重が大きくなり，かえって精度が悪くなってしまいます．関数ｆの性質とｈの取り方によって精度は決まりますが，誤差について次のように検討してみましょう．

より，１次微分係数は近似的に前進差分（１）として求められます．この公式の誤差は｛１／２ｆ”（ｘ）ｈ＋１／６ｆ'''（ｘ）ｈ^2＋・・・｝の部分ですから，ｈにほぼ比例すると考えることができるため，｛・・・｝括弧の部分をＯ（ｈ）と表わします．

となり，中心差分近似の打ち切り誤差はｈ^2に比例します．このことから，中心差分は前進差分・後退差分より平均的な打ち切り誤差が少なく，精度が高いことが理解されます．ただし，関数値の計算回数は余計にかかります．

　次に，２次微分係数ですが，前進差分近似，後退差分近似，中心差分近似のそれぞれに対応した数値微分を作ることが可能です．例えば，中心差分近似式（３）を用いると

　一階中心差分近似式，二階中心差分近似式にはそれぞれに係数［１，－１］，［１，－２，１］がはいっていましたが，それ以上の高階差分近似式にも，［１，－３，３，－１］，［１，－４，６，－４，１］など二項係数 nＣk が符号を交代させた形で入っていることが示されます．これは数学的帰納法を用いて簡単に証明できます．一般式の形で書くこともできますが，読者の演習問題にしておきます．

　また，中心差分近似式ｆ”（ｘ）＝｛ｆ（ｘ＋ｈ）－２ｆ（ｘ）＋ｆ（ｘ－ｈ）｝／ｈ^2を採用した場合の３次微分係数と４次微分係数の差分近似式は，

　　　　　　　　＝１／４ｈ^4｛ｆ（ｘ＋３ｈ）－２ｆ（ｘ＋２ｈ）－ｆ（ｘ＋ｈ）＋４ｆ（ｘ）－ｆ（ｘ－ｈ）－２ｆ（ｘ－２ｈ）＋ｆ（ｘ＋３ｈ）｝＋Ｏ（ｈ^2）

　なお，導関数を近似する差分式で，必要とする精度をもつ公式は原理的には無数に導くことができます．たとえば，Ｏ（ｈ^2）の精度をもつ３点差分公式

　　ｆ’（ｘ）＝１／６０ｈ｛ｆ（ｘ＋３ｈ）－９ｆ（ｘ＋２ｈ）＋４５ｆ（ｘ＋ｈ）－４５ｆ（ｘ－ｈ）＋９ｆ（ｘ－２ｈ）－ｆ（ｘ－３ｈ）｝＋Ｏ（ｈ^6）

　導関数の計算をコンピュ－タで行う場合，連続な変化は扱えないので，小さな変化を決めて離散的に解く数値解法が主体になります．差分近似式（数値微分式）では元の方程式の形が複雑であっても確実に勾配の求まる利点があり，コンピュータを用いた数値解析では，解析的な微分式よりも差分近似式が多用されます．そのため，差分法（ＦＤＭ）は微分方程式の数値解法，有限要素法（ＦＥＭ），境界要素法（ＢＥＭ）にとって欠かせない知識になっています．

　コンピュータ計算では多くの場合，微分方程式や積分を差分解法によって解きますが，この手法は詰まるところ，アナログのデジタル変換による解法といい換えることができます．そのため，微分方程式を正しく近似している差分方程式を使ったとしても，その解は微分方程式の解の近似になるとは限らないという矛盾を生じ，微分方程式の解とはまったく違う性質がでてくることがあります．したがって，無数にある差分方程式から安定なものを選び出し，元の微分方程式の解の性質を保存するように差分化することが必要になってきます．

　２分法や割線法では最初に解をはさむ２つの近似値が必要です．ｆ（ｘ0）＝０となる真の値ｘ0をはさむような近傍の値ｘ1，ｘ2を入力すると，ｆ（ｘ1）とｆ（ｘ2）の符号が異なるので，根ｘ0は区間（ｘ1，ｘ2）の中に必ず存在します．２分法では中点ｘ3＝（ｘ1＋ｘ2）／２，割線法ではｘ軸との交点ｘ3＝｛ｘ1ｆ（ｘ2）－ｘ2ｆ（ｘ1）｝／｛ｆ（ｘ2）－ｆ（ｘ1）｝でのｆ（ｘ3）の符号を調べてみて，もしｆ（ｘ1）とｆ（ｘ3）が同符号ならばｘ3を新しいｘ1として，そうでなければｘ3を新しいｘ2として計算を繰り返します．

　逐次２分法は１次収束，割線法は１．６１８次収束なのでニュートン法の２次収束，ベイリー法の３次収束に比べ収束の遅い漸近近似法ですが，これを続けていけば根ｃの存在する区間は毎回狭められ絞り込まれていき，必ず真の値に収束することになります．導関数の計算が困難であったり，時間がかかるならば２分法や割線法を使ったほうがよいでしょう．

　以上は１変数の非線形方程式ｆ（ｘ）＝０の解法でしたが，多変数の場合でも考え方は同じで，前節のニュートン・ラプソン法Δｘ＝－ｆ（ｘ）／ｆ’（ｘ）は非線形連立方程式の解を求める場合にも拡張することができます．

で与えられますが，ベクトルを用いるとｆ（ｘ）＝０のように書けますから，多変数関数においてはΔｘとｆ（ｘ）はベクトルに，ｆ’（ｘ）はヤコビ行列と呼ばれる行列になるだけの違いです．

の２つの変数ｘ，ｙに関して近似解（ｘ0，ｙ0）がわかっている場合，その解の精度を上げるニュートン・ラプソン法について紹介します．より改良された近似解をｘ＝ｘ0＋Δｘ，ｙ＝ｙ0＋Δｙとします．この連立方程式ではパラメータｘ，ｙに関して非線形ですので，パラメータの近似値ｘ0，ｙ0の回りで１次までのテイラー展開を行ない線形化します．

　ここで，ｆ（ｘ，ｙ）はｆ（ｘ0，ｙ0）に比べて，ｇ（ｘ，ｙ）はｇ（ｘ0，ｙ0）に比べて小さいと考えることができますから，これを無視して整理するとΔｘ，Δｙに関しての線形方程式

と書くことができます．ヤコビアンはヤコビの名をとどめる行列式です．韻をふんだ名前ではガリレオ・ガリレイが有名ですが，ヤコビ・ヤコブというのもその類です．ヤコビが先鞭をつけた関数行列式はヘッセなどに引き継がれ，解析幾何学の面でたびたび利用され発展しました．ヘッセにもヘシアンという彼の名をとどめる行列式があり，２変数の場合は，