セルに=NORM.DIST(115,100,15,TRUE)と入力すると、Excelは即座に0.8413447を返します。この呼び出しは単なる表参照のように見えますが、実際はそうではありません。その1つの数値の背後には、閉形式を持たない積分である累積正規分布が存在します。また、CHISQ.INV.RTやBETA.DISTの背後には、適当な近似値で済ませるのではなく、ライブラリが厳密に評価しなければならない特殊関数が存在します。Excel互換を謳うスプレッドシートコンポーネントは、Excelが表示する最後の1桁までこれらの値を再現する必要があり、それは関数の名前だけでなく、背後にある数値計算手法そのものを再現することを意味します
HotXLSはこれら50以上の統計関数を実装しており、それらを正確に機能させるための内部処理は数式バーからは完全に見えません。この記事では、エンジンがこれらを計算する仕組み（共通の特殊関数コア、計算の安定性を維持するための分岐処理、および一般的なケースでは問題のコード行が実行されなかったために長い間検出されなかった正規分布の逆関数に関するバグ）について紹介します
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などがあります。呼び出し側から見れば、それぞれが単一の数式です。セルに入力を設定し、ワークブックに評価を依頼し、結果を読み取るだけです
これらの数式がすでに含まれているワークブックを計算するには、ワークブックを開き、Recalculateを1回呼び出して、結果が出力されたセルを読み取ります。TXLSXWorkbook.Recalculateは依存関係グラフを走査し、トポロジカル順序で各ダーティ数式を評価するため、数式は1回だけ計算され、その結果がセルに直接書き込まれます
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;ワークブックがエンジンに搭載されていない関数を必要とする場合、同じ評価エンジンを使用して、カスタムワークシート関数のフックを介して独自のロジックをプラグインできます。カスタム関数は値を自由に返せるため、組み込み関数と全く同じように動作します
すべての基礎となる2つの関数
このセットに含まれるほとんどの累積分布は、その定義式を足し合わせたり積分したりして計算されるわけではありません。これらは2つの特殊関数、すなわち「正則化下側不完全ガンマ関数」P(a, x)と「正則化不完全ベータ関数」Ix(a, b)から計算されます。内部的には、これらがディスパッチャが依存するヘルパー関数であり、その処理チェーンは短いものです。カイ二乗累積分布関数（CDF）は、形状df/2およびスケール2のガンマ累積分布関数です。ガンマ累積分布関数はP(a, x)そのものです。t分布、F分布、および二項分布の累積分布関数はすべて、適切な引数を指定した正則化不完全ベータ関数の値です。ポアソン分布の累積分布関数は上側不完全ガンマ関数Qです。ガンマ関数とベータ関数を正しく実装するだけで、多くの分布がその精度をそのまま受け継ぐことができます
対角線より下は級数展開、上は連分数展開
不完全ガンマ関数のルーチンは、計算を行う前に1つの判断を下します。xをa + 1と比較するのです。この境界線は恣意的なものではありません。P(a, x)のベキ級数展開は、aに対してxが小さい場合は速やかに収束しますが、xが大きい場合は収束が遅くなり、最終的には使い物にならなくなります。連分数展開はこれとは逆の特性を持っています。そのため、エンジンは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つの収束する数式を評価する方が優れているためです
階乗のオーバーフローを防ぐ離散確率の計算
離散分布には別の問題があります。二項確率質量は二項係数を必要としますが、「52個から26個を選ぶ」組み合わせの数は巨大な整数になります。これを直接計算しようとすると、適切な確率に引き戻すための除算を行う前に、分子が倍精度浮動小数点数（double）の上限をオーバーフローします。エンジンはこれを直接計算しません。ログガンマ関数を介して階乗を対数空間で計算し、その対数を加減算し、成功および失敗確率の対数を組み込み、最後に一度だけ指数関数を適用します
対数ガンマ関数自体はLanczos近似であり、正の実数軸全体にわたって正確であり、計算負荷も低いものです。最後のExpが実行されるまですべての巨大な値は対数として保持されるため、ルーチンが生成する最大値は、高々1である確率値そのものになります。ポアソン質量関数も同じ手法に従い、分母の階乗の代わりに対数ガンマ項を使用します。pがちょうど0または1である極端な境界ケースは個別処理されるため、コードがLn(0)を呼び出すことはありません。HotXLSはBINOM.DIST(5,10,0.5,FALSE)に対して0.2460938を返し、累積のPOISSON.DIST(2,2,TRUE)に対して0.6766764を返し、Excelが表示する桁数と一致します
順関数の曲線挟み込みによる逆関数の算出
逆分布関数は逆の問いを投げかけます。すなわち、確率値が与えられたとき、累積分布関数（CDF）の値がそれに等しくなるような変数を探すことです。このセットの中で、高速な直接計算式を持つ逆関数は1つだけです。標準正規分布の逆関数であるNORM.S.INVは、Acklamの有理式近似（中央領域と2つの尾部に分割され、全範囲にわたりほぼ倍精度の精度を持つ一対の多項式比）を使用します。これは反復処理を必要としない閉形式の評価です
他の逆関数にはそのような直接的な数式がないため、エンジンは数値的に逆関数を計算します。分布のサポート範囲から選択された下限値と上限値で答えを挟み込み、二分法を適用します。中点で順方向の累積分布関数（CDF）を評価し、目標の確率を内包し続けるようにいずれかの境界を移動させ、区間が十分に狭くなるまで繰り返します。ガンマおよびカイ二乗の逆関数の場合、探索範囲はゼロから開始し、形状パラメータとスケールパラメータから構築された余裕のある上限値を使用し、確率がまだ内包されていない場合は上限値を2倍に増やします。t分布の逆関数は外側に広がる対称な境界で挟み込み、F分布の逆関数は非負の区間で二分法を実行します。呼び出しごとに数十回のCDF評価が必要になりますが、スプレッドシートの処理速度から見れば無視できるコストであり、すべての逆関数が反転させる順関数と全く同じ精度を持つというメリットがあります。そのため、CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE)のようなラウンドトリップ計算が、極めて高い精度で0.7を返すことができます
尾部に隠れていた常用対数のバグ
これは長い間存続しやすいタイプのバグであるため、紹介する価値があります。Acklamの正規分布逆関数のルーチンには3つの分岐があります。確率が約0.025から0.975の間にある場合に使用される広い中央分岐は、ログ関数を一切含まない多項式比に値を流します。確率が極めて小さいか極めて大きい場合の2つの尾部分岐は、尾部が「pの自然対数のマイナス平方根」のように振る舞うため、まず入力値の対数を取得します
尾部分岐の初期バージョンでは、自然対数（ln）を使用すべき場所に常用対数（log10）を使用していました。この2つは約2.30の定数倍の差があるため、尾部の計算結果は常に大幅にズレていました。それでも、日常的なチェックはすべて中央付近の範囲で行われるため、簡易テストでは正常に見えました。NORM.S.INV(0.5)はゼロになり、NORM.S.INV(0.975)は教科書通りの1.959964になり、どちらも対数関数を一切呼び出さない中央の多項式によって処理されます。このエラーは、確率が尾部に入った場合（たとえば、本来は-3.0902323を返すべきNORM.S.INV(0.001)）に初めて現れ、自然対数と常用対数の比率だけズレた値が返されました。信頼区間ヘルパーなど、その尾部で正規分布逆関数に依存するすべての関数が、同じ偏りを引き継いでいました。この教訓は平凡ですが重いものです。分岐構造を持つ関数は、すべての分岐の内部にテストポイントを設ける必要があります。なぜなら、正しい共通パスは、バグのあるまれなパスを隠蔽してしまうからです。修正は常用対数から自然対数への1つのトークンの変更のみであり、尾部の値は即座にExcelの値と一致しました
xの符号がt分布の尾部を決定する
スチューデントのt累積分布関数には、間違えやすい微妙な性質があります。その値はdf / (df + x²)で評価された正則化不完全ベータ関数から得られますが、そのベータ値はxの絶対値を超える尾部の確率を示しており、xまでの累積確率ではありません。t分布の対称な形状は、変換処理がxがゼロのどちら側に位置するかに依存することを意味します
xがゼロを超える場合、累積確率は「1から対称尾部の半分を引いたもの」になります。xがゼロを下回る場合はその尾部の半分になります。ゼロのときは正確に2分の1です。ベータ値を直接返すと、分布の誤った側を報告することになり、ゼロ以外のxに対して曲線全体がズレてしまいます。右尾部および両側分布のバリエーションも同じ分岐の上に構築されているため、T.DIST.2T(1,1)は0.5を返し、T.DIST(1,1,TRUE)は0.75を返します。逆関数のT.INVはこの修正された累積分布関数（CDF）に対して二分法を実行するため、ラウンドトリップ計算が成り立ちます
セルからは何も見えませんが、それが期待される結果です。数式を入力すると、Excelと一致する値が出力されます。エンジンを独自のロジックで拡張する方法については、数式エンジンとカスタム関数のチュートリアルでカバーされており、数式が他のシートや名前定義された範囲を参照する仕組みについては、定義名とシート間数式に関する記事でカバーされています。これらはすべて、DelphiおよびC++Builder用のHotXLS spreadsheet componentに含まれており、他の記事で紹介されている読み書き、グラフ作成、およびフォーマット用APIと並んで提供されています