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

　統計解析では，特定の確率分布において，確率変数がある値以上または以下の値をとる確率（片側確率）と逆に片側確率を与えたときに対応する確率変数の値（パーセント点）を計算する必要が生じます．

　一昔前までは，これらの値を求めるのに統計数値表をひくことが必要になり，検定や推定は大変手間のかかる作業でしたし，数表が利用できないときはまったくお手上げ状態でした．現在では，統計プログラムがこの作業を肩代わりしてくれるので，検定・推定は大変身近なものになっています．

　ところで，統計プログラムでは，統計数値表をパソコンに記憶させておいて，片側確率やパーセント点を求めているわけではありません．統計数値表をパソコンに記憶させておいて使う方法は非常に多くの数値をデータとして格納しなければならず，メモリを占有するだけでなくプログラムが大きくなりすぎて実用的ではないからです．

　正規分布・ｔ分布・Ｆ分布・χ2分布では，片側確率やそのパーセント点を計算する簡便さと適当な精度をもつ近似関数が Hastings, Wallace, Paulson, Wilson-Hilferty などにより求められています．これらの近似関数を利用して計算に必要な数値はすべてプログラム内部で求めるようにできればわざわざ数表を引く必要がなくなり大変便利ですし，手間もかからず効率的に片側確率とパーセント点の値を得ることができます．近似式とはいっても，相対誤差は％オーダーですから，実用的には十分な有効数字があります．実際，これらの近似式は有意性の判定等に利用されていて，多くの場合，近似式で所期の目的を十分に達成することができます．

　ところが，近似式を使って必要な数値を計算する方法には欠点があります．ｔ分布・Ｆ分布・χ2分布では，任意の自由度の全範囲にわたってよい精度を保証してくれる近似式はないのです．多くの統計計算パッケージでは，パーセント点の計算に関して，近似式を採用していますが，一般的にいって，近似式では，自由度νの小さいところでｔ（ν，α），自由度ν2の小さいところでＦ（ν1，ν2，α）の誤差がかなり大きいのですが，たとえば，後述する Wallace のアルゴリズムは４よりも小さい自由度に関しては信頼できません．統計解析ではあまり使わないところですから，実用上ほとんど問題ないといって妥協したくはありません．

　一方，片側確率の計算に関しては，部分積分によって求められた展開式に基づいて正確な値が得られるアルゴリズムを，多くの統計計算パッケージが採用しています．厳密な級数を足し上げる方法ですから，誤差は演算誤差だけになります．それらは，計算精度のみならず計算速度の点でも満足のゆくものですが，今度は逆に整数自由度でないと使えないという欠点が生じます．

　ウェルチのｔ検定や広津の累積χ2検定等では小数付き自由度での片側確率・パーセント点を求める必要がありますが，その場合，整数自由度で求めた２つの片側確率・パーセント点より，補間で求めることが必要になり，結構面倒ですし，なによりせっかくの精度が台なしです．そのため，小生は，精度よりも小数付き自由度でも使えるという条件を優先させて，これまで近似式のほうを好んで用いてきました．

　小数付き自由度の有意点について，自由度が大きければ Hastings, Wallace, Paulson, Wilson-Hilferty などの近似式で求めればよいのですが，自由度が小さい場合，正確な値が得られません．そこで，今回のコラムでは，整数・小数任意の自由度で，かつ，正確な値を保証してくれる超幾何関数を用いた確率分布の計算プログラムを取り上げてみることにします．

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

　そこで，まず，近似式を使って，片側確率からそのパーセント点を求めるプログラムとその逆関数，すなわち，パーセント点を与えて片側確率を求めるプログラムを紹介します．以下では，下側確率をP0，上側確率をPP，正規分布，ｔ分布，Ｆ分布，χ2分布の各パーセント点をそれぞれUUN,UUT,UUF,UUX，自由度をDFで表しています．

　正規分布の確率密度関数の指数部分ををテイラー展開して積分すると片側確率の漸近展開式が得られますが，収束が遅く，計算式としてはあまり実用的ではありません．その点，ラプラスやシェントンによる連分数展開式は収束が速くかつ精度も高い式といえましょう．

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

　χ2分布の密度関数はガンマ関数を，分布関数は不完全ガンマ関数を含んでいます．ガンマ関数は階乗に密接な関係があり，引数が整数や半整数ならば，近似解ではなくて，解析的な解を得ることができますが，任意の自由度では，部分積分による正確なアルゴリズムは使えません．

　自由度νをもったχ2について，Wilson-Hilferty は（χ2/ν）^(1/3)が，平均値(1-2/9ν)，分散2/9νの正規分布で近似できることを示しました．この近似はＦ分布についての近似値を導くのにも応用されますが，ここでは，Wilson-Hilferty 近似に基づくアルゴリズムを紹介します．

　ｔ分布の場合も，正確な片側確率を求めるアルゴリズムは部分積分から作成できるのですが，小数付き自由度では使えないという欠点があります．そこで，Wallace は任意自由度で使える次の変換と逆変換を提案しました．このように，厳密に正規分布ではないもののそれに近い分布に対しては正規分布を用いて確率の近似計算を行なうことができます．

　ｔ分布やＦ分布は不完全ベータ関数と密接に関係しています．不完全ベータ関数の部分積分に基づくＦ分布の正確な確率計算アルゴリズムは，厳密な級数を足し上げる方法ですが，小数付き自由度では使えませんから，χ2 分布の Wilson-Hilferty 近似に基づいた Paulson の結果を利用します．

　５１３０行では，ν2が３以下の場合，精度をよくするための補正を行っています．この補正によって，自由度の全範囲にわたって数％の相対誤差という精度で計算してくれます．

　以上，小数付き自由度の場合でも使える確率分布の近似計算プログラムを掲載してきました．各々のサブルーチンは非常にコンパクトにまとまっています．それに対して，整数自由度の場合の厳密な級数和の計算プログラムは，自由度が偶数の場合と奇数の場合に分けて計算する必要がでてきますから，プログラムは意外に長くなります．補間プログラムまで含めるとさらに長くなるので，割愛することにしました．要するに，整数自由度の場合にしか使えない片側確率の計算プログラムはダサイと感じられるのです．

　これらの近似式では，自由度νの小さいところでｔ（ν，α），自由度ν2の小さいところでＦ（ν1，ν2，α）の誤差がかなり大きいのですが，統計解析ではあまり使わないところであるから，実用上ほとんど問題ないといって済ませるわけにはゆきません．たとえば，広津の累積χ2検定では小さい自由度でかつ小数付き自由度に対するχ2値が必要になるからです．

　　Ｆ（α，β，γ：ｘ）＝１＋αβ／γｘ＋1/2!α（α＋１）β（β＋１）／γ（γ＋１）ｘ^2＋1/3!α（α＋１）（α＋２）β（β＋１）（β＋２）／γ（γ＋１）（γ＋２）ｘ^3＋・・・

であるような級数にほかなりません．指数関数，対数関数，三角関数，２項関数，ベッセル関数，直交多項式列，不完全ガンマ関数，指数積分，ガウスの誤差関数なども超幾何級数であって，超幾何関数は一般に収束半径１をもちます．

　正規分布は誤差関数と，χ2分布は不完全ガンマ関数，ｔ分布やＦ分布は不完全ベータ関数と密接に関係していて，正規分布，χ2分布，ベータ分布の分布関数はそれぞれ以下のように表現されます．

　このことさえ知っておけば，正規分布・χ2 分布・ｔ分布・Ｆ分布の片側確率が計算できます．以下に掲げるプログラムは，前述の近似式よりは長くなりますが，小数自由度にも対応していて，しかも，整数自由度の場合の計算プログラムよりも短く済んでいます．

　１６９０行～２０１０行が超幾何関数の級数和を求めるサブルーチンです．PC-9801の場合，単精度計算は１０進法で７桁の精度で格納され，６桁までの桁数で印刷されます（倍精度計算は格納・表示とも１６桁）．そこでプログラムでは，要求する精度が1E-9以下になったら反復を終了させるようにしています（１７３０行・１８９０行）．このように，超幾何関数を用いれば，確率計算を任意の精度で行うことができるという精度管理上のメリットが得られます．また，必要な精度があらかじめ定まっている場合には，その精度に見合った値を選んで使うことになります．

　このプログラムの１２１０行ではガンマ関数値，１３７０行～１３９０行，１５８０行～１６００行ではベータ関数値の計算をしています．２０２０行から２２１０行がガンマ関数のサブルーチンです．

したがって，ｘが正の整数ｎのときには，Γ（ｎ＋１）＝ｎ！が成り立ち，ガンマ関数は階乗の一般形となっていることがわかります．階乗の解析的補間をしている関数がガンマ関数なのです．

です．引数の値が整数（ｎ＋１）でも半整数（ｎ＋１／２）でもないガンマ関数に対する正確な計算は複雑ですが，物理・化学・工学などの諸分野の問題の解を求めるときには，どうしてもそれらの解を数値を用いて表す必要がでてきます．現在では膨大な数表を用いることは稀であって，コンピュータが広く使用されていますから，必要に応じて計算する方法がとられています．

を用いています．この近似式は有効範囲が０≦ｘ≦１に限られていますが，ガンマ関数にはΓ（ｘ＋１）＝ｘΓ（ｘ）の関係があり，この漸化式を繰り返し適用して，０≦ｘ≦１の範囲になるように再帰的な方法を用いればｘ≧１にも拡張することができます．また，ガンマ関数Γ（ｘ）はｘが大きくなるとすぐにオーバーフローエラーを起こしてしまいますので，プログラムではｌｎ（Γ（ｘ））＝ｌｎ（Γ（ｘ＋１））－ｌｎ（ｘ）によってその対数ｌｎ（Γ（ｘ））を求めて飽和現象を回避しています．プログラムの２１８０行目では，いわゆるホーナーの方法によって８次多項式の値は求めています．ホーナーの方法は計算回数を少なくし，計算誤差の点からも有効な方法です．なお，ガンマ関数Γ（ｘ）はｘ＞０について微分可能で，ｘ＝１．４６１６３２１４４９・・・で最小となります．

　ガンマ関数は，Γ（ｎ＋１）＝ｎ！より，階乗の補間関数であり，初等的でない関数の中で最も簡単かつ重要な数学的関数といえましたが，ベータ関数とガンマ関数との間には

　１１００行～１１３０行，１２４０行～１２７０行，１４４０行～１４７０行，１６４０行～１６７０行では，超幾何関数の微分を行っていますが，超幾何関数には微分・積分してもふたたび超幾何関数になるという特性があります．

　以下に掲げるパーセント点を求めるアルゴリズムでは，この性質をうまく利用しています．近似式で求まるパーセント点を初期値として，これを３３２０行～３６１０行の１次の微分まで用いたニュートン法によって補正すると，精度の良いパーセント点が得られます．

　ニュートン法の幾何学的な意味は「初期値ｘ0における関数の勾配を求めて接線とｘ軸との交点を求める．この点において同様の作業を反復させて行なうと，ｘは順次根に近づいてゆく．」と解釈されます．そして，ｘの変化があらかじめ与えた許容誤差以内になったとき，そのｘをもって根と判定します．

　超幾何関数の利点を活かした確率分布の計算法を紹介しました．このアルゴリズムから導かれた値と真値を比較した結果は示しませんでしたが，単精度の範囲で数値表の値とよく一致し，非常に精確な値が得られました．小数付き自由度を含む任意の自由度についての片側確率とパーセント点が簡便に求められるうえに，精度が高く，その実用上の意味は大きいと考えられました．

　さて，超幾何関数の値を求める場合，級数を収束するまで加えればよいのですが，引数が負のときは交代級数ですから，桁落ちが起こるという問題があります．

　また，せっかく超幾何関数の高精度計算をしたわけですから，次の課題として，ガンマ関数を高精度で計算するサブルーチンを用意したいところです．ガンマ関数Γ（ｘ）の値を求めるには，いろいろな近似式，例えば，ヘイスティングス，コリンジ，ハートらによる多項式近似やスターリングの漸近展開式などがあります．

およびその逐次導関数φ’（ｘ），φ”（ｘ），・・・，φ(k)（ｘ）（トリガンマ関数，テトラガンマ関数，ペンタガンマ関数など，総称してポリガンマ関数と呼ぶ）の値を求めるときにも応用可能で，有理式であるハートの近似式と比べても大変便利な近似式になっています．

　ヘイスティングスの近似式がΓ（ｘ）の１≦ｘ≦２の範囲の近似式であったのに対し，コリンジの近似式はΓ（ｘ）の２≦ｘ≦３の範囲の近似式です．ヘイスティングスの近似式と同様の再帰的な方法で，ｘ＜２，ｘ＞３にも拡張できます．

　１／Γ（ｘ）を区間１．５≦ｘ≦２．５内でチェビシェフ補間して得られた近似式で，引数を倍精度にすればガンマ関数を倍精度で計算してくれるサブルーチンです．

　このサブルーチンでは４１５０行で連分数展開を用いています．スターリングの近似式はｘが大きいほど精確なので，ｘがある値（変数名ｎｚ）より小さいときは漸化式

　超幾何関数が統計計算に現れることは少なく，多くのひとにとって超幾何関数は馴染みの薄い関数であることは間違いありません．ここでは，小生がこのアイディアをどのようにして得たのかについて，簡単に紹介したいと思います．

　このように，超幾何関数は非心χ分布の積率と関係しています．非心χ分布の積率計算では必然的に超幾何関数を用いることになりますが，そのような中から超幾何関数による確率計算を着想したというわけです．

　なお，χ分布を拡張する方向としては，ひとつには自由度を増すこと，もうひとつは非心パラメータをつけることが考えられます．χ分布の自由度を増すと，半正規分布→レイリー分布→マクスウェル分布→・・・となりますが，非心χ分布の自由度を増すと，折り重ね正規分布→ライス分布→・・・が得られます．すなわち，非心χ分布において，自由度１の場合が折り重ねられた正規分布（folded normal distribution）で，期待値が０でない正規分布をｙ軸（ｘ=０）で折り返して得られた分布になっています．また，自由度２の非心χ分布はライス分布です．

　通信では電波が互いに干渉しあって，うねりのような強弱がついて非常に聴取しづらい現象が起こります．受信点では様々な伝搬経路を通ってきた多数の電波が合成されるからであり，この干渉をフェージング（fading）といいます．地表波と空間波の干渉を近距離フェージング，また，経路の異なる空間波どうしの干渉によるものを遠距離フェージングといいます．

　フェージングはその強度が等しいとき最も激しく起こりますが，合成された電波の位相ならびに振幅が不規則に変化するので，そのデータはランダムな性格を有し，明確な数式で記述されるよりは確率的記述と統計的平均値で表されなければなりません．その際必要とされるのがライス分布です．ライス分布はマーカムのＱ関数と関連していて，電気通信分野ではｎ分布の名で通っています．ライス分布において，μ＝０ならばレイリー分布に一致します．また，μ→∞のときの極限は正規分布になります．

　ライス分布は第２次大戦中，米国ベル研のライスと本邦の仲上稔教授によって独立に研究されたものですが，とくに，仲上氏のフェージングの研究は戦時中のおそらく実験設備も十分でない環境にもかかわらず，先進国米国の研究に先立つオリジナルな成果として発表されています．しかし，originally made in Japan であるにもかかわらず，本邦ではその重要性が認められませんでした．ライス分布と名づけられていますが，この機会に仲上・ライス分布として紹介したいと考えています．

　非心分布は，元来，検定の検出力を調べるために導入されたものですが，ここまで読んで頂いた方への特別付録として，非心ｔ・Ｆ・χ2分布の片側確率とパーセント点を近似計算するプログラムを掲載します．