在儲存格中輸入 =NORM.DIST(115,100,15,TRUE),Excel 便會不動聲色地回傳 0.8413447。這個呼叫看起來像查表,但其實不是。那個數字背後是常態分布的累積分布,一個沒有封閉解的積分;而 CHISQ.INV.RT 和 BETA.DIST 背後則是特殊函式,細心的函式庫必須真正去計算,而不是手工近似。宣稱與 Excel 相容的試算表元件,必須把這些值算到 Excel 顯示的最後一位都一致,這表示要重現的是數值方法,而不只是函式名稱
HotXLS 實作了五十多個這類統計函式,而讓它們保持正確的工作,在公式列裡幾乎完全看不出來。這是一趟檢視引擎如何計算它們的導覽:共用的特殊函式核心、讓算術保持穩定的分支判斷,以及一個長期躲在尾端的反常態分布 bug,因為常見情況從來沒碰到那條壞掉的程式碼
一個工作表呼叫,背後五十種分布
這些函式涵蓋統計工作簿會用到的家族。有常態家族,包含 NORM.DIST、NORM.S.DIST 及其反函式;有 gamma 與卡方家族,包含 GAMMA.DIST、CHISQ.DIST、CHISQ.DIST.RT、CHISQ.INV.RT;有 beta 家族,包含 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 標準的逗號。下面所有統計函式都由同一個評估器分派,所以只要其中一個在 Calculate 裡能運作,其餘都走同一條路
支撐一切的兩個函式
這組函式裡的大多數累積分布,不是透過把自己的定義加總或積分來算的。它們是從兩個特殊函式算出來的:正規化下不完全 gamma,寫作 P(a, x),以及正規化不完全 beta,寫作 Ix(a, b)。在內部,這些就是分派器依賴的輔助工具,而鏈條很短。卡方 CDF 就是形狀為 df/2、尺度為 2 的 gamma CDF。gamma CDF 直接就是 P(a, x)。t、F 和二項累積函式都是在正確參數上求得的正規化不完全 beta 值。Poisson 的 CDF 則是上不完全 gamma Q。把 gamma 和 beta 函式做好,十幾種分布就能自動繼承它們的精確度
「正規化」這個字就是重點。原始的不完全 gamma 會像階乘一樣成長,而原始的 beta 積分在答案出現之前很久就可能發生下溢或上溢。正規化形式會除以完整 gamma 或 beta,所以它們完全落在 0 到 1 的區間內,這正是機率所在的範圍。這種標準化讓同一個常式可以同時服務自由度為 2 的卡方和自由度為 200 的卡方,而中間項不會先把雙精確度浮點數逼到極限。它也說明了為什麼您不會把一長串密度項相加來算 CDF:每一項都帶著自己的捨入誤差,級數跑得越長,誤差就越累積,而正規化特殊函式則直接用快速收斂的級數或連分數來繞過整個求和
對角線下用級數,對角線上用連分數
不完全 gamma 常式在計算任何東西之前先做一個決定:把 x 和 a + 1 比較。這條界線不是任意的。當 x 相對於 a 很小時,P(a, x) 的冪級數展開收斂很快;當 x 很大時,收斂很慢,最後甚至沒有用。連分數則剛好相反。所以引擎在 x 低於 a + 1 時用冪級數,當 x 等於或高於 a + 1 時用 Lentz 連分數,讓每個分支只做它擅長的工作
連分數需要一個防護。Lentz 方法會一路帶著分子和分母前進,並在每一步把分母取倒數;如果其中任何一個接近零,取倒數就會炸掉。修正方式是一個很小的下限:只要中間項的大小降到大約 1e-30 以下,就把它夾到 1e-30,這樣可以讓遞迴保持有限,而不會干擾收斂值。不完全 beta 的連分數也出現同樣的夾限,原因一樣。這是一個小常數在做承重工作,差別就在穩定評估和除以幾乎等於零的東西之間
上尾 Q(a, x) 就是 1 減去 P(a, x),Poisson 的累積分支也是這樣算的:平均值為 λ、最多 k 個事件的機率就是 Q(k + 1, λ)。把它送進上不完全 gamma,而不是把 k + 1 個 Poisson 項一個一個加起來,仍然是在選擇只評估一個收斂式,而不是累積很多小數
沒有階乘溢位的離散質量
離散分布帶來另一種危險。二項機率質量會牽涉二項式係數,而五十二選二十六的係數是一個巨大的整數。若直接形成它,在分母把它拉回合理機率之前,分子就會先溢出雙精確度浮點數。引擎從來不直接形成它。它透過 log-gamma 函式把階乘算在對數空間中,對對數做加減,併入成功與失敗機率的對數,然後直到最後才做一次指數運算
// 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。Poisson 的質量函式遵循同樣的配方,只是用單一 log-gamma 項取代分母裡的階乘。邊界上的閉式會被特別處理,也就是 p 正好是 0 或 1 的情況,這樣程式碼就不會去呼叫 Ln(0)。HotXLS 對 BINOM.DIST(5,10,0.5,FALSE) 回傳 0.2460938,對累積的 POISSON.DIST(2,2,TRUE) 回傳 0.6766764,和 Excel 列出的數字一致
用夾擠正向曲線來求反函式
反分布函式問的是相反的問題:給定一個機率,找出那個 CDF 等於該機率的值。這組裡只有一個反函式有快速直接公式。NORM.S.INV,也就是反標準常態分布,使用 Acklam 有理近似,一對多項式比值,在整個範圍內大致精確到雙精確度浮點數的精度,並分成中央區域和兩個尾端。這是一個沒有迭代的閉式評估
其他反函式沒有這種公式,所以引擎用數值方式反解。它先用分布支援範圍內選出的低界與高界把答案夾住,然後做二分:在中點評估正向 CDF,把能維持目標機率仍被包住的那一側界線移過去,重複直到區間夠窄。對 gamma 和卡方反函式,夾擠從零開始,再配上一個根據形狀和尺度算出的寬鬆上界,如果機率還沒被夾住,就把上界加倍。t 反函式會用向外擴張的對稱界線來夾擠;F 反函式則在非負區間上二分。代價是每次呼叫要做幾十次 CDF 評估,但在試算表速度下幾乎看不出來;好處則是每個反函式都和它反解的正向函式一樣精確。這就是為什麼像 CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) 這種往返,會幾乎精準回到 0.7
藏在尾端的底數十對數
這裡有個值得一提的 bug,因為它能活很久。Acklam 反常態分布常式有三個分支。寬廣的中央分支,凡機率介於大約 0.025 與 0.975 之間時都會用,它把輸入送進一個多項式比值,裡面完全沒有對數。兩個尾端分支則分別先對非常小或非常大的機率取對數,因為尾端的行為類似於負的自然對數再開根號
早期的尾端分支版本在應該用自然對數的地方用了底數 10 的對數。兩者差了一個大約 2.30 的常數因子,所以尾端結果會以一致而明顯的幅度出錯。可是在任何隨手檢查中,這個函式看起來都沒問題,因為隨手檢查都發生在中間。NORM.S.INV(0.5) 是 0,NORM.S.INV(0.975) 是教科書上的 1.959964,這兩個都走中央多項式,根本不會呼叫對數。錯誤只有在機率進入尾端時才出現,例如 NORM.S.INV(0.001),理論上應該回傳 -3.0902323,結果卻因自然對數和底數 10 對數的比率而偏掉。任何在尾端依賴反常態分布的函式,包括信賴區間輔助函式,都繼承了同樣的偏斜。教訓平凡卻昂貴:有分支結構的函式,必須在每個分支都放測試點,因為一條正確的常見路徑會很樂意掩蓋一條壞掉的罕見路徑。修正只是一個標記的更動,從底數 10 對數改成自然對數,尾端值就和 Excel 對齊了
x 的正負決定 t 分布的尾端
Student t 累積函式有個很容易弄反的細節。它的值來自在 df / (df + x²) 處評估的正規化不完全 beta,但那個 beta 值其實是 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,則是尾端的一半;在零點時正好是二分之一。直接回傳 beta 值就會報到分布的錯邊,對任何非零 x 都會偏掉整條曲線主體。右尾與雙尾變體都建立在同一個分支上,所以 T.DIST.2T(1,1) 會回到 0.5,而 T.DIST(1,1,TRUE) 會回到 0.75;反函式 T.INV 也會以這個修正後的 CDF 進行二分,因此往返就能閉合
從儲存格裡看不出這些,而這正是目的。您寫下一個公式,讀到一個和 Excel 一致的數字。如果您要用自己的邏輯擴充引擎,註冊函式的機制在我們對公式引擎與自訂函式的逐步解說中有說明,而公式如何跨工作表與具名範圍存取,則在關於定義名稱與跨工作表公式的文章中說明。所有這些都包含在適用於 Delphi 和 C++Builder 的 HotXLS 試算表元件 中,並與本部落格其他地方介紹的讀取、寫入、圖表和格式化 API 一起提供