Upišite =NORM.DIST(115,100,15,TRUE) u ćeliju i Excel vraća 0.8413447 bez ikakvih ceremonija. Poziv izgleda kao obično pretraživanje. Ali nije. Iza tog jednog broja krije se kumulativna normalna razdioba, integral bez zatvorenog oblika, a iza CHISQ.INV.RT i BETA.DIST nalaze se posebne funkcije koje pažljiva knjižnica mora izračunati, a ne aproksimirati ručno. Proračunska komponenta koja tvrdi da je kompatibilna s Excelom mora reproducirati ove vrijednosti do zadnje znamenke koju Excel prikazuje, što znači reproducirati numeričke metode, a ne samo nazive funkcija.
HotXLS implementira više od pedeset ovih statističkih funkcija, a rad koji ih čini točnima gotovo je u potpunosti nevidljiv iz trake formula. Ovo je obilazak načina na koji ih motor izračunava: zajednička jezgra posebnih funkcija, odluke o granama koje održavaju aritmetiku stabilnom, te jedna pogreška inverzne normale koja se dugo skrivala u repu jer uobičajeni slučaj nikada nije dodirnuo neispravnu liniju.
Jedan poziv tablice, pedeset distribucija iza njega
Funkcije obuhvaćaju obitelji za kojima poseže statistička knjiga. Tu je normalna obitelj, NORM.DIST i NORM.S.DIST s njihovim inverzima; obitelj gama i hi-kvadrat, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; obitelj beta, BETA.DIST i BETA.INV; distribucije uzorkovanja T.DIST, T.DIST.2T, F.DIST i F.INV; diskretni par BINOM.DIST i POISSON.DIST; te pomoćnici za zaključivanje poput CONFIDENCE.T i CONFIDENCE.NORM. Iz perspektive pozivatelja svaka je jedna formula. Postavite ulaze u ćelije, zatražite od tablice procjenu 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 knjizi proračunske tablice prevodi i procjenjuje ad-hoc formulu u odnosu na aktivni list i vraća Variant. Jedan detalj zbunjuje ljude pri prvom pokušaju: parser formula iza Calculate uzima točka-zarez kao separator argumenata, tako da je to =SUM(A1;B1), a ne =SUM(A1,B1). Pohranjene formule ćelija zadržavaju Excel-standardni zarez. Isti evaluator otprema svaku statističku funkciju u nastavku, pa čim jedna od njih proradi u Calculate, ostale prate isti put.
Dvije funkcije na kojima se sve ostalo gradi
Većina kumulativnih razdioba u ovom skupu ne izračunava se zbrajanjem ili integriranjem vlastitih definicija. One se računaju iz dviju posebnih funkcija: regularizirane donje nepotpune gama funkcije, napisane kao P(a, x), i regularizirane nepotpune beta funkcije, napisane kao Ix(a, b). Interno su to pomoćnici na koje se dispečeri oslanjaju, a lanac je kratak. Hi-kvadrat CDF je gama CDF s oblikom df/2 i skalom 2. Gama CDF je izravno P(a, x). Kumulativne funkcije t, F i binomna sve su vrijednosti regularizirane nepotpune beta funkcije pri ispravnim argumentima. Poissonov CDF je gornja nepotpuna gama Q. Implementirajte dobro gama i beta funkcije i desetak distribucija naslijedit će njihovu točnost besplatno.
Riječ "regularizirano" je cijela poanta. Sirova nepotpuna gama raste poput faktorijela, a sirovi beta integral može doživjeti podljev ili preljev puno prije nego što to učini sam odgovor. Regularizirani oblici podijeljeni su s potpunom gama ili beta funkcijom, tako da žive u potpunosti u intervalu od nula do jedan, što je točno raspon koji vjerojatnost zauzima. Ta normalizacija je ono što omogućuje istoj rutini da služi hi-kvadratu s dva stupnja slobode i onom s dvjesto bez da se srednji članovi preliju izvan granica dvostruke preciznosti (double). To također objašnjava zašto ne izračunavate CDF zbrajanjem dugog repa članova gustoće: svaki član nosi vlastitu pogrešku zaokruživanja, pogreške se akumuliraju kako se serija odvija, a regularizirana posebna funkcija u potpunosti zaobilazi zbroj procjenjujući umjesto toga brzo konvergirajući red ili verižni razlomak.
Redovi ispod dijagonale, verižni razlomak iznad nje
Rutina nepotpune game donosi jednu odluku prije nego što bilo što izračuna: uspoređuje x s a + 1. Ta granica nije proizvoljna. Razvoj P(a, x) u potencijski red konvergira brzo kada je x malen u odnosu na a, a polako, na kraju i beskorisno, kada je x velik. Verižni razlomak ima suprotan karakter. Tako motor koristi potencijski red za x ispod a + 1 i Lentzov verižni razlomak za x na ili iznad a + 1, te se od svake grane traži da obavi samo onaj posao u kojem je dobra.
Verižni razlomak treba jednu zaštitu. Lentzova metoda radi tako da nosi tekući brojnik i nazivnik i invertira nazivnik na svakom koraku, a ako se bilo koji približi nuli, inverzija propada. Rješenje je sićušni donji prag (floor): kad god srednji član padne ispod otprilike 1e-30 po veličini, steže se na 1e-30, što održava rekurziju konačnom bez ometanja konvergirane vrijednosti. Isti stisak se pojavljuje u verižnom razlomku nepotpune bete iz istog razloga. To je mala konstanta koja radi težak posao, razlika između stabilne evaluacije i dijeljenja s nečim što se ne razlikuje od nule.
Gornji rep, Q(a, x), je jednostavno 1 minus P(a, x), i to je način na koji se izračunava Poissonova kumulativna grana: vjerojatnost od najviše k događaja sa srednjom vrijednošću λ je Q(k + 1, λ). Usmjeravanje kroz gornju nepotpunu gamu umjesto zbrajanja k + 1 Poissonovih članova je, opet, izbor da se procijeni jedan konvergirajući izraz umjesto akumuliranja mnogih malih.
Diskretne mase bez prelijevanja faktorijela
Diskretne distribucije donose drugačiju opasnost. Binomna masa vjerojatnosti uključuje binomni koeficijent, a koeficijent za pedeset-dva-odaberi-dvadeset-šest je ogroman cijeli broj. Formirajte ga izravno i brojnik će preliti double prije dijeljenja koje bi ga vratilo na razumnu vjerojatnost. Motor ga nikada ne formira. On izračunava faktorijele u logaritamskom prostoru putem log-gama funkcije, zbraja i oduzima logaritme, pribraja logaritme vjerojatnosti uspjeha i neuspjeha te potencira 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 Lanczoseva aproksimacija, točna na cijeloj pozitivnoj osi i jeftina za evaluaciju. Budući da se svaka velika količina drži kao svoj logaritam do konačnog Exp, najveći broj koji rutina ikada stvara je sama vjerojatnost, koja je najviše jedan. Poissonova funkcija mase slijedi isti recept, s jednim log-gama članom koji mijenja faktorijel u nazivniku. Zatvoreni oblici se posebno tretiraju na rubovima, gdje je p toč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 POISSON.DIST(2,2,TRUE), podudarajući se s Excelom u svim prikazanim znamenkama.
Inverzne vrijednosti uokvirivanjem izravne krivulje
Funkcija inverzne razdiobe postavlja suprotno pitanje: za zadani vjerojatnost pronaći vrijednost čiji je CDF jednak njoj. Samo jedan inverz u ovom skupu ima brzu izravnu formulu. NORM.S.INV, inverzna standardna normalna, koristi Acklamovu racionalnu aproksimaciju, par omjera polinoma točnih na otprilike preciznost double-a u cijelom rasponu, podijeljen na središnje područje i dva repa. To je evaluacija zatvorenog oblika bez iteracije.
Ostali inverzi nemaju takvu formulu, pa ih motor invertira numerički. On uokviruje odgovor donjom i gornjom granicom odabranom iz podrške distribucije, a zatim vrši bisekciju: procjenjuje izravni CDF na središnjoj točki, pomiče granicu koja drži ciljnu vjerojatnost zatvorenom i ponavlja postupak dok interval ne postane uzak. Za inverze game i hi-kvadrata, okvir počinje od nule i velikodušne gornje procjene izgrađene na temelju oblika i skale, udvostručujući gornju granicu ako vjerojatnost još nije zatvorena. Inverz t-distribucije uokviruje simetrične granice koje se šire prema van; inverz F-distribucije vrši bisekciju na nenegativnom intervalu. Cijena je nekoliko desetaka evaluacija CDF-a po pozivu, što je nevidljivo pri brzini proračunske tablice, a prednost je što je svaki inverz točno onoliko točan koliko i izravna funkcija koju invertira. Zato povratno putovanje poput CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) returns 0.7 to within a hair.
Logaritam s bazom 10 koji se skrivao u repu
Evo buga koji vrijedi ispričati, jer je to vrsta koja dugo preživljava. Acklamova inverzno-normalna rutina ima tri grane. Široka središnja grana, koja se koristi kad god vjerojatnost leži između otprilike 0.025 i 0.975, provlači ulaz kroz omjer polinoma bez ikakvog logaritma u sebi. Dvije grane repa, za vrlo male ili vrlo velike vjerojatnosti, svaka najprije uzimaju logaritam ulaza, jer se rep ponaša kao kvadratni korijen minus prirodnog logaritma od p.
Rana verzija grane repa uzimala je logaritam s bazom 10 tamo gdje je pripadao prirodni logaritam. Njih dva se razlikuju za konstantni faktor od oko 2.30, pa su rezultati u repu bili pogrešni za dosljednu, znatnu marginu. Pa ipak, funkcija je izgledala dobro u svakoj usputnoj provjeri, jer usputne provjere žive u sredini. NORM.S.INV(0.5) je nula, NORM.S.INV(0.975) je školski 1.959964, i oba prolaze kroz središnji polinom koji uopće ne poziva logaritam. Pogreška se pojavila tek kada bi vjerojatnost prešla u rep, recimo NORM.S.INV(0.001), koji mora vratiti -3.0902323, a umjesto toga se vratio s odstupanjem zbog omjera prirodnog logaritma i onog s bazom 10. Svaka funkcija koja ovisi o inverznoj normali u svom repu, uključujući pomoćnike za intervale pouzdanosti, naslijedila je istu iskrivljenost. Lekcija je jednostavna i skupa: funkcija s razgranatom strukturom treba testne točke unutar svake grane, jer će ispravna uobičajena putanja veselo prikriti pokvarenu rijetku. Rješenje je bila promjena jednog tokena iz logaritma s bazom 10 u prirodni logaritam, i vrijednosti u repu su se uskladile s Excelovima.
Predznak od x odlučuje o repu t-razdiobe
Studentova kumulativna t-funkcija nosi suptilnost koju je lako shvatiti obrnuto. Njezina vrijednost dolazi iz regularizirane nepotpune bete procijenjene na df / (df + x²), no ta beta vrijednost je vjerojatnost u repu izvan veličine x, a ne kumulativna vjerojatnost do x. Simetrični oblik t-razdiobe znači da konverzija ovisi o tome na koju stranu od nule pada 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 meanZa x iznad nule kumulativna vjerojatnost je jedan minus polovica simetričnog repa; za x ispod nule to je polovica tog repa; na nuli je točno jedna polovica. Vratite beta vrijednost izravno i prijavit ćete pogrešnu stranu distribucije, s odstupanjem za cijelo tijelo krivulje za bilo koji nonzero x. Varijante desnog repa i dvaju repova 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 inverz T.INV vrši bisekciju u odnosu na ovaj ispravljeni CDF tako da se krug zatvara.
Ništa od ovoga nije vidljivo iz ćelije, i to je željeni ishod. Upišete formulu i pročitate broj koji se slaže s Excelom. Ako proširujete motor vlastitom logikom, mehanika registracije funkcije pokrivena je u našem vodiču kroz motor formula i prilagođene funkcije, a način na koji formule dosežu preko listova i imenovanih raspona pokriven je u članku o definiranim nazivima i formulama između listova. Sve se to isporučuje unutar proračunske komponente HotXLS za Delphi i C++Builder, uz API-je za čitanje, pisanje, grafikone i oblikovanje koji su pokriveni drugdje na ovom blogu.