Excelin tilastolliset funktiot Delphissä: NORM, CHISQ, BETA

Kirjoita soluun =NORM.DIST(115,100,15,TRUE), ja Excel palauttaa luvun 0.8413447 sen kummemmitta seremonioitta. Kutsu vaikuttaa hakutoiminnolta. Se ei ole sitä. Tämän yhden luvun takana on kumulatiivinen normaalijakauma, integraali ilman suljettua muotoa, ja funktioiden CHISQ.INV.RT ja BETA.DIST takana on erityisfunktioita, jotka huolellisen kirjaston on arvioitava, ei arvioitava käsin. Taulukkolaskentakomponentin, joka väittää olevansa Excel-yhteensopiva, on tuotettava nämä arvot viimeistä Excelin näyttämää numeroa myöten, mikä tarkoittaa numeeristen menetelmien, ei vain funktioiden nimien, toistamista.

HotXLS toteuttaa yli viisikymmentä näistä tilastollisista funktioista, ja työ, joka tekee niistä oikeita, on lähes täysin näkymätöntä kaavariviltä. Tämä on esittely siitä, miten moottori laskee ne: jaettu erityisfunktioperusta, haarautumispäätökset, jotka pitävät aritmetiikan vakaana, ja yksi käänteinen normaalijakaumavirhe, joka piileskeli häntäosassa pitkään, koska yleinen tapaus ei koskaan koskettanut rikkinäistä koodiriviä.

Yksi taulukkokutsu, viisikymmentä jakaumaa sen takana

Funktiot kattavat ne ryhmät, joita tilastollinen työkirja tarvitsee. Siellä on normaalijakaumaperhe, NORM.DIST ja NORM.S.DIST käänteisfunktioineen; gamma- ja khiin neliö -perhe, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; betajakaumaperhe, BETA.DIST ja BETA.INV; otanta-jakaumat T.DIST, T.DIST.2T, F.DIST ja F.INV; diskreetti pari BINOM.DIST ja POISSON.DIST; sekä päättelyapulaiset, kuten CONFIDENCE.T ja CONFIDENCE.NORM. Kutsujan näkökulmasta kukin on yksittäinen kaava. Asetat syötteet soluihin, pyydät työkirjaa arvioimaan ja luet tuloksen.

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;

Työkirjan Calculate-metodi kääntää ja arvioi ad-hoc-kaavan reaaliaikaista taulukkoa vasten ja palauttaa Variant-arvon. Yksi yksityiskohta voi aiheuttaa päänvaivaa ensimmäisellä kerralla: Calculate-metodin takana oleva kaavajäsentäjä käyttää puolipistettä argumenttien erottimena, joten se on muotoa =SUM(A1;B1), ei =SUM(A1,B1). Tallennetut solukaavat säilyttävät Excel-standardin mukaisen pilkun. Sama arvioija lähettää jokaisen alla olevan tilastollisen funktion, joten kun yksi näistä toimii Calculate-metodissa, muut seuraavat samaa polkua.

Kaksi funktiota, joille kaikki muu on rakennettu

Suurinta osaa tämän joukon kumulatiivisista jakaumista ei lasketa laskemalla yhteen tai integroimalla niiden omia määritelmiä. Ne lasketaan kahdesta erityisfunktiosta: regularisoidusta alemmasta epätäydellisestä gammafunktiosta, joka kirjoitetaan P(a, x), ja regularisoidusta epätäydellisestä betafunktiosta, joka kirjoitetaan I_x(a, b). Sisäisesti nämä ovat apureita, joihin jakelijat tukeutuvat, ja ketju on lyhyt. Khiin neliön CDF on gamma-jakauman CDF, jonka muoto on df/2 ja asteikko 2. Gammakauman CDF on suoraan P(a, x). T-, F- ja binomijakauman kumulatiiviset funktiot ovat kaikki regularisoidun epätäydellisen betafunktion arvoja oikeilla argumenteilla. Poissonin CDF on ylempi epätäydellinen gammafunktio Q. Toteuta gamma- ja betafunktiot hyvin, niin tusina jakaumaa perii niiden tarkkuuden ilmaiseksi.

Sana "regularisoitu" on koko homman ydin. Raaka epätäydellinen gamma kasvaa kuin kertoma ja raaka betaintegraali voi alivuotaa tai ylivuotaa kauan ennen kuin vastaus tekee niin. Regularisoidut muodot jaetaan kokonaisella gammalla tai betalla, joten ne elävät täysin välillä nollasta yhteen, mikä on juuri se alue, jonka todennäköisyys täyttää. Kyseinen normalisointi on se, mikä antaa saman rutiinin palvella khiin neliötä kahdella vapaudenasteella ja sellaista, jolla on kaksisataa, ilman että välitermit karkaavat doublen ulkopuolelle. Se myös selittää, miksi et laske CDF:ää laskemalla yhteen pitkää tiheystermien häntää: jokainen termi kantaa oman pyöristysvirheensä, virheet kertyvät sarjan käydessä, ja regularisoitu erityisfunktio ohittaa summan kokonaan arvioimalla sen sijaan nopeasti konvergoituvan sarjan tai ketjumurtoluvun.

Sarja diagonaalin alapuolella, ketjumurtoluku sen yläpuolella

Epätäydellisen gamman rutiini tekee yhden päätöksen ennen kuin se laskee mitään: se vertaa x-arvoa arvoon a + 1. Kyseinen raja ei ole mielivaltainen. P(a, x):n potenssisarjan laajennus konvergoituu nopeasti, kun x on pieni suhteessa a-arvoon, ja hitaasti, lopulta hyödyttömästi, kun x on suuri. Ketjumurtoluvulla on päinvastainen luonne. Joten moottori käyttää potenssisarjaa, kun x on alle a + 1, ja Lentzin ketjumurtolukua, kun x on vähintään a + 1, ja kumpaakin haaraa pyydetään tekemään vain sitä työtä, jossa se on hyvä.

Ketjumurtoluku tarvitsee yhden suojan. Lentzin menetelmä toimii kuljettamalla juoksevaa osoittajaa ja nimittäjää ja kääntämällä nimittäjän jokaisella vaiheella, ja jos jompikumpi lähestyy nollaa, kääntäminen räjähtää. Korjaus on pieni lattia: aina kun välitermi putoaa suuruusluokaltaan alle arvon 1e-30, se puristetaan arvoon 1e-30, mikä pitää toistuvuuden rajallisena häiritsemättä konvergoitunutta arvoa. Sama puristus esiintyy epätäydellisen betan ketjumurtoluvussa samasta syystä. Se on pieni vakio, joka tekee kantavaa työtä – ero vakaan arvioinnin ja sellaisen nollalla jakamisen välillä, jota ei voi erottaa nollasta.

Ylempi häntä, Q(a, x), on yksinkertaisesti 1 miinus P(a, x), ja näin Poissonin kumulatiivinen haara lasketaan: enintään k tapahtuman todennäköisyys keskiarvolla λ on Q(k + 1, λ). Sen reitittäminen ylemmän epätäydellisen gamman kautta k-kappaleen Poisson-termin summaamisen sijaan on jälleen valinta arvioida yksi konvergoituva lauseke sen sijaan, että kerrytettäisiin monia pieniä.

Diskreetit massat ilman kertoman ylivuotoa

Diskreetit jakaumat nostavat esiin toisen vaaran. Binomijakauman todennäköisyysmassa sisältää binomialikertoimen, ja kerroin viidellekymmenellekahdelle-valitse-kaksikymmentäkuusi on valtava kokonaisluku. Muodosta se suoraan, niin osoittaja ylivuotaa doublen ennen jakoa, joka toisi sen takaisin järkevään todennäköisyyteen. Moottori ei koskaan muodosta sitä suoraan. Se laskee kertomat logaritmitilassa log-gamma-funktion kautta, laskee yhteen ja vähentää logaritmit, lisää mukaan onnistumis- ja epäonnistumistodennäköisyyksien logaritmit ja suorittaa eksponentiaation vasta aivan lopussa.

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

Itse log-gamma-funktio on Lanczosin approksimaatio, joka on tarkka koko positiivisella akselilla ja halpa arvioida. Koska jokainen suuri määrä pidetään logaritmina lopulliseen Exp-kutsuun asti, suurin luku, jonka rutiini koskaan realisoi, on itse todennäköisyys, joka on enintään yksi. Poissonin massafunktio noudattaa samaa kaavaa nimittäjän kertoman korvaavalla yksittäisellä log-gamma-termillä. Suljetut muodot käsitellään erikseen reunoilla, joissa p on tarkalleen nolla tai yksi, joten koodi ei koskaan kutsu Ln(0)-funktiota. HotXLS palauttaa arvon 0.2460938 kaavalle BINOM.DIST(5,10,0.5,FALSE) ja arvon 0.6766764 kumulatiiviselle kaavalle POISSON.DIST(2,2,TRUE) vastaten Exceliä sen tulostamien numeroiden kautta.

Käänteisarvot sulkemalla eteenpäin suuntautuva käyrä

Käänteinen jakaumafunktio kysyy päinvastaista kysymystä: etsi annetulle todennäköisyydelle arvo, jonka CDF vastaa sitä. Vain yhdellä tämän joukon käänteisfunktiolla on nopea suora kaava. Standardinormaalijakauman käänteisfunktio NORM.S.INV käyttää Acklamin rationaalista approksimaatiota, joka on polynomisuhteiden pari, joka on tarkka noin doublen tarkkuudelle koko alueella jaettuna keskialueeseen ja kahteen häntään. Se on suljetun muodon arviointi ilman iteraatiota.

Muilla käänteisfunktioilla ei ole tällaista kaavaa, joten moottori kääntää ne numeerisesti. Se rajaa vastauksen jakauman tuesta valitulla alarajalla ja ylärajalla ja puolittaa: arvioi eteenpäin suuntautuvan CDF:n keskipisteessä, siirtää sitä rajaa, joka pitää kohdetodennäköisyyden suljettuna, ja toistaa, kunnes väli on kapea. Gammakauman ja khiin neliön käänteisarvoille rajaus alkaa nollasta ja anteliaasta ylärajasta, joka on rakennettu muodosta ja asteikosta, tuplaten ylärajan, jos todennäköisyys ei vielä ole rajattu. T-käänteisarvo rajaa symmetriset rajat, jotka laajenevat ulospäin; F-käänteisarvo puolittaa ei-negatiivisella välillä. Kustannus on muutaman kymmenen CDF-arviointia kutsua kohden, mikä on näkymätön taulukkolaskentanopeudella, ja etuna on, että jokainen käänteisfunktio on täsmälleen yhtä tarkka kuin eteenpäin suuntautuva funktio, jota se kääntää. Siksi edestakainen matka, kuten CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE), palauttaa arvon 0.7 hiuksenhienolla tarkkuudella.

Kymmenkantainen logaritmi, joka piiloutui häntään

Tässä on mainitsemisen arvoinen virhe, koska se on sellainen, joka selviää pitkään. Acklamin käänteisen normaalijakauman rutiinilla on kolme haaraa. Laaja keskushara, jota käytetään aina, kun todennäköisyys on välillä 0.025–0.975, ajaa syötteen polynomisuhteen läpi ilman logaritmia missään kohdassa koodia. Kaksi häntähaaraa hyvin pienille tai erittäin suurille todennäköisyyksille ottavat kumpikin ensin logaritmin syötteestä, koska häntä käyttäytyy kuin neliöjuuri miinus p:n luonnollisesta logaritmista.

Häntähaaran varhainen versio otti kymmenkantaisen logaritmin sinne, minne luonnollinen logaritmi kuului. Nämä kaksi eroavat toisistaan vakiokertoimella, joka on noin 2.30, joten hännän tulokset olivat väärin tasaisella, huomattavalla marginaalilla. Ja silti funktio näytti hyvältä jokaisessa satunnaisessa tarkistuksessa, koska satunnaiset tarkistukset elävät keskellä. NORM.S.INV(0.5) on nolla, NORM.S.INV(0.975) on oppikirjamainen 1.959964, ja molemmat kulkevat keskimmäisen polynomin läpi, joka ei koskaan kutsu logaritmia lainkaan. Virhe tuli näkyviin vasta, kun todennäköisyys siirtyi häntään, esimerkiksi NORM.S.INV(0.001), jonka on palautettava -3.0902323 ja joka sen sijaan palautti virheellisen arvon luonnollisen ja kymmenkantaisen logaritmin suhteen vuoksi. Mikä tahansa funktio, joka riippuu käänteisestä normaalijakaumasta hännässään, mukaan lukien luottamusväliapulaiset, peri saman vinouman. Opetus on arkinen ja kallis: funktio, jolla on haararakenne, tarvitsee testipisteitä jokaisen haaran sisälle, koska oikea yleinen polku peittää iloisesti rikkinäisen harvinaisen polun. Korjaus oli yhden merkin muutos kymmenkantaisesta logaritmista luonnolliseen logaritmiin, ja hännän arvot napsahtivat Excelin arvoihin.

x-arvon etumerkki ratkaisee t-jakauman hännän

Studentin t-jakauman kumulatiivinen funktio kantaa hienovaraisuutta, joka on helppo saada väärinpäin. Sen arvo saadaan regularisoidusta epätäydellisestä betasta arvioituna kohdassa df / (df + x²), mutta tämä beta-arvo on todennäköisyys hännässä x-arvon suuruuden tuolla puolen, ei kumulatiivinen todennäköisyys x-arvoon asti. T-jakauman symmetrinen muoto tarkoittaa, että muunnos riippuu siitä, kummalle puolelle nollaa x putoaa.

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

Nollan yläpuolella oleville x-arvoille kumulatiivinen todennäköisyys on yksi miinus puolet symmetrisestä hännästä; nollan alapuolella oleville x-arvoille se on puolet tuosta hännästä; tarkalleen nollassa se on täsmälleen puolet. Palauta beta-arvo suoraan, niin ilmoitat väärän puolen jakaumasta, mikä heittää koko käyrän rungolla minkä tahansa nollasta poikkeavan x-arvon kohdalla. Oikea häntä- ja kaksihäntävariantit rakentuvat samalle haaralle, minkä vuoksi T.DIST.2T(1,1) palauttaa arvon 0.5 ja T.DIST(1,1,TRUE) arvon 0.75, ja käänteinen T.INV suorittaa puolitushaun tätä korjattua CDF:ää vasten, jotta edestakainen matka sulkeutuu.

Mikään tästä ei näy soluun, ja se on tarkoitettu lopputulos. Kirjoitat kaavan ja luet luvun, joka vastaa Exceliä. Jos laajennat moottoria omalla logiikallasi, funktion rekisteröinnin mekaniikka käsitellään ohjeessamme kaavamoottorista ja mukautetuista funktioista, ja tapa, jolla kaavat tavoittavat taulukoita ja nimettyjä alueita, käsitellään artikkelissa määritellyistä nimistä ja ristiintaulukkokaavoista. Kaikki tämä toimitetaan Delphin ja C++Builderin HotXLS-taulukkolaskentakomponentin sisällä luonti-, kirjoitus-, kaavio- ja muotoilu-API:en rinnalla, joita käsitellään muualla tässä blogissa.