I 2009 publicerede Pieter Vermeesch en lille note med den provokerende titel Lies, damned lies, and statistics (in geology). I noten betragter forfatteren 118415 jordskælv af styrke 4 eller over på Richterskalaen i perioden 1/1-1999 til 1/1-2009 (data fra earthquake.usgs.cov) og fordeler disse på ugedag. Billedet her viser, hvor jordskælv optræder. Data kan ses i den følgende tabel og findes i filen JordskaelvDag.csv fra opgave 1.4. $$\begin{array}{lccccccc} \hline \text{Data} & \text{Mandag} & \text{Tirsdag} & \text{Onsdag} & \text{Torsdag} & \text{Fredag} & \text{Lørdag} & \text{Søndag} \\ \hline \text{Observeret} & 16853&16553&16490&17399&16348&17019&17753 \\ \hline \end{array}$$

1. Opstil multinomialmodellen for disse data, hvor sandsynlighederne for at falde i de syv kasser er vilkårlige.
2. Opskriv, inden for den opstillede multinomialmodel, hypotesen om ligelig fordeling på de syv ugedage. Udregn de forventede antal, og lav $G$-testet for hypotesen. Hvad bliver konklusionen af testet ?

Den overraskende konklusion er, at jordskælvene ikke fordeler sig ligeligt på de syv ugedage. Vermeesch ser dette resultat som udtryk for en svaghed i de statistiske metoder. Har Vermeesch ret, eller bruger han metoderne forkert ? Når vi laver goodness of fit test for en ligelig fordeling på de syv ugedage, sammenligner vi med uafhængige kast af en syvsidet terning. Som påpeget af Sornette og Pisarenko i artiklen On the correct use of statistical tests: Reply to "Lies, damned lies and statistics (in geology)", er det helt store problem mangel på uafhængighed på grund af efterskælv og klynger af små skælv. I datasættet er cirka $2/3$ af skælvene efterskælv! Et andet stort problem er, at for de små jordskælv kan data indeholde menneskeskabte hændelser, og der kan være forskel i baggrundsstøj på de syv ugedage. Sornette og Pisarenko foreslår at betragte jordskælv med en styrke over 5 på Richterskalaen for at komme ud over det sidstnævnte problem, og foretager også en rensning af disse for at fjerne efterskælv (dette kan ikke gøres fuldstændigt sikkert). Data kan ses i den følgende tabel. $$\begin{array}{lccccccc} \hline \text{Data} & \text{Mandag} & \text{Tirsdag} & \text{Onsdag} & \text{Torsdag} & \text{Fredag} & \text{Lørdag} & \text{Søndag} \\ \hline \text{Over 5} & 2374&2511&2291&2497&2153&2282&2360 \\ \text{Renset} & 780&847&793&831&785&821&779 \\ \hline \end{array}$$

3. Lav et $G$-test for en ligelig fordeling, både for alle jordskælv med en styrke over 5 og for delmængden, hvor efterskælv er fjernet.

Ifølge Wikipedia blev udtrykket lies, damned lies and statistics gjort kendt af Mark Twain, som skrev "There are three kinds of lies: lies, damned lies, and statistics". Mark Twain mente, at udtrykket stammede fra den britiske premierminister Benjamin Disraeli, men dette kan tilsyneladende ikke bekræftes. Vermeesch var utilfreds med, at han fik en lav $p$-værdi for en hypotese, der synes oplagt sand. Dette er dog en mulighed, vi må være indstillet på, jævnfør omtalen af fejl af type 1 i afsnit 1.4.

Luk igen

Data i denne opgave vedrører positionen af kratere på Mars. Hvis kratere er spredt tilfældigt ud over Mars, vil antallet af kratere i et område følge en poissonfordeling. Data er hentet i Mars Crater Database Search, og vi betragter kun kratere med en diameter over 4 km. Området mellem 0 og 15 graders bredde og 0 og 45 graders længde er delt op i $15\times 15$ områder af lige stort areal. Dette er vist i den følgende figur, og antallet af kratere i de 225 områder findes i filen KraterAntal.txt.

1. Indlæs de 225 tællinger ved ordren nKrater=scan("KraterAntal.txt").
   Lav et antalshistogram af data med intervalendepunkter endePkt=c(1:19)-0.5 (det første interval er fra 0.5 til 1.5, svarende til at den mindste værdi i nKrater er 1, og det sidste interval er fra 17.5 til 18.5, svarende til at den største værdi i nKrater er 18). Indsæt titler på akserne i figuren ved at benytte xlab og ylab i kaldet til hist. Placer antallet af observationer i hvert interval i en vektor antal. Vælg et af intervallerne ud, og eftervis antallet i antal ved en direkte optælling blandt de 225 dataværdier.
2. Opskriv multinomialmodellen for den stokastiske antalsvektor Antal, hvor sandsynligheden for at falde i de forskellige kasser er vilkårlig.

I denne opgave skal I lave et goodness of fit test for at antallet af kratere i et område er poissonfordelt. Lad $\lambda$ være rateparameteren i poissonfordelingen (enhed: forventede antal per område). Da den første "kasse" i multinomialmodellen indeholder tælletal mindre end eller lig med 1, vil hypotesen om en poissonfordeling betyde at sandsynligheden for at falde i den første kasse er $\sum_{j=0}^1(\lambda^j/j!)e^{-\lambda}.$ Den anden kasse indeholder alle tælletal med værdien 2, og sandsynligheden for at falde i den anden kasse er $(\lambda^2/2!)e^{-\lambda},$ og så videre op til kasse 17. Sandsynligheden for at falde i den sidste kasse (kasse nummer 18) er 1 minus summen af sandsynlighederne for de første 17 kasser. Som skøn over $\lambda$ bruges gennemsnittet af de 225 observationer, se underafsnit 2.5.1.

3. Opskriv, inden for din multinomialmodel, hypotesen om, at antallet af kratere i et område er poissonfordelt.
   Beregn de forventede antal under hypotesen. Hertil kan du benytte koden nedenfor. I R beregnes punktsandsynligheder i poissonfordelingen med dpois(x,lambda), og sandsynligheden for en værdi mindre end eller lig med $x$ beregnes med ppois(x,lambda). Forklar, at koden giver de forventede værdier.
   Indtegn de forventede antal i histogrammet fra spørgsmål (a) som en rød kurve med kommandoen lines(c(1:18),forvent,col=2), hvor forvent er vektoren med de forventede antal.
4. Lav $G$-testet for hypotesen, at antal kratere i et område er poissonfordelt. Slå kasser sammen, hvis de forventede ikke er større end 5 (slå kasser sammen fra hver sin ende, indtil det forventede antal er større end 5).
   Hvad bliver konklusionen af dit goodness of fit test ?

Ovenstående undersøgelse var baseret på kratere med en diameter over 4 km. Hvis man tager kratere med, som har en diameter mellem 1 og 4 km, viser det sig, at antallet i et rektangel ikke længere kan beskrives med en poissonfordeling. Kan du se nogle forklaringer på dette ?

Luk igen

Data i denne opgave er mængden af krom i 61 prøver af vulkanske lavabjergarter fra Ontario i Canada. Mængden af krom måles ved spektrokemiske metoder. Data er fra artiklen Minor element content of Ontario diabase. I artiklen er 9 værdier blot angivet som større end 0.077. I de data I skal bruge, er disse 9 værdier erstattet af simulerede værdier. Data ligger i filen Ontario.txt. Opgaven går ud på at lave et goodness of fit test for, at kromindholdet kan beskrives med en normalfordeling.

1. Indlæs de 61 mængder af krom (enheden er koncentrationsprocent) fra filen Ontario.txt, placer disse i vektoren krom, og dan vektoren logKrom=log(krom) med logaritmen til mængden. Lav et tæthedshistogram af data med intervalinddelingen endePkt=-8+c(0:13)*0.5. Placer antallet af observationer i hvert af de 13 intervaller i en vektor antal.
   Hvis logaritmen til mængden af krom skal beskrives med en normalfordeling, er det bedste valg af middelværdi $\hat\mu=-4.8961,$ og det bedste valg af spredning er $\hat\sigma=1.6555.$ Indtegn normalfordelingstætheden i histogrammet med kommandoen

   curve(dnorm(x,-4.8961,1.6555),from=-8,to=-1,add=TRUE)
2. Opskriv multinomialmodellen for den stokastiske vektor Antal, hvor sandsynlighederne for at falde i de forskellige intervaller er vilkårlige.
   Opskriv dernæst hypotesen, at sandsynlighederne for at falde i de 13 intervaller er givet ved sandsynlighederne for intervallerne i en normalfordeling med middelværdi $\mu$ og spredning $\sigma.$ Husk at i denne sammenhæng skal det første interval opfattes som intervallet fra minus uendelig til -7.5, og det sidste interval skal opfattes som intervallet fra -2.0 til uendelig.
   Beregn de forventede antal under hypotesen. Hertil kan du benytte koden nedenfor. I R beregnes sandsynligheden for en værdi mindre end eller lig med $x$ i en normalfordeling med kommandoen pnorm(x,$\mu$,$\sigma$). Forklar, at koden giver de forventede værdier.
3. Lav $G$-testet for hypotesen, at kromindholdet er normalfordelt. Kan disse data beskrives med en normalfordeling ?

Luk igen

Data i denne opgave stammer fra rapporten Stratigraphy, lithology, and sedimentary features of quaternary alluvial deposits of the South Platte River and some of its tributaries east of the Front Range, Colorado. Der er indsamlet 100 sten fire steder i South Platte River flodløbet. (Overvej, hvordan man indsamler 100 tilfældige sten!) Stenene er klassificeret efter litologi, og i denne opgave vil vi kun betragte grupperne Gneiss, Pegmatite, Vein quartz og Quartzite (hvorfor totaltallene i den følgende tabel ikke er 100). $$\begin{array}{lccccc } \hline \text{Område} & \text{Gneiss} & \text{Pegmatite} & \text{Vein quartz} & \text{Quartzite} & \text{Total} \\ \hline \text{Rocky Flats 47} & 21 & 9 & 3 & 64 & 97 \\ \text{Rocky Flats 50} & 16 & 5 & 8 & 69 & 98 \\ \hline \text{Big Dry Creek 62} & 24 & 30 & 8 & 7 & 69 \\ \text{Big Dry Creek 64} & 20 & 33 & 9 & 7 & 69 \\ \hline \end{array}$$ En eventuel forskel i fordelingen på litologi kan indgå i en forståelse af de geologiske processer, der har formet floden.

1. Opstil den statistiske model, hvor tælletallene for hvert område i Rocky Flats følger sin egen multinomialfordeling. Angiv inden for den opstillede model hypotesen, at der er samme sandsynlighedsvektor for kategorierne (Gneiss, Pegmatite, Vein quartz, Quartzite) for de to områder.
2. Undersøg, om data er i overensstemmelse med hypotesen om samme sandsynlighedsvektor for kategorierne (Gneiss, Pegmatite, Vein quartz, Quartzite) for de to Rocky Flats områder (benyt eventuelt R-koden fra eksempel 3.7.2).

En tilsvarende undersøgelse for de to Big Dry Crek områder viser, at det også her kan antages at fordelingen på de fire kategorier (Gneiss, Pegmatite, Vein quartz, Quartzite) er den samme. I skal nu sammenligne Rocky Flats områderne med Big Dry Creek områderne. Baseret på resultaterne ovenfor vil det være naturligt at se på sum-tallene fra de to dele af flodsystemet. Disse tal er vist i den følgende tabel. $$\begin{array}{lccccc } \hline \text{Område} & \text{Gneiss} & \text{Pegmatite} & \text{Vein quartz} & \text{Quartzite} & \text{Total} \\ \hline \text{Rocky Flats 47+50} & 37 & 14 & 11 & 133 & 195 \\ \text{Big Dry Creek 62+64} & 44 & 63 & 17 & 14 & 138 \\ \hline \end{array}$$

3. Opstil en statstisk model for disse sum-data. Opskriv hypotesen, at der er samme fordeling på litologi for de to dele af flodsystemet, og lav et test for denne hypotese. Hvad bliver resultatet af dette test ?

Luk igen

Lad os vende tilbage til data i opgave 2.4 med tælletal for antallet af jordskælv i tre styrkeintervaller. I opgaven her betragtes kun de to første styrkeintervaller, 6.0-6.3 og 6.3-6.6. Idet vi vil bruge data til at vurdere holdbarheden af Gutenberg-Richter loven, skrives raterne i de to intervaller som $86\lambda_1$ og $0.5\cdot 86\lambda_2,$ eftersom en beregning viser, at under Gutenberg-Richter loven (med "b-value" lig med 1) er raten i det andet styrkeinterval halvt så stor som raten i det første styrkeinterval. Enheden på $\lambda$ er antal forventede jordskælv per år. Gutenberg-Richter loven svarer således til hypotesen $\lambda_1=\lambda_2.$ Lad os formulere situationen generelt gennem modellen $$Y_1\sim\text{pois}(t_1\lambda_1),\quad Y_2\sim\text{pois}(t_2\lambda_2),\quad t_1=86,\enspace t_2=0.5\cdot 86,$$ hvor $Y_1$ og $Y_2$ er de stokastiske tælletal svarende til de to styrkeintervaller. For at teste hypotesen $\lambda_1=\lambda_2$ kan man benytte følgende teoretiske resultat. Hvis vi forestiller os, at $Y_1+Y_2$ er fast (vi "betinger" med summen), så vil $Y_1$ være binomialfordelt: $$\big(Y_1\enspace\text{givet at}\enspace Y_1+Y_2=n\big)\sim \text{binom}\big(n,p \big), \quad p=\frac{t_1\lambda_1}{t_1\lambda_1+t_2\lambda_2}. \tag{3.9.1}$$ Hvis $\lambda_1=\lambda_2$ bliver $p=t_1\big/(t_1+t_2).$ Et test for hypotesen $\lambda_1=\lambda_2$ kan derfor laves som et test i binomialfordelingen for hypotesen, at sandsynlighedsparameteren $p$ har værdien $p=t_1\big/(t_1+t_2).$ I vores tilfælde bliver dette hypotesen, at $p=2/3.$

1. Find $p$-værdien for et test af hypotesen $\lambda_1=\lambda_2$ for data i opgave 2.4 ved at teste $p=2/3$ i modellen (3.9.1). Hvad bliver konklusionen af dette test ?

Bemærkning:$\enspace$Hvis data strider mod hypotesen $\lambda_1=\lambda_2,$ vil vi være interesseret i at indføre en parameter $\theta,$ således at $\lambda_2=\theta\lambda_1.$ Parameteren $\theta$ angiver, hvor mange gange større $\lambda_2$ er i forhold til $\lambda_1.$ Hvis vi lader $p=t_1\big/(t_1+\theta t_2),$ kan vi løse for $\theta,$ og får $\theta=(1-p)t_1\big/(pt_2).$ Konfidensinterval for sandsynlighedsparameteren $p$ i binomialmodellen (3.9.1) kan derfor oversættes til et konfidensinterval for forholdet $\theta.$ Hvis denne metode benyttes på data i første og tredje styrkeinterval i opgave 2.4, får vi intervallet $[0.65,\,0.82]$ for 95%-konfidensintervallet af sandsynlighedsparameteren $p.$ Oversat til forholdet $\theta$ mellem rateparametrene i de to styrkeintervaller giver dette intervallet $[0.22,\,0.55].$ Dette interval indeholder værdien $\frac{1}{4},$ som svarer til Gutenberg-Richter loven.

Luk igen

DMI vedligeholder en side med alle storme i Danmark fra 1891 og fremefter. Stormene kalssificeres i fire styrkekategorier ud fra vindstyrken. I nedenstående tabel har jeg optalt antallet af storme i de forskellige kategorier for fire 30-års perioder. $$\begin{array}{lccc } \hline \text{Periode} & \text{Stormstyrke 1} & \text{Stormstyrke 2} & \text{Stormstyrke 3 og 4} \\ \hline \text{1891-1920} & 39 & 16 & 4 \\ \text{1921-1950} & 21 & 8 & 8 \\ \text{1951-1980} & 14 & 12 & 5 \\ \text{1981-2010} & 18 & 12 & 10 \\ \hline \end{array}$$

1. Opstil den statistiske model, hvor antallet af storme for hver periode følger sin egen multinomialfordeling på de tre kategorier 1, 2 og 3+4. Angiv inden for den opstillede model hypotesen, at der er samme sandsynlighedsvektor for de tre styrkekategorier 1, 2 og 3+4 for de fire tidsperioder.
2. Undersøg, om data er i overensstemmelse med hypotesen om samme sandsynlighedsvektor for kategorierne 1, 2 og 3+4 for de fire tidsperioder.

For storme i Danmark finder I således ikke en ændring over tid i sandsynlighedsvektoren for de tre styrkekategorier. Den næste tabel viser klassifikationen af hurricanes fra verdenshavene på styrkekategorierne 1-3 og 4-5. Der er data fra to tidsperioder: 1975-1989 og 1990-2004. Data er fra artiklen Changes in Tropical Cyclone Number, Duration, and Intensity in a Warming Environment. $$\begin{array}{lcc } \hline \text{Periode} & \text{Hurricanes 1-3} & \text{Hurricanes 4-5} \\ \hline \text{1975-1989} & 617 & 171 \\ \text{1990-2004} & 508 & 269 \\ \hline \end{array}$$

3. Undersøg om sandsynlighedsvektoren for de to styrkekategorier er den samme for de to tidsperioder.

For hurricanes ser det således ud til, at der er sket en ændring i sandsynlighedsvektoren, således at der er kommet flere kategori 4-5 hurricanes relativt til kategori 1-3 hurricanes. Vi er derfor interesserede i, om denne øgede mængde af kategori 4-5 hurricanes kun ses i visse verdenshave eller ses overalt. Du skal derfor se på hvordan kategori 4-5 hurricanes, for en given tidsperiode, fordeles ud på 5 verdenshave, og sammenligne de to tidsperioder. Den næste tabel viser fordelingen af de største hurricanes (kategori 4-5) på fem verdenshave for de to tidsperioder. Data er i den næste tabel, hvor søjleoverskrifterne er NA: North Atlantic, EPO: East Pacific Ocean, WPO: West Pacific Ocean, SwP: Southwestern Pacific og In: Indian. $$\begin{array}{lccccc } \hline \text{Verdenshav:} & \text{NA} & \text{EPO} & \text{WPO} & \text{SwP} & \text{In} \\ \hline \text{1975-1989} & 16 & 36 & 85 & 10 & 24 \\ \text{1990-2004} & 25 & 49 & 116 & 22 & 57 \\ \hline \end{array}$$

4. Opstil model for data, og undersøg, om sandsynlighedsvektoren for fordelingen af kategori 4-5 hurricanes på de fem verdenshave er den samme i de to tidsperioder.

Luk igen