Tehnični članak

Statistične funkcije Excel v Delphi: NORM, CHISQ, BETA

V celico vpišite =NORM.DIST(115,100,15,TRUE) in Excel brez posebnega slovesa vrne 0.8413447. Klic se bere kot iskanje. Vendar ni tako. Za to številko stoji kumulativna normalna porazdelitev, integral brez zaključene oblike, za CHISQ.INV.RT in BETA.DIST pa stojijo posebne funkcije, ki jih mora natančna knjižnica ovrednotiti in ne le oceniti na roko. Komponenta delovnega lista, ki trdi, da je združljiva z Excelom, mora poustvariti te vrednosti do zadnje številke, ki jo prikazuje Excel, kar pomeni poustvarjanje številskih metod in ne le imen funkcij

HotXLS izvaja več kot petdeset teh statističnih funkcij, delo, zaradi katerega so pravilne, pa je skoraj v celoti nevidno iz vrstice za formule. To je pregled, kako jih motor izračuna: skupno jedro posebnih funkcij, vejitve, ki ohranjajo aritmetiko stabilno, ter ena napaka inverzne normalne porazdelitve, ki se je dolgo skrivala v repu, ker pogost primer nikoli ni prizadel pokvarjene vrstice

En klic delovnega lista, za njim pa petdeset porazdelitev

Funkcije obsegajo družine, po katerih posega statistični delovni zvezek. Na voljo je normalna družina, NORM.DIST in NORM.S.DIST z njihovimi inverzi; družina gama in hi-kvadrat, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; družina beta, BETA.DIST in BETA.INV; vzorčne porazdelitve T.DIST, T.DIST.2T, F.DIST in F.INV; diskretni par BINOM.DIST in POISSON.DIST; ter pomočniki za sklepanje, kot sta CONFIDENCE.T in CONFIDENCE.NORM. Z mesta klicatelja je vsaka od njih posamezna formula. Vhode nastavite v celicah, zahtevate od delovnega zvezka, da jih ovrednoti, in preberete 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 delovnem zvezku prevede in ovrednoti ad-hoc formulo glede na živi list ter vrne Variant. Ena podrobnost ljudi zmoti ob prvem poskusu: parser formul za Calculate vzame podpičje kot ločilo argumentov, zato je =SUM(A1;B1) in ne =SUM(A1,B1). Shranjene formule v celicah ohranjajo standardno Excelovo vejico. Isti vrednotilnik odpremi vsako spodnjo statistično funkcijo, tako da ko ena od teh deluje v Calculate, ostale sledijo isti poti

Dve funkciji, na katerih sloni vse ostalo

Večina kumulativnih porazdelitev v tem naboru se ne izračuna s seštevanjem ali integriranjem lastnih definicij. Izračunajo se iz dveh posebnih funkcij: regularizirane spodnje nepopolne funkcije gama, zapisane kot P(a, x), in regularizirane nepopolne funkcije beta, zapisane kot Ix(a, b). Interno sta to pomočnika, na katera se opirajo odpremniki, veriga pa je kratka. CDF hi-kvadrat je CDF gama z obliko df/2 in merilom 2. CDF gama je neposredno P(a, x). Kumulativne funkcije t, F in binomska kumulativna funkcija so vse vrednosti regularizirane nepopolne funkcije beta pri pravih argumentih. Poissonov CDF je zgornja nepopolna gama Q. Dobro izvedite funkciji gama in beta in ducat porazdelitev bo brezplačno podedoval njihovo natančnost

Beseda "regularizirano" je bistvo vsega. Surova nepopolna gama raste kot faktorska vrednost, surov integral beta pa lahko doživi podpretok ali prekoračitev veliko prej kot sam odgovor. Regularizirane oblike so deljene s celotno gamo ali beto, zato živijo v celoti v intervalu od nič do ena, kaj je natanko obseg, ki ga zaseda verjetnost. Ta normalizacija omogoča, da ista rutina služi hi-kvadratu z dvema stopnjama prostosti in tistemu z dvesto stopnjami, ne da bi vmesni členi pobegnili izven obsega double. Pojasnjuje tudi, zakaj CDF-ja ne izračunate s seštevanjem dolgega repa gostotnih členov: vsak člen prinaša lastno zaokrožitveno napako, napake se med tekom vrste kopičijo, regularizirana posebna funkcija pa se seštevku v celoti izogne tako, da namesto tega ovrednoti hitro konvergentno vrsto ali verižni ulomek

Vrsta pod diagonalo, verižni ulomek nad njo

Vrsta pod diagonalo, verižni ulomek nad njo

Rutina za nepopolno gamo sprejme eno odločitev, preden sploh kaj izračuna: primerja x z a + 1. Ta meja ni poljubna. Razvoj v potenčno vrsto za P(a, x) hitro konvergira, ko je x majhen glede na a, in počasi, na koncu neuporabno, ko je x velik. Verižni ulomek ima nasprotne lastnosti. Zato motor uporablja potenčno vrsto za x pod a + 1 in Lentzov verižni ulomek za x pri ali nad a + 1, od vsake veje pa se zahteva le delo, v katerem je dobra

Verižni ulomek potrebuje eno zaščito. Lentzova metoda deluje tako, da prenaša tekoči števec in imenovalec ter na vsakem koraku invertira imenovalec; če se kateri od njiju približa ničli, inverzija eksplodira. Popravek je drobno dno: kadar koli vmesni člen pade pod približno 1e-30 po velikosti, se vpnne na 1e-30, kar ohranja rekurzivnost končno, ne da bi motilo konvergirano vrednost. Enako vpenjanje se iz istega razloga pojavi v verižnem ulomku nepopolne bete. To je majhna konstanta, ki opravlja nosilno delo, razlika med stabilnim vrednotenjem in deljenjem z nečim, kar je nerazpoznavno od ničle

Zgornji rep, Q(a, x), je preprosto 1 minus P(a, x) in tako se izračuna Poissonova kumulativna veja: verjetnost največ k dogodkov s srednjo vrednostjo λ je Q(k + 1, λ). Usmerjanje skozi zgornjo nepopolno gamo namesto seštevanja k + 1 Poissonovih členov je znova izbira ovrednotenja enega konvergentnega izraza namesto kopičenja več majhnih

Diskretne mase brez prekoračitve faktoriala

Diskretne porazdelitve prinašajo drugačno nevarnost. Binomska masa verjetnosti vključuje binomski koeficient, koeficient za petindvajset-izmed-dvainpetdeset pa je ogromno celo število. Če ga ustvarite neposredno, števec prekorači double še pred deljenjem, ki bi ga vrnilo na razumno verjetnost. Motor ga nikoli ne ustvari. Faktoriale izračuna v logaritemskem prostoru prek funkcije log-gama, sešteje in odšteje logaritme, zloži logaritme verjetnosti uspeha in neuspeha ter eksponentira le enkrat na samem koncu

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

Funkcija log-gama sama je Lanczosev približek, natančna čez celotno pozitivno os in poceni za vrednotenje. Ker se vsaka velika količina hrani kot njen logaritem do končnega Exp, je največje število, ki ga rutina sploh materializira, sama verjetnost, ki je največ ena. Poissonova masna funkcija sledi istemu receptu, z enim členom log-gama, ki nadomešča faktorial v imenovalcu. Zaključene oblike so posebej obravnavane na robovih, kjer je p natanko nič ali ena, tako da koda nikoli ne pokliče Ln(0). HotXLS vrne 0.2460938 za BINOM.DIST(5,10,0.5,FALSE) in 0.6766764 za kumulativni POISSON.DIST(2,2,TRUE), kar se ujema z Excelom skozi vse izpisane številke

Inverzi z oklepanjem usmerjene krivulje

Drugi inverzi nimajo takšne formule, zato jih motor invertira številsko. Odgovor obda z spodnjo in zgornjo mejo, izbranima iz podpore porazdelitve, nato pa bisektira: ovrednoti usmerjeni CDF na srednji točki, premakne tisto mejo, ki ohranja ciljno verjetnost zaprto, in ponavlja, dokler ni interval ozek. Za inverza gama in hi-kvadrat se oklepaj začne pri ničli in radodarni zgornji oceni, zgrajeni iz oblike in merila, pri čemer podvoji zgornjo mejo, če verjetnost še ni zaprta. Inverz t oklepa simetrični meji, ki se širita navzven; inverz F bisektira na nenegativnem intervalu. Strošek je nekaj ducat vrednotenj CDF na klic, kar je pri hitrosti delovnega lista neopazno, prednost pa je, da je vsak inverz natanko tako natančen kot usmerjena funkcija, ki jo invertira. Zato round-trip, kot je CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE), vrne 0.7 do lasu natančno

Logaritem z osnovo 10, ki se je skrival v repu

Zgodnja različica veje repa je vzela logaritem z osnovo 10 tam, kjer je bil mesto naravnemu logaritmu. Razlikujeta se za konstantni faktor približno 2.30, so so bili rezultati v repu napačni za skladen, precejšen odmik. Pa vendar je bila funkcija ob vsakem priložnostnem preverjanju videti v redu, ker priložnostna preverjanja živijo v sredini. NORM.S.INV(0.5) je nič, NORM.S.INV(0.975) je učbeniški 1.959964 in oba tečeta skozi osrednji polinom, ki sploh ne kliče logaritma. Napaka se je pojavila šele, ko je verjetnost prešla v rep, na primer NORM.S.INV(0.001), ki mora vrniti -3.0902323, namesto tega pa se je vrnil z odmikom razmerja naravnega logaritma in logaritma z osnovo 10. Vsaka funkcija, ki je v svojem repu odvisna od inverzne normale, vključno s pomočniki za intervale zaupanja, je podedovala isto izkrivljenost. Nauk je vsakodneven in drag: funkcija z vejitveno strukturo potrebuje testne točke znotraj vsake veje, saj bo pravilna skupna pot z veseljem prikrila pokvarjeno redko. Popravek je bil sprememba enega žetona iz logaritma z osnovo 10 v naravni logaritem, vrednosti v repu pa so se zaskočile z Excelovimi

Znak x določa rep porazdelitve t

Studentova kumulativna funkcija t prinaša prefinjenost, ki jo je enostavno obrniti narobe. Njena vrednost izhaja iz regularizirane nepopolne bete, ovrednotene pri df / (df + x²), vendar je ta vrednost bete verjetnost v repu zunaj velikosti x in ne kumulativna verjetnost do x. Simetrična oblika porazdelitve t pomeni, da je pretvorba odvisna od tega, na kateri strani ničle pade x

Za x nad ničlo je kumulativna verjetnost ena minus polovica simetričnega repa; za x pod ničlo je to polovica tega repa; pri ničli je to natanko ena polovica. Vrnite vrednost bete neposredno in poročali boste o napačni strani porazdelitve, zamaknjeni za celotno telo krivulje za kateri koli neničelni x. Različici z desnim repom in dvema repoma gradita na isti veji, zato se T.DIST.2T(1,1) vrne kot 0.5 in T.DIST(1,1,TRUE) kot 0.75, inverz T.INV bisektira proti temu popravljenemu CDF-ju, tako da se round-trip zapre

Nič od tega ni vidno iz celice in to je načrtovani rezultat. Napišete formulo in preberete številko, ki se ujema z Excelom. Odločitev o razširitvi motorja z lastno logiko je zajeta v našem pregledu motorja formul in uporabniških funkcij, način, kako formule segajo čez liste in imenovana območja, pa je zajet v članku o definiranih imenih in formulah čez liste. Vse to se dobavlja znotraj komponente delovnega lista HotXLS Component za Delphi in C++Builder, skupaj z API-ji za branje, pisanje, grafikone in oblikovanje, ki so obravnavani drugje na tem blogu