Excel statističke funkcije u Delphi-ju: NORM, CHISQ, BETA

Upišite =NORM.DIST(115,100,15,TRUE) u ćeliju i Excel bez ikakve ceremonije vraća 0.8413447. Taj poziv izgleda kao obično pretraživanje. Ali nije. Iza tog jednog broja krije se kumulativna normalna raspodela, integral koji nema zatvoreni oblik (closed form), a iza CHISQ.INV.RT i BETA.DIST stoje specijalne funkcije koje pažljivo napisana biblioteka mora da evaluira, a ne da aproksimira ručno. Komponenta za tabele koja tvrdi da je kompatibilna sa Excel-om mora da reprodukuje ove vrednosti do poslednje cifre koju Excel prikazuje, što znači da mora da reprodukuje numeričke metode, a ne samo nazive funkcija

HotXLS implementira više od pedeset ovih statističkih funkcija, a rad koji ih čini ispravnim je gotovo potpuno nevidljiv iz trake za formule. Ovo je pregled načina na koji ih mehanizam računa: zajedničko jezgro specijalnih funkcija, odluke o granama koje održavaju aritmetiku stabilnom i jedna greška u inverznoj normalnoj raspodeli koja se dugo krila u repu jer uobičajeni slučaj nikada nije dotakao pokvarenu liniju koda

Jedan poziv radnog lista, pedeset raspodela iza njega

Ove funkcije obuhvataju porodice koje se koriste u statističkim radnim sveskama. Tu je porodica normalne raspodele, NORM.DIST i NORM.S.DIST sa svojim inverznim funkcijama; porodica gama i hi-kvadrat raspodele, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; porodica beta raspodele, BETA.DIST i BETA.INV; raspodele uzoraka T.DIST, T.DIST.2T, F.DIST i F.INV; diskretni par BINOM.DIST i POISSON.DIST; kao i pomoćne funkcije za zaključivanje kao što su CONFIDENCE.T i CONFIDENCE.NORM. Iz perspektive pozivaoca, svaka od njih je jedna formula. Podesite ulaze u ćelijama, zatražite od radne sveske da izvrši evaluaciju i pročitate rezultat

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 na radnoj svesci kompajlira i evaluira ad-hoc formulu u odnosu na aktivni list i vraća tip Variant. Jedan detalj zbunjuje ljude pri prvom pokušaju: parser formula iza metode Calculate uzima tačku i zarez (semicolon) kao separator argumenata, tako da se piše =SUM(A1;B1), a ne =SUM(A1,B1). Sačuvane formule ćelija zadržavaju standardni zarez u Excel-u. Isti evaluator otprema svaku statističku funkciju u nastavku, pa kada jedna od njih proradi u metodi Calculate, sve ostale prate istu putanju

Dve funkcije na kojima se sve ostalo gradi

Većina kumulativnih raspodela u ovom skupu se ne računa sabiranjem ili integrisanjem sopstvenih definicija. One se računaju iz dve specijalne funkcije: regularizovane donje nepotpune gama funkcije, koja se piše kao P(a, x), i regularizovane nepotpune beta funkcije, koja se piše kao I_x(a, b). Interno, to su pomoćne funkcije na koje se dispečeri oslanjaju, i taj lanac je kratak. CDF (funkcija raspodele) hi-kvadrata je CDF gama funkcije sa oblikom df/2 i skalom 2. CDF gama funkcije je direktno P(a, x). Kumulativne funkcije za t, F i binomnu raspodelu su sve vrednosti regularizovane nepotpune beta funkcije pri odgovarajućim argumentima. Poisson-ov CDF je gornja nepotpuna gama funkcija Q. Implementirajte dobro gama i beta funkcije i desetak raspodela će besplatno naslediti njihovu tačnost

Reč "regularizovana" (regularized) je suština svega. Sirova nepotpuna gama raste poput faktorijela, a sirovi beta integral može pasti ispod donje granice (underflow) ili premašiti gornju granicu (overflow) davno pre nego što sam rezultat to uradi. Regularizovani oblici su podeljeni sa potpunom gama ili beta funkcijom, tako da u potpunosti žive u intervalu od nule do jedinice, što je upravo opseg koji verovatnoća zauzima. Ta normalizacija je ono što omogućava istoj rutini da usluži hi-kvadrat sa dva stepena slobode i onaj sa dve stotine, a da intermediate (među) termini ne prelete preko granica tipa double. To takođe objašnjava zašto ne računate CDF sabiranjem dugog repa članova gustine: svaki član nosi sopstvenu grešku zaokruživanja, greške se akumuliraju kako se serija odvija, a regularizovana specijalna funkcija u potpunosti zaobilazi sumu tako što umesto toga evaluira brzo konvergirajući red ili verižni razlomak (continued fraction)

Redovi ispod dijagonale, verižni razlomak iznad nje

Rutina nepotpune gama funkcije donosi jednu odluku pre nego što išta izračuna: poredi vrednost x sa izrazom a + 1. Ta granica nije proizvoljna. Razvoj u stepeni red funkcije P(a, x) brzo konvergira kada je x malo u odnosu na a, a sporo, na kraju i beskorisno, kada je x veliko. Verižni razlomak ima suprotan karakter. Zato mehanizam koristi stepeni red za x manje od a + 1, i Lencov (Lentz) verižni razlomak za x veće ili jednako a + 1, i od svake grane se traži da radi samo onaj posao u kojem je dobra

Verižnom razlomku je potrebna jedna zaštita (guard). Lencova metoda radi tako što nosi tekući brojilac i imenilac i invertuje imenilac u svakom koraku, a ako se bilo koji od njih približi nuli, inverzija puca. Rešenje je sićušni prag (floor): kad god međuvrednost padne ispod otprilike 1e-30 po veličini, ona se fiksira (clamp) na 1e-30, što održava rekurentnu vezu konačnom bez narušavanja konvergirane vrednosti. Isti fiksator se iz istog razloga pojavljuje i u verižnom razlomku nepotpune beta funkcije. To je mala konstanta koja radi težak posao, čineći razliku između stabilne evaluacije i deljenja nečim što se ne razlikuje od nule

Gornji rep, Q(a, x), je jednostavno 1 minus P(a, x), i tako se računa kumulativna grana Poisson-ove raspodele: verovatnoća od najviše k događaja sa srednjom vrednošću λ je Q(k + 1, λ). Njegovo usmeravanje kroz gornju nepotpunu gama funkciju umesto sabiranja k + 1 Poisson-ovih članova je, opet, izbor da se evaluira jedan konvergirajući izraz umesto akumuliranja mnogo malih vrednosti

Diskretne mase bez prekoračenja faktorijela

Diskretne raspodele donose drugačiju opasnost. Binomna masa verovatnoće uključuje binomni koeficijent, a koeficijent za izbor dvadeset šest od pedeset dva je ogroman ceo broj. Formirajte ga direktno i brojilac će dovesti do prekoračenja (overflow) tipa double pre deljenja koje bi ga vratilo na razumnu verovatnoću. Mehanizam ga nikada ne formira direktno. On računa faktorijele u logaritamskom prostoru preko log-gama funkcije, sabira i oduzima logaritme, dodaje logaritam verovatnoće uspeha i neuspeha, i vrši stepenovanje (exponentiate) samo jednom na samom kraju

// 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 log-gama funkcija je Lancozova (Lanczos) aproksimacija, tačna na čitavoj pozitivnoj osi i jeftina za evaluaciju. Pošto se svaka velika količina čuva kao njen logaritam sve do konačnog Exp poziva, najveći broj koji rutina ikada stvara u memoriji je sama verovatnoća, koja iznosi najviše jedan. Poisson-ova funkcija mase prati isti recept, sa jednim log-gama članom koji zamenjuje faktorijel u imeniocu. Zatvoreni oblici se tretiraju kao specijalni slučajevi na ivicama, gde je p tačno nula ili jedan, tako da kod nikada ne poziva Ln(0). HotXLS vraća 0.2460938 za BINOM.DIST(5,10,0.5,FALSE) i 0.6766764 za kumulativnu funkciju POISSON.DIST(2,2,TRUE), poklapajući se sa Excel-om kroz sve cifre koje štampa

Inverzne funkcije uokvirivanjem (bracketing) direktne krive

Inverzna funkcija raspodele postavlja suprotno pitanje: za datu verovatnoću, pronađi vrednost čiji je CDF jednak njoj. Samo jedna inverzna funkcija u ovom skupu ima brzu direktnu formulu. NORM.S.INV, inverzna standardna normalna raspodela, koristi Aklamovu (Acklam) racionalnu aproksimaciju, par polinomskih odnosa tačnih otprilike do preciznosti tipa double na čitavom opsegu, podeljen na centralnu regiju i dva repa. To je evaluacija zatvorenog oblika bez iteracije

Ostale inverzne funkcije nemaju takvu formulu, pa ih mehanizam invertuje numerički. On uokviruje odgovor donjom i gornjom granicom izabranim iz nosača (support) raspodele, a zatim vrši bisekciju (bisection): evaluira direktni CDF na srednjoj tački, pomera granicu koja zadržava ciljnu verovatnoću zatvorenom, i ponavlja postupak dok interval ne postane uzak. Za inverznu gama i hi-kvadrat funkciju, okvir počinje od nule i velikodušne gornje procene izgrađene na osnovu oblika i skale, duplirajući gornju granicu ako verovatnoća još uvek nije obuhvaćena. Za inverznu t-raspodelu uokviruju se simetrične granice koje se šire ka spolja; inverzna F-raspodela vrši bisekciju na nenegativnom intervalu. Trošak je nekoliko desetina CDF evaluacija po pozivu, što je neprimetno pri brzini rada sa tabelama, a korist je u tome što je svaka inverzna funkcija tačno onoliko precizna koliko i direktna funkcija koju invertuje. Zato kružni tok (round trip) kao što je CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) vraća 0.7 sa zanemarljivim odstupanjem

Dekadni logaritam koji se krio u repu

Evo baga koji vredi ispričati, jer je to vrsta koja dugo preživljava. Aklamova rutina za inverznu normalnu raspodelu ima tri grane. Široka centralna grana, koja se koristi kad god se verovatnoća nalazi između otprilike 0.025 i 0.975, propušta ulaz kroz polinomski odnos bez ikakvog logaritma u sebi. Dve grane na krajevima (repa), za veoma male ili veoma velike verovatnoće, prvo uzimaju logaritam ulaza, jer se rep ponaša kao kvadratni koren iz minus prirodnog logaritma od p

Rana verzija grane repa uzimala je dekadni logaritam (base-10 logarithm) tamo gde je pripadao prirodni logaritam (natural logarithm). Njih dvoje se razlikuju za konstantni faktor od oko 2.30, tako da su rezultati u repu bili pogrešni za doslednu, priličnu marginu. Pa ipak, funkcija je izgledala u redu u svakoj uobičajenoj proveri, jer uobičajene provere žive u sredini. Poziv NORM.S.INV(0.5) je nula, NORM.S.INV(0.975) je školski primer 1.959964, a oba ta primera prolaze kroz centralni polinom koji uopšte ne poziva logaritam. Greška se pojavljivala tek kada bi verovatnoća prešla u rep, recimo kod NORM.S.INV(0.001), koji mora da vrati -3.0902323, a umesto toga se vraćao sa odstupanjem blizu odnosa prirodnog i dekadnog logaritma. Svaka funkcija koja zavisi od inverzne normalne raspodele u svom repu, uključujući pomoćne funkcije intervala poverenja (confidence-interval), nasledila je isti nagib. Pouka je jednostavna i skupa: funkciji sa razgranatom strukturom potrebne su test tačke unutar svake grane, jer će ispravna uobičajena putanja rado sakriti pokvarenu retku putanju. Popravka je bila promena jednog tokena sa dekadnog logaritma na prirodni, i vrednosti u repu su se uskladile sa Excel-ovim

Znak od x odlučuje o repu t-raspodele

Studentova t kumulativna funkcija nosi suptilnost koju je lako shvatiti pogrešno. Njena vrednost dolazi iz regularizovane nepotpune beta funkcije evaluirane pri df / (df + x²), ali ta beta vrednost je verovatnoća u repu izvan veličine x, a ne kumulativna verovatnoća do x. Simetrični oblik t-raspodele znači da konverzija zavisi od toga sa koje strane nule x pada

// 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

Za x iznad nule, kumulativna verovatnoća je jedan minus polovina simetričnog repa; za x ispod nule, to je polovina tog repa; na nuli, to je tačno jedna polovina. Vratite beta vrednost direktno i prijavićete pogrešnu stranu raspodele, promašivši čitavo telo krive za bilo koje x različito od nule. Varijante sa desnim repom i dva repa grade se na istoj grani, zbog čega se T.DIST.2T(1,1) vraća kao 0.5, a T.DIST(1,1,TRUE) kao 0.75, a inverzna funkcija T.INV vrši bisekciju u odnosu na ovaj ispravljeni CDF tako da se kružni tok zatvara

Ništa od ovoga nije vidljivo iz ćelije, i to je planirani ishod. Vi pišete formulu i čitate broj koji se slaže sa Excel-om. Ako proširujete mehanizam sopstvenom logikom, mehanika registrovanja funkcije je pokrivena u našem pregledu mehanizma formula i prilagođenih funkcija, a način na koji formule dopiru do drugih listova i definisanih opsega pokriven je u članku o definisanim imenima i međulistnim formulama. Sve to se isporučuje unutar HotXLS komponente za tabele za Delphi i C++Builder, zajedno sa API-jima za čitanje, pisanje, grafikone i formatiranje koji su pokriveni na drugim mestima na ovom blogu