Nhập =NORM.DIST(115,100,15,TRUE) vào một ô và Excel trả về 0.8413447 mà không cần nghi lễ nào. Cuộc gọi này đọc lên giống như một tra cứu bảng. Nhưng không phải vậy. Đằng sau con số đơn độc đó là phân phối chuẩn tích lũy, một tích phân không có dạng đóng, và đằng sau CHISQ.INV.RT cùng BETA.DIST là các hàm đặc biệt mà một thư viện cẩn thận phải đánh giá, chứ không phải xấp xỉ bằng tay. Một thành phần bảng tính tuyên bố khả năng tương thích với Excel phải tái tạo các giá trị này đến chữ số cuối cùng mà Excel hiển thị, điều đó có nghĩa là tái tạo các phương pháp số, chứ không chỉ là tên hàm.
HotXLS triển khai hơn năm mươi hàm thống kê này, và công việc giúp chúng chính xác gần như hoàn toàn vô hình từ thanh công thức. Đây là một chuyến tham quan về cách công cụ tính toán chúng: lõi hàm đặc biệt dùng chung, các quyết định rẽ nhánh giúp giữ cho số học ổn định, và một lỗi phân phối chuẩn ngược ẩn náu ở phần đuôi (tail) trong một thời gian dài vì trường hợp phổ biến không bao giờ chạm đến dòng mã bị hỏng đó.
Một cuộc gọi trang tính, năm mươi phân phối đằng sau nó
Các hàm này bao gồm các họ hàm mà một bảng tính thống kê hướng tới. Có họ phân phối chuẩn, NORM.DIST và NORM.S.DIST cùng các hàm ngược của chúng; họ gamma và chi-bình phương, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; họ beta, BETA.DIST và BETA.INV; các phân phối lấy mẫu T.DIST, T.DIST.2T, F.DIST và F.INV; cặp phân phối rời rạc BINOM.DIST và POISSON.DIST; và các trình hỗ trợ suy luận như CONFIDENCE.T và CONFIDENCE.NORM. Từ góc độ người gọi, mỗi cái là một công thức duy nhất. Bạn thiết lập đầu vào trong các ô, yêu cầu bảng tính đánh giá, và đọc kết quả.
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;Phương thức Calculate trên bảng tính sẽ biên dịch và đánh giá một công thức đặc biệt (ad-hoc) đối chiếu với trang tính trực tiếp và trả về một Variant. Một chi tiết gây bối rối cho mọi người trong lần thử đầu tiên: trình phân tích công thức đằng sau Calculate lấy dấu chấm phẩy làm dấu phân cách đối số của nó, vì vậy nó là =SUM(A1;B1), chứ không phải =SUM(A1,B1). Các công thức ô được lưu trữ giữ dấu phẩy chuẩn Excel. Cùng một bộ đánh giá sẽ điều phối mọi hàm thống kê bên dưới, vì vậy một khi một trong các hàm này hoạt động trong Calculate, các hàm còn lại sẽ đi theo cùng một đường dẫn.
Hai hàm mà mọi thứ khác được xây dựng trên đó
Hầu hết các phân phối tích lũy trong tập hợp này không được tính toán bằng cách tổng hợp hoặc tích phân các định nghĩa của chính chúng. Chúng được tính toán từ hai hàm đặc biệt: hàm gamma không đầy đủ dưới chuẩn hóa (regularized lower incomplete gamma), viết là P(a, x), và hàm beta không đầy đủ chuẩn hóa (regularized incomplete beta), viết là Ix(a, b). Về mặt nội bộ, đây là các trình hỗ trợ mà các trình điều phối dựa vào, và chuỗi liên kết rất ngắn. Hàm chi-bình phương CDF là hàm gamma CDF với hình dạng df/2 và tỷ lệ 2. Hàm gamma CDF là P(a, x) trực tiếp. Các hàm tích lũy t, F, và nhị phân đều là các giá trị của hàm beta không đầy đủ chuẩn hóa tại các đối số phù hợp. Hàm Poisson CDF là hàm gamma không đầy đủ trên Q. Triển khai tốt các hàm gamma và beta thì một tá phân phối sẽ kế thừa độ chính xác của chúng miễn phí.
Từ "chuẩn hóa" (regularized) là toàn bộ mấu chốt vấn đề. Hàm gamma không đầy đủ thô tăng lên giống như một giai thừa và tích phân beta thô có thể bị tràn ngược (underflow) hoặc tràn thuận (overflow) từ lâu trước khi có câu trả lời. Các dạng chuẩn hóa được chia cho hàm gamma hoặc beta đầy đủ, vì vậy chúng nằm hoàn toàn trong khoảng từ 0 đến 1, chính xác là phạm vi mà một xác suất chiếm giữ. Sự chuẩn hóa đó là điều cho phép cùng một quy trình phục vụ phân phối chi-bình phương với hai bậc tự do và một phân phối với hai trăm bậc tự do mà không làm các số hạng trung gian chạy ra ngoài giới hạn của một số double. Nó cũng giải thích tại sao bạn không tính toán một CDF bằng cách cộng dồn một phần đuôi dài các số hạng mật độ: mỗi số hạng mang lỗi làm tròn riêng của nó, các lỗi này tích lũy khi chuỗi chạy, và hàm đặc biệt chuẩn hóa tránh hoàn toàn phép tổng bằng cách đánh giá một chuỗi hội tụ nhanh hoặc một phân số liên tục (continued fraction) thay thế.
Chuỗi bên dưới đường chéo, phân số liên tục bên trên nó
Quy trình gamma không đầy đủ đưa ra một quyết định trước khi nó tính toán bất cứ thứ gì: nó so sánh x với a + 1. Ranh giới đó không phải là tùy tiện. Việc mở rộng chuỗi lũy thừa (power-series expansion) của P(a, x) hội tụ nhanh chóng khi x nhỏ so với a, và chậm chạp, cuối cùng là vô dụng, khi x lớn. Phân số liên tục có đặc điểm ngược lại. Vì vậy, công cụ sử dụng chuỗi lũy thừa cho x dưới a + 1 và phân số liên tục Lentz cho x ở hoặc trên a + 1, và mỗi nhánh chỉ được yêu cầu thực hiện phần việc mà nó làm tốt.
Phân số liên tục cần một biện pháp bảo vệ. Phương pháp Lentz hoạt động bằng cách mang theo một tử số và mẫu số lũy tiến và nghịch đảo mẫu số trên mỗi bước, và nếu một trong hai tiến gần đến không, phép nghịch đảo sẽ bùng nổ. Giải pháp là một mức sàn nhỏ: bất cứ khi nào một số hạng trung gian giảm xuống dưới khoảng 1e-30 về độ lớn, nó sẽ bị kẹp lại (clamp) ở mức 1e-30, giúp giữ cho quan hệ đệ quy có hạn mà không làm xáo trộn giá trị hội tụ. Kẹp tương tự cũng xuất hiện trong phân số liên tục của hàm beta không đầy đủ vì cùng một lý do. Đó là một hằng số nhỏ thực hiện công việc chịu tải quan trọng, sự khác biệt giữa một đánh giá ổn định và một phép chia cho thứ không thể phân biệt được với số không.
Phần đuôi trên, Q(a, x), đơn giản là 1 trừ đi P(a, x), và đó là cách nhánh tích lũy của Poisson được tính toán: xác suất của tối đa k sự kiện với giá trị trung bình λ là Q(k + 1, λ). Định tuyến nó qua hàm gamma không đầy đủ trên thay vì tính tổng k + 1 số hạng Poisson, một lần nữa, là lựa chọn để đánh giá một biểu thức hội tụ thay vì tích lũy nhiều biểu thức nhỏ.
Khối lượng rời rạc không bị tràn giai thừa
Các phân phối rời rạc gây ra một mối nguy hại khác. Một khối xác suất nhị phân (binomial probability mass) liên quan đến một hệ số nhị phân, và hệ số cho năm-mươi-hai-chọn-hai-mươi-sáu là một số nguyên khổng lồ. Nếu tạo ra nó trực tiếp, tử số sẽ làm tràn kiểu double trước phép chia vốn đưa nó trở lại một xác suất hợp lý. Công cụ không bao giờ tạo ra nó. Nó tính toán các giai thừa trong không gian log thông qua hàm log-gamma, cộng và trừ các log, tích hợp log của xác suất thành công và thất bại, và thực hiện lũy thừa một lần duy nhất ở phần cuối cùng.
// 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));Bản thân hàm log-gamma là một phép xấp xỉ Lanczos, chính xác trên toàn bộ trục dương và rẻ để đánh giá. Bởi vì mọi lượng lớn đều được giữ dưới dạng logarit của nó cho đến hàm Exp cuối cùng, số lớn nhất mà quy trình từng hiện thực hóa chính là xác suất đó, tối đa là một. Hàm khối lượng Poisson tuân theo cùng một công thức, với thuật ngữ log-gamma đơn độc đại diện cho giai thừa ở mẫu số. Các dạng đóng được xử lý trường hợp đặc biệt tại các biên, nơi p chính xác bằng không hoặc một, vì vậy mã không bao giờ gọi Ln(0). HotXLS trả về 0.2460938 cho BINOM.DIST(5,10,0.5,FALSE) và 0.6766764 cho hàm tích lũy POISSON.DIST(2,2,TRUE), khớp với Excel qua các chữ số mà nó in ra.
Các hàm ngược bằng cách kẹp đường cong thuận
Một hàm phân phối ngược hỏi câu hỏi ngược lại: cho một xác suất, hãy tìm giá trị mà CDF của nó bằng xác suất đó. Chỉ có một hàm ngược trong tập hợp này có công thức trực tiếp nhanh chóng. NORM.S.INV, hàm ngược phân phối chuẩn tiêu chuẩn, sử dụng phép xấp xỉ hữu tỷ Acklam (Acklam rational approximation), một cặp tỷ lệ đa thức chính xác đến khoảng độ chính xác của một số double trên toàn bộ phạm vi, được chia thành một vùng trung tâm và hai vùng đuôi. Nó là một đánh giá dạng đóng không có vòng lặp đệ quy.
Các hàm ngược khác không có công thức như vậy, vì vậy công cụ đảo ngược chúng bằng phương pháp số. Nó kẹp (bracket) câu trả lời bằng một biên thấp và biên cao được chọn từ vùng hỗ trợ của phân phối, sau đó chia đôi (bisect): đánh giá CDF thuận tại điểm giữa, di chuyển biên nào giữ cho xác suất mục tiêu được bao quanh, và lặp lại cho đến khi khoảng hẹp lại. Đối với các hàm ngược gamma và chi-bình phương, kẹp bắt đầu tại số không và một ước lượng trên hào phóng được xây dựng từ hình dạng và tỷ lệ, nhân đôi biên trên nếu xác suất chưa được bao quanh. Hàm ngược t kẹp các biên đối xứng mở rộng ra ngoài; hàm ngược F chia đôi trên một khoảng không âm. Chi phí là một vài tá đánh giá CDF cho mỗi cuộc gọi, điều này vô hình ở tốc độ bảng tính, và lợi ích là mỗi hàm ngược chính xác giống như hàm thuận mà nó đảo ngược. Đó là lý do tại sao một chuyến đi khứ hồi như CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) trả về 0.7 chính xác đến từng sợi tóc.
Logarit cơ số 10 ẩn náu ở phần đuôi
Đây là lỗi đáng được kể lại, bởi vì nó là loại lỗi tồn tại rất lâu. Quy trình phân phối chuẩn ngược Acklam có ba nhánh. Nhánh trung tâm rộng, được sử dụng bất cứ khi nào xác suất nằm trong khoảng từ 0.025 đến 0.975, chạy đầu vào qua một tỷ lệ đa thức không có logarit ở bất kỳ đâu trong đó. Hai nhánh đuôi, cho các xác suất rất nhỏ hoặc rất lớn, mỗi nhánh lấy một logarit của đầu vào trước, bởi vì phần đuôi hoạt động giống như căn bậc hai của trừ logarit tự nhiên của p.
Một phiên bản ban đầu của nhánh đuôi đã lấy logarit cơ số 10 ở nơi mà logarit tự nhiên thuộc về. Cả hai khác nhau bởi một hệ số không đổi khoảng 2.30, do đó các kết quả phần đuôi đã bị sai lệch bởi một biên độ nhất quán và khá lớn. Ấy thế mà hàm trông vẫn ổn trong mọi kiểm tra thông thường, bởi vì các kiểm tra thông thường nằm ở giữa. NORM.S.INV(0.5) bằng không, NORM.S.INV(0.975) là giá trị sách giáo khoa 1.959964, and cả hai đều chạy qua đa thức trung tâm không bao giờ gọi logarit. Lỗi chỉ xuất hiện một khi xác suất đi vào vùng đuôi, ví dụ NORM.S.INV(0.001), vốn phải trả về -3.0902323 nhưng thay vào đó lại quay về mức gần sai lệch tỷ lệ logarit tự nhiên so với cơ số 10. Bất kỳ hàm nào phụ thuộc vào phân phối chuẩn ngược ở vùng đuôi của nó, bao gồm cả các trình hỗ trợ khoảng tin cậy (confidence-interval), đều kế thừa cùng một độ lệch. Bài học này rất tầm thường nhưng đắt giá: một hàm có cấu trúc rẽ nhánh cần các điểm kiểm tra bên trong mỗi nhánh, bởi vì một đường dẫn chung chính xác sẽ che giấu một đường dẫn hiếm gặp bị hỏng. Bản sửa lỗi là thay đổi một token duy nhất từ log cơ số 10 sang log tự nhiên, và các giá trị phần đuôi đã khớp ngay với Excel.
Dấu của x quyết định phần đuôi của phân phối t
Hàm tích lũy t-Student mang một sự tinh tế rất dễ làm ngược. Giá trị của nó đến từ hàm beta không đầy đủ chuẩn hóa được đánh giá tại df / (df + x²), nhưng giá trị beta đó là xác suất ở phần đuôi vượt ra ngoài độ lớn của x, chứ không phải xác suất tích lũy tính đến x. Hình dạng đối xứng của phân phối t có nghĩa là chuyển đổi phụ thuộc vào việc x nằm ở bên nào của số không.
// 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Đối với x lớn hơn không, xác suất tích lũy là một trừ đi một nửa phần đuôi đối xứng; đối với x dưới không, nó bằng một nửa phần đuôi đó; tại số không, nó chính xác bằng một nửa. Trả về giá trị beta trực tiếp và bạn sẽ báo cáo sai phần của phân phối, lệch đi toàn bộ thân đường cong đối với bất kỳ x khác không nào. Các biến thể đuôi phải và hai đuôi xây dựng trên cùng một nhánh, đó là lý do tại sao T.DIST.2T(1,1) trả về 0.5 và T.DIST(1,1,TRUE) là 0.75, và hàm ngược T.INV chia đôi chống lại CDF đã sửa này để chuyến khứ hồi khép lại.
Không điều nào trong số này có thể nhìn thấy từ ô dữ liệu, và đó là kết quả mong muốn. Bạn viết một công thức và đọc một số khớp với Excel. Nếu bạn đang mở rộng công cụ với logic của riêng mình, cơ chế đăng ký một hàm được đề cập trong our walkthrough of the formula engine and custom functions, và cách các công thức tiếp cận các trang tính và phạm vi được đặt tên được đề cập trong the article on defined names and cross-sheet formulas. All of it ships inside the HotXLS spreadsheet component for Delphi and C++Builder, alongside the reading, writing, charting, and formatting APIs covered elsewhere on this blog.