DelphiでのExcel統計関数：NORM、CHISQ、BETA

セルに=NORM.DIST(115,100,15,TRUE)と入力すると、Excelは何の苦もなく0.8413447を返します。この呼び出しは単なるルックアップのように見えますが、そうではありません。その1つの数値の背後には、閉形式（closed form）を持たない積分である累積正規分布が存在し、CHISQ.INV.RTやBETA.DISTの背後には、安易に近似するのではなく慎重なライブラリが正しく評価しなければならない特殊関数が存在します。Excel互換を謳うスプレッドシートコンポーネントは、Excelが表示する最後の桁までこれらの値を再現する必要があり、そのためには関数名だけでなく数値計算手法そのものを再現しなければなりません。

HotXLSは50以上のこれらの統計関数を実装しており、それらを正確にするための処理は数式バーからはほぼ完全に不可視です。これは、エンジンがこれらをどのように計算しているか（共通の特殊関数コア、算術演算の安定性を保つための分岐決定、そして一般的なケースでは通過しないために長い間末尾に隠れていた1つの逆正規分布のバグ）を示すツアーです。

1つのワークシート呼び出し、その背後にある50の分布

関数は、統計ワークブックが使用するファミリーを網羅しています。正規分布ファミリー（NORM.DIST、NORM.S.DISTとそれぞれの逆関数）、ガンマおよびカイ二乗分布ファミリー（GAMMA.DIST、CHISQ.DIST、CHISQ.DIST.RT、CHISQ.INV.RT）、ベータ分布ファミリー（BETA.DIST、BETA.INV）、標本分布（T.DIST、T.DIST.2T、F.DIST、F.INV）、離散分布のペア（BINOM.DIST、POISSON.DIST）、および推測統計のヘルパー（CONFIDENCE.T、CONFIDENCE.NORM）があります。呼び出し側から見れば、それぞれが単一の数式です。セルに入力を設定し、ワークブックに評価を求め、結果を読み取ります。

var
  wb: IXLSWorkbook;
  sh: IXLSWorksheet;
begin
  wb := TXLSWorkbook.Create;
  sh := wb.Sheets.Add;
  sh.Range['A1', 'A1'].Value := 115;   // observation
  sh.Range['A2', 'A2'].Value := 100;   // mean
  sh.Range['A3', 'A3'].Value := 15;    // standard deviation

  // The XLS formula parser uses ';' as the argument separator.
  Writeln(wb.Calculate('=NORM.DIST(A1;A2;A3;TRUE())'));   // 0.8413447
  Writeln(wb.Calculate('=CHISQ.INV.RT(0.05;10)'));        // 18.3070381
  Writeln(wb.Calculate('=BETA.DIST(0.5;2;3;TRUE())'));    // 0.6875
end;

ワークブックのCalculateメソッドは、稼働中のシートに対して臨時の数式をコンパイルして評価し、Variantを返します。最初に試す際に注意すべき点があります。Calculateの背後にある数式パーサーは、引数の区切り文字としてセミコロンを使用するため、=SUM(A1,B1)ではなく=SUM(A1;B1)になります。保存されるセル数式はExcel標準のカンマを維持します。同じ評価器が以下のすべての統計関数をディスパッチするため、1つがCalculateで動作すれば、残りの関数も同じパスに従います。

他のすべての基盤となる2つの関数

このセットにおけるほとんどの累積分布は、それ自体の定義を合計したり積分したりすることで計算されているわけではありません。これらは2つの特殊関数から計算されます。P(a, x)と書かれる正規化下方不完全ガンマ関数（regularized lower incomplete gamma function）と、I_x(a, b)と書かれる正規化不完全ベータ関数（regularized incomplete beta function）です。内部的に、これらはディスパッチャーが依存するヘルパーであり、その連鎖は短いです。カイ二乗CDFは、形状df/2およびスケール2のガンマCDFです。ガンマCDFはP(a, x)そのものです。t分布、F分布、および二項累積関数はすべて、適切な引数における正規化不完全ベータの値です。ポアソン分布のCDFは、上方不完全ガンマ関数Qです。ガンマ関数とベータ関数を適切に実装すれば、多くの分布がその精度を無償で継承します。

「正規化」という言葉がすべての鍵です。生の不完全ガンマは階乗のように増大し、生のベータ積分は結果が確定するはるか前にアンダーフローまたはオーバーフローを起こす可能性があります。正規化された形式は、完全ガンマ関数または完全ベータ関数で除算されているため、完全に0から1の範囲に収まります。これは確率が占める範囲そのものです。この正規化により、中間項がdoubleの限界値を超えて暴走することなく、自由度2のカイ二乗と自由度200のカイ二乗に対して同じルーチンを使用することができます。また、密度項の長い末尾を加算することによってCDFを計算しない理由もこれによって説明されます。各項が独自の丸め誤差を伴い、シリーズが進むにつれて誤差が累積するためです。正規化された特殊関数は、代わりに急速に収束するシリーズまたは連鎖分数を評価することで、その加算処理を完全に回避します。

対角線の下側はシリーズ、上側は連鎖分数

不完全ガンマ関数のルーチンは、計算を行う前に1つの決定を下します。xをa + 1と比較します。この境界は任意ではありません。P(a, x)のべき級数展開（power-series expansion）は、xがaに対して小さい場合に急速に収束し、xが大きい場合は収束が遅く、最終的には使い物にならなくなります。連鎖分数（continued fraction）は逆の性質を持っています。そのため、エンジンはxがa + 1未満の場合はべき級数を使用し、xがa + 1以上の場合はLentzの連鎖分数を使用し、各分岐はその分岐が得意とする処理のみを担当します。

連鎖分数には1つのガードが必要です。Lentzの方法は、実行中の分子と分数を保持し、各ステップで分数を反転させることで機能します。もし一方がゼロに近づくと、反転処理が破綻します。解決策は極小のフロア値です。中間項の絶対値が約1e-30を下回る場合は常に1e-30にクランプし、収束値を乱すことなく再帰を有限に保ちます。同じクランプ処理が、同じ理由で不完全ベータ関数の連鎖分数にも現れます。これは安定した評価と、ゼロと区別がつかない何かによる除算との分岐点となる、重要な負荷を支える小さな定数です。

上方の尾部であるQ(a, x)は単に1マイナスP(a, x)であり、これがポアソンの累積分岐の計算方法です。平均λで最大k回のイベントが発生する確率はQ(k + 1, λ)です。k + 1個のポアソン項を加算するのではなく、上方不完全ガンマを介してルーティングすることは、これもまた多数の小さな値を累積する代わりに、1つの収束式を評価するという選択です。

階乗のオーバーフローを起こさない離散質量

離散分布は異なる危険性をはらんでいます。二項確率質量（binomial probability mass）は二項係数を伴い、52選抜26といった係数は巨大な整数になります。これを直接形成すると、妥当な確率に戻すための除算を行う前に、分子がdoubleをオーバーフローします。エンジンはこれを直接形成しません。log-gamma関数を介して対数（log）空間で階乗を計算し、対数を加減算し、成功と失敗の確率の対数を組み込み、最後に一回だけ指数関数（Exp）を適用します。

// Binomial probability mass, evaluated entirely in log space.
// LnGammaF(n+1) is ln(n!); the three log-factorials form ln(C(n,k)),
// and the whole exponent is built before a single Exp call.
//   ln P(X=k) = ln(n!) - ln(k!) - ln((n-k)!) + k*ln(p) + (n-k)*ln(1-p)
result := Exp(LnGammaF(nt + 1) - LnGammaF(kk + 1) - LnGammaF(nt - kk + 1)
  + kk * Ln(pp) + (nt - kk) * Ln(1 - pp));

log-gamma関数自体はLanczos近似であり、正の軸全体で正確であり、評価コストも低いです。すべての大きな値は、最後のExpまでその対数として保持されるため、ルーチンが具現化する最大の数値は確率そのものであり、それは最大でも1です。ポアソンの質量関数も同じレシピに従い、分母の階乗を表す単一のlog-gamma項を備えています。エッジケース（pが正確にゼロまたは1である場所）では閉形式が特別扱いされるため、コードがLn(0)を呼び出すことはありません。HotXLSは、BINOM.DIST(5,10,0.5,FALSE)に対しては0.2460938を返し、累積のPOISSON.DIST(2,2,TRUE)に対しては0.6766764を返し、Excelが出力する桁数と一致します。

順方向曲線のブラケッティングによる逆関数

逆分布関数（inverse distribution function）は逆の質問をします。確率が与えられた場合、そのCDFがそれに等しくなる値を求めます。このセットの中で、高速な直接式を持つ逆関数は1つだけです。標準正規逆関数であるNORM.S.INVは、Acklamの有理式近似（rational approximation）を使用します。これは、中央領域と2つの尾部に分割され、範囲全体でほぼdouble精度に正確な一対の多項式比率です。反復を伴わない閉形式の評価です。

他の逆関数にはそのような数式がないため、エンジンはそれらを数値的に逆算します。分布のサポート範囲から選択した上限と下限で答えを挟み込み（ブラケッティング）、二分（バイセクト）します。中点で順方向CDFを評価し、ターゲット確率を内包し続けるようにいずれかの境界を移動させ、区間が狭くなるまで繰り返します。ガンマおよびカイ二乗の逆関数の場合、ブラケットはゼロと、形状およびスケールから構築された余裕を持った上限見積もりから開始し、確率がまだ内包されていない場合は上限を2倍にします。t分布の逆関数は、外側に広がる対称的な境界を挟み込みます。F分布の逆関数は非負の区間で二分します。コストは呼び出しあたり数十回のCDF評価ですが、これはスプレッドシートの処理速度では不可視であり、メリットはすべての逆関数がそれが反転させる順方向関数とまったく同じ精度になることです。そのため、CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE)のような往復処理は、誤差の範囲内で0.7を返します。

尾部に隠れていた常用対数のバグ

ここに、長い間残っていたバグの話があります。Acklamの逆正規分布ルーチンには3つの分岐があります。確率が約0.025と0.975の間にある場合に使用される広い中央分岐は、どこにも対数を使用しない多項式比率を実行します。非常に小さい、または非常に大きい確率のための2つの尾部分岐は、まず入力の対数を取ります。尾部はpのマイナスの自然対数の平方根のように動作するためです。

尾部分岐の初期バージョンは、自然対数（natural log）を使用すべき場所に常用対数（base-10 log）を使用していました。この2つは約2.30の定数係数で異なるため、尾部の結果は一貫して大きく間違っていました。それにもかかわらず、カジュアルなチェックでは問題がないように見えました。カジュアルなチェックは中央部分で行われるためです。NORM.S.INV(0.5)はゼロであり、NORM.S.INV(0.975)は教科書通りの1.959964であり、これら両方は対数を一切呼び出さない中央の多項式を通過します。エラーが現れたのは、確率が尾部に入り込んだ場合（例えばNORM.S.INV(0.001)。これは-3.0902323を返さなければならないところ、自然対数対常用対数の比率のズレに近い値が返されていました）だけでした。逆正規分布に依存するすべての関数（信頼区間ヘルパーなど）も、同じ歪みを継承していました。教訓は平凡で高価なものです。分岐構造を持つ関数にはすべての分岐内にテストポイントが必要であるということです。正しい共通パスは、破損した稀なパスを平然と覆い隠してしまうためです。修正は常用対数から自然対数への1トークンの変更であり、尾部の値はExcelの値にスナップしました。

xの符号がt分布の尾部を決定する

Studentのt累積関数は、逆になりやすい微妙な点を備えています。その値は、df / (df + x²) において評価される正規化不完全ベータ関数から得られますが、そのベータ値はxの絶対値を超える尾部の確率を測定したものであり、xまでの累積確率ではありません。t分布の対称的な形状は、変換がxがゼロのどちら側に落ちるかに依存することを意味します。

// Student t CDF. ib is the regularized incomplete beta at df/(df+x*x),
// which measures the symmetric tail. The cumulative value depends on
// the sign of x; returning ib unconverted gives the wrong tail.
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
  result := 1 - 0.5 * ib        // above the mean: one minus half the tail
else if x < 0 then
  result := 0.5 * ib            // below the mean: half the tail
else
  result := 0.5;                // exactly at the mean

xがゼロを超える場合、累積確率は1マイナス対称尾部の半分です。xがゼロ未満の場合、対称尾部の半分です。ゼロの場合はちょうど半分です。ベータ値を直接返してしまうと、任意のゼロ以外のxに対して曲線全体がずれ、間違った側の確率を報告することになります。右尾部（one-tailed right）および両尾部（two-tailed）のバリアントも同じ分岐に基づいて構築されているため、T.DIST.2T(1,1)は0.5を返し、T.DIST(1,1,TRUE)は0.75を返し、逆関数のT.INVはこの補正されたCDFに対して二分処理を行うため、往復が閉じます。

セルからはこれらの詳細は一切見えず、それが意図された結果です。数式を書き、Excelと一致する数値を読み取ります。独自のロジックでエンジンを拡張する場合、関数の登録メカニズムは数式エンジンとカスタム関数に関するチュートリアルでカバーしています。また、数式がシートをまたいで名前付き範囲にアクセスする方法は、定義された名前とシートをまたぐ数式に関する記事でカバーしています。これらすべてが、読み込み、書き込み、グラフ作成、およびフォーマットAPIと並んで、DelphiおよびC++Builder向けのHotXLS spreadsheet componentに含まれています。