Afsnit 1.7: Øvelse 1: Geologi

Denne første uges øvelse skal især gøre jer fortrolige med programpakken R. Derudover skal I vurdere holdbarheden af en hypotese ved hjælp af simulationer (at simulere betyder her at generere tilfældige tal fra en given fordeling). Vi starter dog først ud med en opgave, der skal vise jer nogle eksempler på den type data, I kan støde på i geologi. Alle opgaverne skal være forberedt hjemmefra og gennemgås ved tavlen til øvelserne.

Opgave 1.1: Dataeksempler

I denne opgave skal I ikke lave formelle test, men blot argumentere ud fra jeres "mavefornemmelse".
  1. I en undersøgelse af sammensætningen af en kalkstensprøve er der skåret en tynd skive ud. Denne er delt i 6 områder. I hvert område lægges et gitternet ned over, og der tælles op i 100 knudepunkter, om en bestemt fossil (foraminifer) forefindes. Dette giver følgende resultat:
    Synes I, at disse data tyder på den samme udbredelse af fossilen i de 6 områder?
  2. To steder i Canada er en jordprøve analyseret, og man har talt op, hvor mange guldkorn der er fundet i prøven. Data er:
    Vil du foretrække det ene sted fremfor det andet, hvis du skal grave efter guld?
  3. I opgave 1.5 nedenfor omtales en undersøgelse, hvor 266 mikrojordskælv er klassificeret efter, om de tidsmæssigt falder sammen med stigende tidevand eller faldende tidevand. Af de 266 skælv er 117 i den stigende fase af tidevandet. Er dette udtryk for, at der ikke er en sammenhæng mellem disse mikrojordskælv og tidevandet?
De tre eksempler ovenfor vedrører alle tælledata. Dette er også hovedemnet for øvelserne hørende til de tre første kapitler af denne bog.

Opgave 1.2: Regne i R

Når I skal køre R på jeres egen computer, bør I som hovedregel ikke skrive jeres kommandoer direkte i kommandovinduet. I stedet skal I åbne en editor, skrive kommandoerne her og derefter overføre dem til kommandovinduet. På denne måde er det nemmere at rette fejl og at gentage beregninger med forskelligt input. I windowsversionen af R er der en indbygget editor, hvor I blot skal bruge control R for at overføre markeret tekst (i Mac-versionen skal man i nogle opsætninger bruge command enter i stedet).
Som forberedelse til denne opgave skal I have læst afsnit 1.2, der indeholder en introduktion til R.
  1. Definer en variabel med værdien 4. Udregn kvadratroden af logaritmen af (naturlige logaritme) og eksponentialfunktionen af
    Lad endvidere være en variabel med værdien 10 og en variabel med værdien 1.96. Beregn
    Hvad tror du, at kommandoen n+c(-1,1)*x giver?
  2. Lad være en vektor med indgangene og Beregn og ( er den naturlige logaritme). Numerisk værdi af et tal beregnes i R med funktionen abs.
  3. Lad være en vektor med indgangene og og lad være en vektor med indgangene og Beregn og
  4. Lad være vektoren med indgangene og lad være vektoren med indgangene Beregn og I den sidste sum skal der kun medtages de elementer i hvor det tilhørende element i er lig med 2.
    Lad nu være vektoren med indgangene "Lise", "Lise", "Peter", "Lise", "Peter", "Peter". Beregn
    Prøv at gætte på resultatet af beregningen sum(x[-c(4:6)]).

Opgave 1.3: Figur i R

Denne opgave går ud på, at I skal prøve at lave en figur i R, som beskrevet i afsnit R.5. Til dette skal I bruge data hentet på hjemmesiden quakesearch.geonet.org.nz omkring jordskælv i New Zealand. Tabellen nedenfor indeholder antallet af jordskælv i et givet styrkeinterval (Richterskalaen) i perioden 1930-2015. Styrkeintervallet er i tabellen givet ved midtpunktet og har en bredde på 0.2.
  1. Dan en vektor styrk med midtpunkterne i styrkeintervallerne. Disse værdier er på Richterskalaen, som er en skala af udsving på en seismograf. Dan også vektoren Antal med antallene i de 10 styrkeintervaller.
  2. Benyt kommandoen plot(styrk,Antal) til at lave en figur, hvor antal jordskælv afsættes mod styrken. Prøv at køre kommandoen igen, hvor du indsætter xlab="Richterskala" efter Antal i kaldet af plot.
    Tilføj også en titel til andenaksen i figuren ved at tildele ylab en værdi i kaldet af plot. En overskrift til figuren kan opnås ved at indsætte main="overskrift" i kaldet af plot.
  3. Gutenberg-Richter loven siger, at antallet af jordskælv af en given styrke er eksponentialfunktionen til en lineær funktion af styrken, hvor er styrken. Den bedste tilpasning til data ovenfor fås med og For at indtegne denne kurve i jeres figur kan I benytte kommandoen
    curve(exp(19.97-2.61*x),from=5,to=7,add=TRUE)
    Funktionen curve skal som input have et funktionsudtryk i variablen et startpunkt og et slutpunkt. Tilføjelsen add=TRUE gør, at kurven indtegnes i den allerede eksisterende figur.
    Gentag plot-kommandoen, og prøv at tilføje col=2 til kaldet af curve. Prøv også at tilføje lty=3 til kaldet af curve.

Opgave 1.4: Indlæse datafil i R (dataframe)

I denne opgave skal I prøve at indlæse data fra en fil. I filen JordskaelvDag.csv ligger det observerede antal jordskælv for perioden 1999-2008 fordelt på ugens syv dage. Der er tre versioner, nemlig alle jordskælv over 4 på Richterskalaen, alle jordskælv over 5 på Richterskalaen og en renset version af sidstnævnte, hvor efterskælv er fjernet. Filen har tre søjler. I den første søjle angives hvilken version, der betragtes, kodet som over4, over5 eller renset. I den anden søjle angives ugedagen som og i den tredje søjle står antallet af jordskælv. Data er fundet på hjemmesiden earthquake.usgs.gov.
  1. Indlæs filen, og placer indholdet i Jordskaelv, med kommandoen
    Jordskaelv=read.csv("JordskaelvDag.csv",header=TRUE,stringsAsFactors=TRUE)
    Prøv at skrive head(Jordskaelv) for at se strukturen af denne. Kommandoen head giver de første få rækker i Jordskaelv. Skriv dernæst class(Jordskaelv).
  2. Kommando class fortæller jer, at Jordskaelv er en dataframe. En dataframe er en samling af søjler, der alle har den samme længde, og hvor søjlerne kan være af forskellig type.
    Prøv nu både at skrive Jordskaelv[,3] og JordskaelvAntal. Dette viser, at en søjle kan addresseres på to måder. Prøv også at skrive class(Jordskaelv[,3]).
    Dan en variabel Version med versionsnavnet (søjle 1), en variabel Dag med ugedagen (søjle 2) og en variabel Antal med antal jordskælv (søjle 3).
  3. Lav et datasæt med antal jordskælv for de 7 ugedage for data vedrørende jordskælv med en styrke over 4 ved kommandoen
    AntalOver4=Antal[Vers=="over4"]
    Lav tilsvarende et datasæt med jordskælvsantallene for jordskælv med en styrke over 5 og for rensede jordskælv med en styrke over 5. Find for hvert datasæt summen over de syv ugedage af jordskælvsantallene og gennemsnit per år.
  4. Du skal nu lave en dataframe med dine beregnede værdier ved hjælp af funktionen data.frame. Benyt følgende kommando, hvor du selv må udfylde, hvor der blot står "...".
    Mintabel=data.frame(Sum=...,Gennemsnit=...,
    row.names=c("over4","over5","renset"))
    Skriv dernæst Mintabel for at se resultatet.

Opgave 1.5: Sammenligne observerede med forventede

I artiklen Pulse of the seafloor: Tidal triggering of microearthquakes at 50N East Pacific rise undersøges, om tidspunkter for mikrojordskælv (-0.4 til 2.0 på Richterskalaen) korrelerer med fasen af tidevand (den kombinerede effekt af ocean tidal loading og solid Earth tide). Figuren nedenfor er en skitseagtig fremstilling af figur 2 i artiklen, og viser tidevandsspændingen i 16 dage startende fra 15.december 2003, med tidspunkterne for jordskælv markeret som punkter.
Fra denne figur er antallet af jordskælv i den opadgående og nedadgående fase talt op og angivet i den følgende tabel.
  1. Lad os tænke på dette, som at 266 jordskælv fordeles på de to perioder opadgående og nedadgående. Hvis jordskælvene kommer på tilfældige tidspunkter uafhængigt af tidevandet, vil vi forvente, at de 266 jordskælv fordeles tilfældigt på de to tidevandsfaser. Udtryk dette som en sandsynlighed, for at et jordskælv kommer i den opadgående fase, og angiv de forventede antal svarende til de tomme pladser i ovenstående tabel. Hvad er din umiddelbare fortolkning: tyder data på, at jordskælvene er fordelt ligeligt på de to tidevandsfaser?
For at få en idé om hvorvidt det observerede antal på 117 ud af 266 jordskælv, der forekommer i den opadgående fase, er typisk eller ekstremt under antagelsen om en tilfældig fordeling på de to faser, skal I nu sammenligne udfaldet med en simulering. Vi kan bruge en simulering som den vist i kodevinduet i afsnit 1.1 (med 580 erstattet af 266 og c(0.75,0.25) erstattet af c(0.5,0.5)) til at simulere 266 kast med en mønt (terning med 2 sider), og dernæst tælle op antallet af gange, hvor terningen viser 1, som vi tolker som antallet af gange, at jordskælvet falder i den opadgående fase. Vi kan dog gøre dette nemmere, idet vi ved, at det stokastiske antal af jordskælv, der falder i den opadgående fase, er binomialfordelt. I R kan man simulere fra binomialfordelingen med rbinom:
rbinom(Nsim,266,0.5)
Denne kommando giver Nsim observationer fra en binomialfordeling med antalsværdi 266 og sandsynlighedsparameter
  1. Prøv at køre ovenstående kommando et par gange med Kør dernæst kommandoen med
    Lav en tabel ved hjælp af funktionen data.frame med tre søjler. Første søjle skal indeholde en vektor med de 20 simulerede værdier fra rbinom. Den anden søjle skal indeholde absolutværdien af afstanden mellem det simulerede antal og det forventede antal, Funktionen abs i R giver absolutværdien af et tal. Den tredje søjle skal være enten 0 eller 1, alt efter om er under eller er større end eller lig med denne værdi. Her er afstanden mellem det observerede antal på 117 og det forventede antal. Du kan benytte funktionen ifelse til at lave den tredje søjle: ifelse(D<16,0,1).
  2. Hvor mange af dine 20 simuleringer (som er baseret på en tilfældig placering i forhold til tidevandsfasen) giver anledning til en lige så stor eller større afstand end den, der er observeret i vores faktiske eksperiment? Argumenter for, om observationen på 117 i den opadgående fase er typisk eller atypisk, for hvad man observerer ved en tilfældig fordeling på de to faser.
    Hvis du ikke allerede har gjort det, så læs nu Definition 1.1.1 på en -værdi, og se Resultat 1.4.1 omkring hvordan en -værdi bruges.

ForegåendeNæste