Funkcje statystyczne Excela w Delphi: NORM, CHISQ, BETA

Wpisz =NORM.DIST(115,100,15,TRUE) do komórki, a Excel bez zbędnych ceremonii zwróci wartość 0.8413447. Wywołanie to wygląda jak proste wyszukanie. Tak jednak nie jest. Za tą jedną liczbą kryje się dystrybuanta rozkładu normalnego — całka niemająca postaci zamkniętej, a za funkcjami CHISQ.INV.RT i BETA.DIST stoją funkcje specjalne, które dokładna biblioteka musi rzetelnie obliczyć, a nie tylko przybliżać na oko. Komponent arkusza kalkulacyjnego deklarujący zgodność z Excelem musi odtwarzać te wartości do ostatniej cyfry wyświetlanej przez Excel, co oznacza konieczność odtworzenia metod numerycznych, a nie tylko samych nazw funkcji.

HotXLS implementuje ponad pięćdziesiąt takich funkcji statystycznych, a praca zapewniająca ich poprawność jest prawie całkowicie niewidoczna z poziomu paska formuły. Oto przegląd sposobu, w jaki silnik je oblicza: omówienie wspólnego rdzenia funkcji specjalnych, decyzji o rozgałęzieniach zapewniających stabilność obliczeń oraz jednego błędu odwrotnego rozkładu normalnego, który przez długi czas ukrywał się w ogonie rozkładu, ponieważ typowe przypadki nigdy nie trafiały na wadliwą linię kodu.

Jedno wywołanie arkusza, pięćdziesiąt rozkładów pod spodem

Funkcje te obejmują rodziny rozkładów, po które najczęściej sięga się w arkuszach statystycznych. Jest tu rodzina rozkładu normalnego (NORM.DIST i NORM.S.DIST wraz z ich odwrotnościami); rodzina rozkładów gamma i chi-kwadrat (GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT); rodzina rozkładu beta (BETA.DIST i BETA.INV); rozkłady z próby (T.DIST, T.DIST.2T, F.DIST i F.INV); para rozkładów dyskretnych (BINOM.DIST i POISSON.DIST); a także pomocnicze funkcje wnioskowania statystycznego, takie jak CONFIDENCE.T i CONFIDENCE.NORM. Z perspektywy wywołującego każda z nich to pojedyncza formuła. Ustawiasz dane wejściowe w komórkach, zlecasz przeliczenie arkusza i odczytujesz wynik.

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;

Metoda Calculate skoroszytu kompiluje i oblicza formułę ad-hoc dla aktywnego arkusza, zwracając wartość typu Variant. Jeden szczegół potrafi zaskoczyć przy pierwszej próbie: parser formuł stojący za Calculate przyjmuje średnik jako separator argumentów, dlatego zapisuje się =SUM(A1;B1), a nie =SUM(A1,B1). Formuły zapisane w komórkach zachowują standardowy dla Excela przecinek. Ten sam ewaluator obsługuje wszystkie funkcje statystyczne opisane poniżej, więc jeśli jedna z nich działa w Calculate, reszta podąża tą samą ścieżką.

Dwie funkcje, na których opiera się cała reszta

Większość dystrybuant w tym zestawie nie jest obliczana poprzez sumowanie lub całkowanie ich własnych definicji. Są one obliczane przy użyciu dwóch funkcji specjalnych: znormalizowanej dolnej niepełnej funkcji gamma, zapisywanej jako P(a, x), oraz znormalizowanej niepełnej funkcji beta, oznaczanej jako I_x(a, b). Wewnętrznie są to pomocnicy, na których opierają się dyspozytorzy wywołań, a łańcuch zależności jest krótki. Dystrybuanta chi-kwadrat to dystrybuanta gamma z parametrem kształtu df/2 i skali 2. Dystrybuanta gamma to bezpośrednio P(a, x). Dystrybuanty t-Studenta, F-Snedecora oraz dwumianowe to wartości znormalizowanej niepełnej funkcji beta przy odpowiednich argumentach. Dystrybuanta Poissona to górna niepełna funkcja gamma Q. Dobra implementacja funkcji gamma i beta sprawia, że kilkanaście rozkładów dziedziczy ich dokładność za darmo.

Słowo "znormalizowana" (regularized) ma tu kluczowe znaczenie. Surowa niepełna funkcja gamma rośnie jak silnia, a surowa całka beta może ulec niedomiarowi (underflow) lub nadmiarowi (overflow) na długo przed uzyskaniem ostatecznego wyniku. Formy znormalizowane są dzielone przez pełną funkcję gamma lub beta, dzięki czemu ich wartości mieszczą się całkowicie w przedziale od zera do jednego, czyli dokładnie tam, gdzie leży wartość prawdopodobieństwa. Ta normalizacja sprawia, że ta sama procedura może obsłużyć rozkład chi-kwadrat z dwoma stopniami swobody oraz z dwustoma, bez ryzyka, że wyrażenia pośrednie wyjdą poza zakres typu double. Wyjaśnia to również, dlaczego nie oblicza się dystrybuanty (CDF) poprzez sumowanie długiego ogona czynników gęstości prawdopodobieństwa: każdy czynnik niesie własny błąd zaokrąglenia, błędy te kumulują się w miarę postępu serii, a znormalizowana funkcja specjalna całkowicie omija sumowanie, szacując szybko zbieżny szereg lub ułamek łańcuchowy.

Szereg poniżej przekątnej, ułamek łańcuchowy powyżej

Procedura niepełnej funkcji gamma podejmuje jedną decyzję przed wykonaniem jakichkolwiek obliczeń: porównuje x z wartością a + 1. Ta granica nie jest przypadkowa. Rozwinięcie w szereg potęgowy funkcji P(a, x) zbiega szybko, gdy x jest małe w stosunku do a, a wolno — ostatecznie bezużytecznie — gdy x jest duże. Ułamek łańcuchowy ma odwrotną charakterystykę. Silnik używa więc szeregu potęgowego dla x mniejszego od a + 1 oraz ułamka łańcuchowego Lentza dla x równego bądź większego od a + 1, dzięki czemu każda gałąź kodu wykonuje tylko te zadania, do których się najlepiej nadaje.

Ułamek łańcuchowy wymaga jednego zabezpieczenia. Metoda Lentza opiera się na utrzymywaniu bieżącego licznika i mianownika oraz odwracaniu mianownika na każdym kroku, a jeśli którykolwiek z nich zbliża się do zera, operacja odwracania ulega uszkodzeniu. Rozwiązaniem jest minimalny próg dolny: gdy czynnik pośredni spada poniżej około 1e-30, jest on dociskany (clamped) do wartości 1e-30, co utrzymuje skończoność rekurencji bez zakłócania zbieżnej wartości. To samo zabezpieczenie pojawia się w ułamku łańcuchowym niepełnej funkcji beta z tego samego powodu. Jest to mała stała wykonująca ogromną pracę — decyduje o różnicy między stabilnym szacowaniem a dzieleniem przez wartość nieodróżnialną od zera.

Górny ogon rozkładu, Q(a, x), to po prostu 1 minus P(a, x) i tak właśnie oblicza się dystrybuantę Poissona: prawdopodobieństwo wystąpienia co najwyżej k zdarzeń ze średnią λ wynosi Q(k + 1, λ). Skierowanie obliczeń przez górną niepełną funkcję gamma zamiast sumowania k + 1 czynników Poissona to kolejny wybór polegający na oszacowaniu jednego zbieżnego wyrażenia zamiast kumulowania wielu małych składników.

Masy dyskretne bez przepełnienia silni

Rozkłady dyskretne niosą inne zagrożenie. Prawdopodobieństwo dwumianowe wymaga współczynnika dwumianowego, a współczynnik dla wyboru dwudziestu sześciu elementów z pięćdziesięciu dwóch to ogromna liczba całkowita. Bezpośrednie jej wyznaczenie spowoduje przepełnienie licznika w typie double przed dokonaniem dzielenia, które sprowadziłoby wartość do rozsądnego prawdopodobieństwa. Silnik nigdy nie oblicza tego bezpośrednio. Wyznacza silnie w przestrzeni logarytmicznej przy użyciu funkcji logarytmu gamma, dodaje i odejmuje logarytmy, uwzględnia logarytmy prawdopodobieństwa sukcesu i porażki, a na samym końcu wykonuje operację potęgowania (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));

Sama funkcja logarytmu gamma to przybliżenie Lanczosa, dokładne na całej dodatniej półosi i tanie w obliczeniach. Ponieważ każda duża wartość jest przechowywana jako logarytm aż do ostatecznego wywołania Exp, największą liczbą rzeczywistą, jaką tworzy ta procedura, jest samo prawdopodobieństwo, czyli wartość wynosząca co najwyżej jeden. Funkcja masy prawdopodobieństwa Poissona opiera się na tym samym schemacie, gdzie pojedyncze wyrażenie logarytmu gamma zastępuje silnię w mianowniku. Postacie zamknięte są traktowane jako przypadki szczególne na granicach przedziału, gdy p wynosi dokładnie zero lub jeden, dzięki czemu kod nigdy nie wywołuje Ln(0). HotXLS zwraca wartość 0.2460938 dla BINOM.DIST(5,10,0.5,FALSE) oraz 0.6766764 dla dystrybuanty POISSON.DIST(2,2,TRUE), co odpowiada wynikom z Excela na wszystkich wyświetlanych pozycjach dziesiętnych.

Odwrotności poprzez osaczanie krzywej w przód

Odwrotna funkcja rozkładu zadaje przeciwne pytanie: dla danego prawdopodobieństwa znajdź wartość, dla której dystrybuanta (CDF) jest mu równa. Tylko jedna odwrotność w tym zestawie posiada szybki bezpośredni wzór. Funkcja NORM.S.INV, czyli odwrotny standardowy rozkład normalny, wykorzystuje przybliżenie wymierne Acklama — parę ilorazów wielomianów dokładną z grubsza do precyzji typu double w całym zakresie, podzieloną na obszar centralny i dwa ogony. Jest to wyznaczenie w postaci zamkniętej bez iteracji.

Inne odwrotności nie mają takiego wzoru, dlatego silnik wyznacza je numerycznie. Ogranicza odpowiedź dolną i górną granicą wybraną z nośnika rozkładu, a następnie stosuje metodę bisekcji: oblicza dystrybuantę w punkcie środkowym, przesuwa granicę, która utrzymuje docelowe prawdopodobieństwo wewnątrz przedziału, i powtarza to do uzyskania odpowiednio wąskiego przedziału. Dla odwrotności gamma i chi-kwadrat przedział zaczyna się od zera i bezpiecznego górnego szacunku wyznaczonego z parametrów kształtu i skali, podwajając górną granicę, jeśli prawdopodobieństwo nie zostało jeszcze ujęte. Odwrotność rozkładu t-Studenta ogranicza symetryczne przedziały rozszerzające się na zewnątrz; odwrotność rozkładu F-Snedecora stosuje bisekcję na przedziale nieujemnym. Koszt to kilkadziesiąt obliczeń dystrybuanty na wywołanie, co jest niezauważalne przy szybkości działania arkusza, a korzyścią jest to, że każda odwrotność jest dokładnie tak dokładna, jak funkcja pierwotna, którą odwraca. Dlatego operacja dwukierunkowa, taka jak CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE), zwraca wartość 0.7 z minimalnym błędem.

Logarytm o podstawie 10, który ukrył się w ogonie

Wczesna wersja gałęzi ogona obliczała logarytm o podstawie 10 tam, gdzie powinien znajdować się logarytm naturalny. Różnią się one o stały współczynnik wynoszący około 2.30, przez co wyniki w ogonach były obarczone stałym, znaczącym błędem. Mimo to funkcja wyglądała poprawnie przy każdym pobieżnym teście, ponieważ takie testy dotyczą wartości środkowych. Wywołanie NORM.S.INV(0.5) daje zero, NORM.S.INV(0.975) zwraca podręcznikowe 1.959964, a oba te wywołania przechodzą przez centralny wielomian, który w ogóle nie korzysta z logarytmu. Błąd ujawniał się dopiero po przejściu prawdopodobieństwa w obszar ogona, na przykład przy NORM.S.INV(0.001), które powinno zwrócić -3.0902323, a dawało wynik przesunięty o stosunek logarytmu naturalnego do logarytmu o podstawie 10. Każda funkcja zależna od odwrotności rozkładu normalnego w jego ogonach, w tym funkcje przedziału ufności, dziedziczyła to samo zniekształcenie. Wniosek jest prosty i kosztowny: funkcja o strukturze rozgałęzionej wymaga punktów testowych w każdej gałęzi, ponieważ poprawna powszechna ścieżka z łatwością zamaskuje uszkodzoną ścieżkę rzadką. Rozwiązaniem była zamiana jednego tokenu z logarytmu dziesiętnego na naturalny, co spowodowało powrót wartości w ogonach do zgodności z Excelowymi.

Znak x decyduje o ogonie rozkładu t-Studenta

Dystrybuanta rozkładu t-Studenta niesie subtelność, którą łatwo odwrócić. Jej wartość pochodzi ze znormalizowanej niepełnej funkcji beta wyznaczonej dla df / (df + x²), jednak ta wartość beta reprezentuje prawdopodobieństwo w ogonie poza zakresem x, a nie dystrybuantę do punktu x. Symetryczny kształt rozkładu t-Studenta sprawia, że konwersja zależy od tego, po której stronie zera leży wartość 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

Dla x powyżej zera dystrybuanta wynosi jeden minus połowa symetrycznego ogona; dla x poniżej zera to połowa tego ogona; w zerze wynosi dokładnie jedną drugą. Zwrócenie wartości beta bezpośrednio dałoby błędną stronę rozkładu, przesuniętą o cały korpus krzywej dla każdego niezerowego x. Warianty prawostronny i dwustronny bazują na tej samej gałęzi kodu, dlatego T.DIST.2T(1,1) zwraca 0.5, a T.DIST(1,1,TRUE) daje 0.75, a odwrotność T.INV stosuje bisekcję wobec tej poprawionej dystrybuanty, co zamyka dwukierunkową operację.

Nic z tego nie jest widoczne z poziomu komórki i taki właśnie jest cel. Zapisujesz formułę i odczytujesz liczbę zgodną z Excelem. Jeśli rozbudowujesz silnik o własną logikę, mechanizmy rejestrowania funkcji opisano w naszym przewodniku po silniku formuł i własnych funkcjach, a sposób, w jaki formuły sięgają do innych arkuszy i nazwanych zakresów, omówiono w artykule na temat nazw zdefiniowanych i formuł międzyarkuszowych. Wszystko to jest dostarczane w pakiecie komponentu arkusza kalkulacyjnego HotXLS dla Delphi i C++Builder obok interfejsów API do odczytu, zapisu, tworzenia wykresów i formatowania omówionych w innych miejscach na tym blogu.