Hvorfor simulering?

Fra og med 2024 inneholder hvert eneste TMA4245-sett minst én oppgave i stokastisk simulering, og de fleste krever at du både beskriver framgangsmåten med ord og implementerer den i Python. Simulering brukes når fordelingen til en størrelse er vanskelig eller umulig å finne analytisk: eksakte p-verdier, eksakte dekningsgrader for tilnærmede konfidensintervall, og kvantiler i sammensatte fordelinger.

Inversjonsmetoden

Alt starter med uniforme tilfeldige tall. Er $F$ en kontinuerlig, strengt voksende kumulativ fordelingsfunksjon og $U\sim\text{Uniform}(0,1)$, så har $X=F^{-1}(U)$ nettopp fordelingen $F$:

$$P(F^{-1}(U)\le x)=P(U\le F(x))=F(x).$$

Oppskrift:

Sett $u=F(x)$ og løs for $x$. For $F(x)=1-e^{-\lambda x}$ gir dette $x=-\ln(1-u)/\lambda$. Siden $1-U$ også er uniform, kan man like gjerne bruke $-\ln(U)/\lambda$.

Eksempel: har $X$ tetthet $f(x)=\frac{2x}{\theta}e^{-x^2/\theta}$, er $F(x)=1-e^{-x^2/\theta}$ og $X=\sqrt{-\theta\ln(1-U)}$. For diskrete fordelinger brukes den generaliserte inversen: velg den minste $x_k$ med $F(x_k)\ge U$, i kode x = np.sum(F < u) når F er vektoren av kumulative sannsynligheter.

Monte Carlo-estimering

Skal vi estimere $p=P(A)$, gjentar vi eksperimentet $m$ ganger og teller. Med indikatorene $I_1,\dots,I_m$ er $\hat{p}=\frac{1}{m}\sum I_i=Z/m$ med $Z\sim\text{bin}(m,p)$, altså

$$E[\hat{p}]=p,\qquad \text{Var}(\hat{p})=\frac{p(1-p)}{m},\qquad \text{SE}(\hat{p})=\sqrt{\frac{\hat{p}(1-\hat{p})}{m}}.$$

Estimatoren er EKSAKT forventningsrett for enhver $m$ — dette er argumentet sensor ser etter. Tilsvarende estimeres $E[g(X)]$ med $\frac{1}{m}\sum g(X_i)$. Merk at $\left(\frac{1}{m}\sum X_i\right)^2$ IKKE estimerer $E[X^2]$, fordi $g$ ikke er lineær.

Monte Carlo-feilen avtar som $1/\sqrt{m}$: skal standardfeilen halveres, må $m$ firedobles. Med $m=10^4$ og $\hat{p}=0{,}04$ er $\text{SE}=\sqrt{0{,}04\cdot0{,}96/10^4}=0{,}0020$.

Eksakt p-verdi

Regn først ut observert testobservator $t_0$. Simuler så $m$ datasett UNDER $H_0$, regn ut $T_j$ for hvert, og sett $I_j=1$ når $T_j$ er minst så ekstrem som $t_0$. Da er

$$\hat{p}_0=\frac{1}{m}\sum_{j=1}^m I_j,\qquad E[\hat{p}_0]=P(T\ \text{minst så ekstrem}\mid H_0).$$

Fordi vi trekker fra den KORREKTE nullfordelingen, er estimatet forventningsrett for den eksakte p-verdien — i motsetning til en normaltilnærming, som har en systematisk feil som ikke forsvinner uansett hvor stor $m$ er.

Eksakt dekningsgrad

Et intervall utledet fra sentralgrenseteoremet har nominelt nivå $1-\alpha$, men bare asymptotisk. Dekningsgraden er $P(L(X)\le\theta\le U(X))$ der sannsynligheten er over DATAENE, med $\theta$ fast. Algoritmen er derfor: velg en sann verdi $\theta_0$, simuler et nytt datasett av størrelse $n$, regn ut intervallet, og tell om $\theta_0$ ligger inni. Andelen treff er forventningsrett for den eksakte dekningen.

Vanlig eksamensvri:

Simuleringen gir f.eks. $0{,}906$ mot nominelt $0{,}95$ med $m=10^4$. Standardfeilen er $\sqrt{0{,}906\cdot0{,}094/10^4}=0{,}0029$, så avviket er 15 standardfeil: underdekningen er reell, og normaltilnærmingen er for grov for denne $n$.

Kvantiler ved simulering

Trenger du kvantiler i en pivotal som ikke står i tabellsamlingen, simulerer du mange realisasjoner og leser av de empiriske kvantilene med np.quantile(y, 1-alpha). Dette er den vanlige avslutningen på prediksjonsintervall-oppgaver der pivotalen er et forhold mellom to kjikvadratstørrelser.

Python-verktøykassen

np.random.uniform(0,1,n), np.random.normal(mu, sigma, n) (merk: STANDARDAVVIK, ikke varians), np.random.exponential(beta, n) (skala = forventning), np.random.gamma(alfa, beta, n), np.mean, np.quantile. Boolske vektorer castes til 0/1, så np.mean(x > c) er direkte et estimat av $P(X>c)$, mens np.sum gir antallet.

Studieguide: Stokastisk simulering og Python