Afsnit 4.2: Normal-qqplot

Eksempel 4.2.1. (Jordens massetæthed)

Henry Cavendish lavede i 1797 et eksperiment til bestemmelse af Jordens massetæthed. Eksperimentet går ud på at bestemme den universelle gravitationskonstant $G$ i Newtons lov, som siger, at tiltrækningskraften mellem to legemer er $F=GmM/r^2$ . Her er $m$ og $M$ masserne af de to legemer og $r$ er afstanden mellem dem. Når vi ser på kraften hvormed jorden trækker i et legeme med masse $m$ kan vi skrive denne både som $mg$ , hvor $g$ er tyngdeaccelerationen og som $GmM_J/R_J^2,$ hvor $M_J$ er Jordens masse og $R_J$ er Jordens radius. Hvis vi lader $\rho$ være Jordens massetæthed, har vi $\rho=M_J/(4\pi R_J^3/3,$ og kan derfor udlede, at $\rho=3g/(4\pi R_JG)$ , således at bestemmelse af $G$ giver en bestemmelse af $\rho.$ Cavendish's eksperiment til bestemmelse af $G$ er illustreret i følgende figur, hvor tiltrækningskraften mellem to kugler bestemmes ud stivheden i en wire der holder den ene kugle. Cavendish's eksperiment er lavet med to forskellige wires og I afsnit 4.5 betragter jeg de 23 målinger hørende til den ene wire. Jeg vil undersøge om Cavendish er nået frem til den samme massetæthed for Jorden som i moderne eksperimenter, det vil sige værdien 5.517 relativt til tætheden for vand (enhed brugt af Cavendish). Jeg er også interesseret i variationen i data som udtryk for kvaliteten af Cavendish's eksperiment.

Jeg vil gerne bruge normalfordelingen til at beskrive mine data, men er dette rimeligt ? Jeg har ikke nok data til at lave et goodness of fit test som beskrevet i afsnit 3.4, og som I prøvede at lave i Opgave 3 i Øvelse 2. I stedet vil jeg her beskrive en grafisk undersøgelse, der kan give en indikation af, om det er rimeligt at bruge en normalfordeling. I den grafiske metode laves en figur, hvor punkterne bør "sno sig" omkring en ret linje, i fald data stammer fra en normalfordeling. Med kun ganske få datapunkter, som i mit eksempel ovenfor, kan det være svært at afgøre, om data afviger fra at "sno sig" omkring en ret linje. Den grafiske undersøgelse er således af større værdi, hvor man har flere datasæt og kan se, om de alle viser den samme form for afvigelse fra "sno sig" egenskaben.

Resultat 4.2.2. (Brugen af normal-qqplot)

I python konstrueres et (normal-) qqplot med kommandoen qqplot. Punkterne der tegnes vil altid være stigende ( $y_2\geq y_1$ for to på hinanden følgende punkter), og som beskrevet til sidst i dette afsnit, forventer vi, at hvis data stammer fra en normalfordeling, så vil punkterne sno sig omkring en ret linje.

Det følgende kodevindue laver den grafiske undersøgelse for data i eksemplet ovenfor.

4.2.3 QQplot i python

I python hedder den relevante funktion qqplot. For at figuren bliver som beskrevet nedenfor, skal man altid tilføje a=0.5 i kaldet. Man kan få en linje indsat, der går gennem første og tredje kvartil af data ved at tilføje line='q' i kaldet til qqplot (python bruger en algoritme til beregning af kvartiler, der ikke er optimal, når antallet af observationer er lille). Figuren, der dannes i kodevinduet her, understøtter, at Cavendish's målinger kan beskrives med en normalfordeling. Prøv at finde værdien 5.517 på andenaksen i figuren, og vurder om denne ligger centralt blandt målingerne.

Jeg beskriver nu den grafiske undersøgelse, lavet i kodevinduet ovenfor, som går under navnet normal-qqplot. Her står "q" for quantile, som på dansk er fraktil, og på dansk taler man om en fraktilsammenligning. For nemhed i notationen vil jeg fremover blot omtale metoden som et qqplot. For at beskrive metoden lader jeg $u_p$ være $p$ -fraktilen i en standard normalfordeling, $N(0,1)$ -fordelingen, det vil sige, at $N_{\text{cdf}}(u_p,0,1)=p.$ I python beregnes $p$ -fraktilen som st.norm.ppf(p).

Vi betragter $n$ datapunkter $x_1,x_2,\ldots,x_n$ og ordner disse efter størrelse, $x_{[1]}$ betegner den mindste, $x_{[2]}$ den næstmindste, og så videre op til $x_{[n]}$ som er den største: $x_{[1]}\leq x_{[2]}\leq\cdots\leq x_{[n]}.$ Et qqplot består i at tegne punkterne

$(u_{(i-0.5)/n},x_{[i]}),\quad i=1,2,\ldots,n.$ I python produceres denne figur, når der i kaldet til qqplot tilføjes $a=0.5.$

QQplots af simulerede data

I nedenstående kodevinduer vises nogle eksempler på qqplots. Først kommer et kodevindue, hvor data er normalfordelt. Der laves en figur med fire qqplots, alle med det samme antal observationer. Prøv at køre koden et par gange. Prøv dernæst at ændre $n$ fra 10 til 40, og dernæst til 100. Kommandoen norm.rvs(0,1,n) simulerer $n$ observationer fra en standard normalfordeling.

Nu følger et kodevindue, der danner en figur med fire qqplots, og hvor data ikke er normalfordelte for de sidste to qqplots. Prøv at køre koden et par gange. Prøv dernæst at ændre $n$ fra 10 til 40, og dernæst til 100.

I det øvre højre delplot betragtes data fra en normalfordeling. Hvad er middelværdi og spredning i denne normalfordeling ?
I det nedre venstre delplot betragtes data fra en stokastisk variabel, der kun kan antage positive værdier. Hvad er fordelingen af logaritmen til den stokastiske variabel ?

Svar: Fordelinger

Der simuleres først data data fra en $N(0,1)$ fordeling, hvorefter disse ganges med 2 og der lægges 3 til. Dette betyder, at de nye data stammer fra en $N(3,2^2)$ -fordeling, middelværdi er 3 og spredning er 2.
Data kommer fra en stokastisk variabel $Y=e^X,$ hvor $X$ er $N(0,1)$ -fordelt. Hvis vi tager logaritmen, får vi $\log(Y)=\log(e^X)=X,$ som er normalfordelt. Man siger, at $Y$ er log-normalfordelt.

Hvorfor giver et qqplot en figur, hvor normalfordelte data snor sig omkring en ret linje ? Her er kort den tekniske ide bag et qqplot. Hvis data følger en $N(\mu,\sigma^2)$ -fordeling, så bør der gælde, at for ethvert $p$ mellem 0 og 1 vil $p\text{-}$ fraktilen beregnet ud fra data $x_1,\ldots,x_n$ ligne $p$ -fraktilen i en $N(\mu,\sigma^2)\text{-}$ fordeling. Hvis $X\sim N(\mu,\sigma^2),$ kan vi skrive $X=\mu+\sigma U,$ hvor $U\sim N(0,1),$ hvorfor $p$ -fraktilen for $X$ kan skrives som $\mu+\sigma u_p.$ Hvad mener jeg med fraktiler beregnet ud fra data $x_1,\ldots,x_n$ ? Per definition af de ordnede værdier vides, at i punktet $x_{[i]}$ er andelen af data mindre end eller lig med denne værdi givet ved $i/n,$ men hvis vi betragter en værdi $x$ lidt mindre end $x_{[i]}$ (men større end $x_{[i-1]}$ ), er andelen af data mindre end eller lig med denne værdi i stedet $(i-1)/n.$ Vi vælger derfor at sige, at $x_{[i]}$ er et skøn over $((i-\frac{1}{2})/n)$ -fraktilen. Vores argument er derfor, at hvis data er $N(\mu,\sigma^2)\text{-}$ fordelt, så bør

$x_{[i]} \approx \mu+\sigma u_{(i-\frac{1}{2})/n},\quad i=1,\ldots,n,$ hvor " $\approx$ " skal læses som "cirka lig med". I qqplottet bør punkterne derfor sno sig om en linje med hældning $\sigma.$

Shapiro-Wilks test som supplement

Et QQplot giver et godt visuelt indblik i data. Hvis antallet af observationer $n$ er lille, kan det dog være svært at vurdere, om data kan siges at følge en normalfordeling. Mange statistikpakker indeholder muligheden for at lave et mere formelt test. En af de populære metoder er Shapiro-Wilks test. Jeg vil ikke gå ind på en nærmere beskrivelse her, men blot nævne at testet kan udføres med kommandoen st.shapiro(x), hvor $x$ er en vector med data. Dette vises i kodevinduet nedenfor med målingerne af jordens massetæthed foretaget af Cavendish.

Foregående Næste