Afsnit 2.5: Poissonmodellen

Jeg slutter dette kapitel af med at sige noget om poissonfordelingen, som I også kender fra sandsynlighedsdelen af calculus. Jeg vil beskrive, hvordan poissonfordelingen kommer frem som en model for tilfældige ankomster. Dernæst vil jeg estimere parameteren i poissonmodellen og lave et konfidensinterval for denne.
Statistisk Model 2.5.1. (Poissonmodellen)
Lad XX være et stokastisk antal, som kan antage alle mulige heltallige værdier større end eller lig med nul. Poissonmodellen med parameter λ\lambda til beskrivelse af XX skrives som
X∼poisson(λ),λ≥0. X\sim\text{poisson}(\lambda),\enspace \lambda\geq 0.
For poissonmodellen gælder der følgende resultater:
P(X=x)=λxx!e−λ,x=0,1,2,…,E(X)=λ,Var(X)=λ.\begin{aligned} & P(X=x)=\frac{\lambda^x}{x!}e^{-\lambda},\enspace x=0,1,2,\ldots, \\ & E(X)=\lambda,\quad \text{Var}(X)=\lambda. \end{aligned}
Desuden har vi egenskaben, at hvis X1,…,XkX_1,\ldots,X_k er uafhængige og poissonfordelte, så er også summen X∙=X1+⋯+XkX_\bullet=X_1+\cdots+X_k poissonfordelt. Den præcise formulering er, at hvis Xi∼poisson(λi),X_i\sim\text{poisson}(\lambda_i), i=1,…,k,i=1,\ldots,k, og de kk stokastiske variable er uafhængige, så er X∙∼poisson(λ∙)X_\bullet\sim\text{poisson}(\lambda_\bullet), hvor X∙X_\bullet er summen af XiX_i-erne og λ∙\lambda_\bullet er summen af λi\lambda_i-erne. Poissonfordelingen bruges til at beskrive fordelingen af et antal, der for eksempel kommer tilfældigt i tid eller tilfældigt over område. Berømte eksempler er antal opkald til en telefoncentral inden for et bestemt tidsinterval (blev studeret af den danske matematiker Erlang), dødsfald ved hestespark i preussiske regimenter, hvor data består af dødsfald per år i 10 regimenter der følges i 20 år (1875‐1894) (billedet nedenfor er fra bogen The art of taming and educating the horse), og fordeling af V1- og V2-bomber over London under anden verdenskrig.

Figur med sandsynligheder

Følgende kode producerer en figur med poissonsandsynligheder. Kør koden, og prøv at ændre på værdien af Lambda i koden (advarsel: benyt aldrig "lambda" som variabelnavn i python). I python beregnes poissonsandsynligheder med st.poisson.pmf.
Kan du se, at "kurveformen" af fordelingen ser ud til at nærme sig en bestemt form, hvis du prøver at gøre λ\lambda større og større ?

Tilfældige ankomster

Jeg beskriver her en model, der kan forklare, hvorfor mange data kan være poissonfordelte. Modellen skal beskrive situationen med ankomster tilfældigt i tid. Betragt tidsintervallet fra 0 til t.t. Intervallet deles op i nn lige store stykker, hvor nn er stor. Uafhængigt af hinanden kaster vi i hvert af de nn små intervaller en "mønt", hvor sandsynligheden for krone er λtn.\lambda\frac{t}{n}. Hvis mønten viser krone, siger vi, at der er en ankomst i det tilsvarende interval. På denne måde er sandsynligheden for en ankomst i et lille interval proportional med længden af intervallet. Proportionalitetskonstanten er λ,\lambda, som får fortolkningen som rate, det vil sige det forventede antal ankomster per tid. Uafhængigheden mellem de små intervaller modellerer, at ankomsterne er tilfældige i tid.
Ideen er nu at lade nn blive større og større, svarende til en finere og finere inddeling af intervallet [0,t][0,t] i små intervaller. Lad XX være den stokastiske variabel, der angiver antal ankomster i intervallet [0,t].[0,t]. Når nn er valgt, står vi i en binomialsituation og har X∼binom(n,λt/n).X\sim\text{binom}(n,\lambda t/n). Lader vi nn blive større og større, kan man matematisk vise, at
P(X=x)=(nx)(λt/n)x(1−λt/n)n−x→(tλ)xx!e−tλ. P(X=x)=\binom{n}{x}(\lambda t/n)^x(1-\lambda t/n)^{n-x} \rightarrow \frac{(t\lambda)^x}{x!}e^{-t\lambda}.
Her står, at i grænsen hvor nn bliver meget stor, vil XX være poissonfordelt med parameter tλ,t\lambda, X∼poisson(tλ).X\sim\text{poisson}(t\lambda).
Man udtrykker ovenstående argument på den måde, at en binomialfordeling binom(n,p)\text{binom}(n,p) med nn stor og pp lille ligner en poissonfordeling med parameter np.np. Nedenstående kommandovindue laver en figur til illustration af dette. I figuren er poissonsandsynlighederne angivet ved søjler og binomialsandsynlighederne er angivet ved røde plustegn. Kør koden, og prøv dernæst at sætte nn til 100. Kan du se hvad der sker ? Prøv også en anden værdi af Lambda.

Test dig selv: valg af statistikmodel

Denne lille testopgave går ud på at træne jer i at kunne vurdere, om data skal beskrives med en binomialmodel eller med en poissonmodel.
Quiz
I løbet af en måned er der indsamlet 20 vandprøver fra forskellige vandboringer, og 5 af disse har pesticidrester over det tilladte niveau. Klik på de rigtige udsagn nedenfor, hvor Pest står for den bagvedliggende stokastiske variabel.
De 5 prøver kommer tilfældigt over tid i løbet af en måned, hvorfor det er en observation fra en poissonfordeling, Pest∼poisson(λt),\text{Pest}\sim\text{poisson}(\lambda t), hvor tt er en måned.
Hver af de 20 prøver kan have pesticidrester enten over eller under grænsen, hvorfor 5 er en observation fra en binomialfordeling, Pest∼binom(20,p).\text{Pest}\sim\text{binom}(20,p).
Prøv også at forklare, hvorfor det forkerte valg er forkert.

Svar: Forkerte statistikmodel

Der er ikke tale om, at der kommer prøver med højt indhold af pesticidrester tilfældigt over tid. Tidsaspektet i eksperimentet har ikke noget at gøre med indeholdet af pesticidrester, men siger blot noget om, hvornår man vælger at tage ud og indsamle vandprøven.

2.5.1 Estimation i poissonmodellen

Betragt uafhængige poissonfordelte variable X1,…,Xn,X_1,\ldots,X_n, Xi∼poisson(tλ),X_i\sim\text{poisson}(t\lambda), hvor tt er kendt og λ>0\lambda>0 er en ukendt parameter, som vi ønsker at estimere. Bemærk her, at parameteren i poissonfordelingen er tλt\lambda, hvor vi tænker på tt som et tidsrum, og λ\lambda fortolkes som den forventede rate per tid. På grund af uafhængigheden bliver likelihoodfunktionen et produkt af poissonsandsynligheder
L(λ)=∏i=1n(tλ)xixi!e−tλ=(1∏ixi!)(tλ)x1+⋯+xne−ntλ. L(\lambda)=\prod_{i=1}^n \frac{(t\lambda)^{x_i}}{x_i!}e^{-t\lambda}= \Big(\frac{1}{\prod_ix_i!} \Big) (t\lambda)^{x_1+\cdots+x_n}e^{-nt\lambda}.
For at finde den værdi af λ\lambda som giver maksimum af likelihoodfunktionen, tager vi logaritmen, differentierer med hensyn til λ\lambda og sætter den afledede lig med nul. Dette giver
x1+⋯+xnλ−nt=0ellerλ^=xˉt,xˉ=x1+⋯+xnn. \frac{x_1+\cdots+x_n}{\lambda}-nt=0\quad\text{eller}\quad \hat\lambda=\frac{\bar x}{t},\quad \bar x=\frac{x_1+\cdots+x_n}{n}.
Ligesom i binomialmodellen giver dette skøn god mening: λ\lambda er den forventede rate per tid, og λ^\hat\lambda er den observerede rate per tid.
Eksempel 2.5.2. (Rutherford og Geiger)
I en berømt artikel fra 1910 giver Rutherford og Geiger resultaterne fra 2608 målinger af antal henfald fra en poloniumsmasse i 72 sekunder. Vi kan beskrive situationen med Xi∼poisson(72⋅λ),X_i\sim\text{poisson}(72\cdot\lambda), i=1,…,2608.i=1,\ldots,2608. I denne situation bliver skønnet over raten per sekund λ^=10097/(72⋅2608)=0.05377,\hat\lambda=10097/(72\cdot 2608)=0.05377, idet summen af de 2608 målinger er 10097. Et histogram med de 2608 målinger er vist i den følgende figur.

2.5.2 Konfidensinterval i poissonmodellen

I poissonmodellen X∼poisson(tλ)X\sim\text{poisson}(t\lambda) er det ikke muligt at lave et konfidensinterval for λ,\lambda, som opfylder betingelsen i Definition 2.2.1 eksakt. I stedet benytter man følgende approksimative konfidensinterval.
Resultat 2.5.3. (Konfidensinterval i poissonmodellen)
For modellen X∼poisson(tλ)X\sim\text{poisson}(t\lambda) er et approksimativt 95%-konfidensinterval for λ\lambda baseret på observationen xx givet ved
[1t(x+u22−ux+u24), 1t(x+u2/2+ux+u2/4)],u=1.96. \Big[\frac{1}{t}\Big(x+\frac{u^2}{2}-u\sqrt{x+\frac{u^2}{4}}\Big),\, \frac{1}{t}\Big(x+u^2/2+u\sqrt{x+u^2/4}\Big)\Big],\quad u=1.96.
Hvis værdien af uu ændres til 1.645, får man i stedet et approksimativt 90%-konfidensinterval. Ligesom i binomialmodellen bygger konfidensintervallet her på den centrale grænseværdisætning.
Eksempel 2.5.4. (Rutherford og Geiger)
I Eksempel 2.5.2 omkring henfald af polonium er der 2608 uafhængige målinger: Xi∼poisson(72⋅λ),X_i\sim\text{poisson}(72\cdot\lambda), i=1,…,2608.i=1,\ldots,2608. Som nævnt ovenfor er summen af XiX_i-erne også poissonfordelt. Vi har derfor, at X=∑iXi∼poisson(72⋅2608⋅λ).X=\sum_i X_i\sim\text{poisson}(72\cdot 2608\cdot\lambda). Den målte værdi af XX er x=10097.x=10097. I kodevinduet nedenfor udregnes det approksimative 95%-konfidensinterval.

2.5.5 Beregning i python af konfidensinterval for rate

Den følgende kode kan bruges generelt til beregning af konfidensinterval for en rate, hvis man selv indskriver det observerede antal (xx) og tidsparameteren (tt).

ForegåendeNæste