在儲存格中輸入 =NORM.DIST(115,100,15,TRUE),Excel 便會直接傳回 0.8413447。該呼叫讀起來像是一次查閱。但事實並非如此。在該單個數字的背後是累加常態分佈(一個沒有封閉形式的積分),而在 CHISQ.INV.RT 和 BETA.DIST 背後是仔細的函式庫必須評估的特殊函式,而不是手動近似。自稱與 Excel 相容的試算表元件必須重現這些值到 Excel 顯示的最後一位數字,這意味著要重現數值方法,而不僅僅是函式名稱。
HotXLS 實現了五十多個此類統計函式,而使它們保持正確的工作在公式列中幾乎完全看不見。這是對引擎如何計算它們的介紹:共享的特殊函式核心、保持算術穩定的分支決策,以及一個長期隱藏在尾部(因為常見情況從未觸及損壞程式碼行)的反常態分佈錯誤。
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 中運作,其餘的都會遵循相同的路徑。
其他一切都建立在之上的兩個函式
此集合中的大多數累加分佈不是透過對其自身定義求和或求積分來計算的。它們是從兩個特殊函式計算出來的:寫作 P(a, x) 的正規化下限不完全 gamma 函式,以及寫作 Ix(a, b) 的正規化不完全 beta 函式。在內部,這些是分派器所依賴的協助工具,且鏈條很短。卡方 CDF 是形狀為 df/2 且比例為 2 的 gamma CDF。gamma CDF 直接是 P(a, x)。t、F 和二項累加函式都是正規化不完全 beta 在正確引數處的值。Poisson 的 CDF 是上限不完全 gamma 函式 Q。將 gamma 和 beta 函式實現好,十幾個分佈將免費繼承其精確度。
「正規化」(regularized) 一詞是整個核心。原始的不完全 gamma 會像階乘一樣增長,且原始的 beta 積分在答案出現之前很久就可能發生下限溢位或溢位。正規化形式被完全 gamma 或 beta 除以,因此它們完全存在於零到一的區間內,這正是機率所佔據的範圍。這種標準化使得同一個常式可以為具有兩個自由度的卡方以及具有兩百個自由度的卡方提供服務,而不會使中間項超出雙精確度浮點數的末端。它也解釋了為什麼您不透過將密度項的長尾相加來計算 CDF:每一項都攜帶其自身的捨入誤差,誤差隨著級數的執行而累積,而正規化的特殊函式則藉由評估快速收斂的級數或連分數來完全避開求和。
對角線下方為級數,對角線上方為連分數
不完全 gamma 常式在計算任何內容之前會做出一個決策:它將 x 與 a + 1 進行比較。該邊界並非任意。當 x 相對於 a 較小時,P(a, x) 的冪級數展開收斂很快,而當 x 較大時,收斂很慢,最終毫無用處。連分數具有相反的特性。因此,引擎對於低於 a + 1 的 x 使用冪級數,對於等於或高於 a + 1 的 x 使用 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 正好為零或一,因此程式碼從不呼叫 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 毫釐不差。
隱藏在尾部的常用對數
這是一個值得講述的錯誤,因為它是那種能存活很長時間的錯誤。Acklam 反常態分佈常式有三個分支。寬闊的中央分支(每當機率介於約 0.025 和 0.975 之間時使用)在輸入端執行多項式比率,其中沒有任何對數。兩個尾部分支(用於非常小或非常大的機率)各首先對輸入取對數,因為尾部的行為類似於減去 p 的自然對數的平方根。
尾部分支的早期版本在屬於自然對數的地方使用了常用對數(底數為 10 的對數)。這兩者相差大約 2.30 的常數因子,因此尾部結果以一致且相當大的幅度出錯。然而,在每次隨意檢查中,該函式看起來都很正常,因為隨意檢查都發生在中間。NORM.S.INV(0.5) 是零,NORM.S.INV(0.975) 是教科書上的 1.959964,這兩者都執行完全不呼叫對數的中央多項式。只有在機率跨入尾部時(例如 NORM.S.INV(0.001),其必須傳回 -3.0902323,卻因自然對數與常用對數的比率偏差而傳回錯誤值),錯誤才會顯現。任何在其尾部相依於反常態分佈的函式(包括信賴區間協助工具)都繼承了相同的偏斜。教訓是平凡且昂貴的:具有分支結構的函式需要在每個分支內進行測試點測試,因為正確的常見路徑會欣然掩蓋損壞的罕見路徑。修正方法是將常用對數改為自然對數(只需變更一個標記),且尾部值立即與 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。