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.
I denne opgave skal I ikke lave formelle test, men blot argumentere ud fra
jeres "mavefornemmelse".
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 ?
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 ?
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.
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.
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 ?
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.
Lad være en vektor med indgangene
og og lad være en vektor med indgangene
og
Beregn
og
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)]).
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.
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.
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.
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.
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.
Indlæs filen, og placer indholdet i Jordskaelv, med kommandoen
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).
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).
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.
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 "...".
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.
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
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).
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.