Download Report

Lunds univrsitet
Matematikcentrum
Matematisk statistik
Biostatistisk grundkurs, MASB11
Laboration 3 HT-2014, 141219
Fördelningsanpassning och Centrala Gränsvärdes Satsen
Introduktion
Syftet med laborationen är dels att vi skall bekanta oss med vad som sker med sannolikhetsfördelning när man bildar summor av slumpvariabler (Centrala gränsvärdessatsen) och
dels med några metoder och funktioner som finns i R för att kunna skatta, studera och jämföra
fördelningen hos ett eller flera stickprov, empirisk fördelningsfunktion, F ∗ (x) och olika typer av
fördelningsplottar. Datamaterialet och script finns på kursens material-hemsida under namnet
www.maths.lth.se/matstat/kurser/masb11/htm4/
• Laborationen redovisas med en skriftlig rapport som ska vara inlämnad senast 150109.
Simulering av slumpvariabler i R
Simulering i R görs genom färdiga funktioner unika för respektive fördelning. Exempel på
några av dessa funktioner finns i efterföljande tabell.
Fördelning
Binomial
Poisson
Normal
Likformig (Rektangel)
Exponential
Funktion i R
rbinom(antal,n,p)
rpois(antal,mean)
rnorm(antal,mean,stddev)
runif(antal,min,max)
rexp(antal,scale) (yr)
Exempel i R
rbinom(100,10,0.5)
rpois(100,2)
rnorm(100,10,15)
runif(100,10,20)
rexp(100,1)
mean = väntevärdet, µ, i fördelningen
stddev = standardavvikelsen, σ, i fördelningen
scale = 1/väntevärdet i exponentialfördelningen
1
Fortsättning laboration 2
Uppgift 5
I de tidigare uppgifterna har vi simulerat vad som händer om vi tar stickprov av olika storlekar
från olika kända fördelningar. Vi skall nu gå vidare och undersöka vad som händer om vi
bildar olika storheter i stickprovet. Vilka egenskaper får då dessa storheter? Det är framförallt
väntevärdet E[X], variansen V [X] och fördelningen FX (x) som vi intresserar oss för.
Vi börjar med att undersöka vilken fördelningen summan av två observationer från en normalfördelning med väntevärde 10 och standardavvikelse 2 har. Börja med generera två nya
stickprov om 1000 observationer norm1 och norm2:
norm1 <- rnorm(1000,10,2)
norm2 <- rnorm(1000,10,2)
När vi kör dessa kommandon kommer det att bildas två nya variabler som heter norm1 och
norm2 och som innehåller 1000 slumptal var.
Bilda nu summan (sum12) av de två kolumnerna norm1 och norm2. Undersök vilken fördelning
summan har genom att göra ett histogram och en Q-Q plot.
sum12 <- norm1+norm2
Vilken fördelning har summan? Vad bör väntevärdet bli? Standardavvikelsen? (Använd gärna
√
R, x beräknas med sqrt(x)) Beräkna också medelvärdet av sum12 med mean(sum12) och
stickprovsstandardavvikelsen sd(sum12). Hur passar de med de teoretiska värdena?
Centrala gränsvärdessatsen
Lägger man ihop, adderar, (eller beräknar medelvärdet) av flera oberoende normalfördelade
slumpvariabler är summan också normalfördelad. Men vad händer om man lägger ihop flera
variabler som alla är rektangelfördelade? Vilken fördelning fås om man adderar exponentialfördelade variabler?
Centrala gränsvärdessatsen säger att om man adderar ett stort antal oberoende variabler
från en godtycklig fördelning blir summan (eller medelvärdet) normalfördelad. Detta märkliga
faktum ska du i denna uppgift undersöka med hjälp av den interaktiva rutinen cgs().
Konkret kan vi tänka oss att du gör ett antal mätningar av en intressant (bio)variabel, du
bildar summan av mätningarna (eller medelvärdet). Det du ska undersöka är hur summan
kommer att variera från mätserie till mätserie? Beror det på ursprungsfördelningen hos den
uppmätta variabeln?
Så här använder du rutinen cgs i RStudio
• När du skriver cgs() får du möjlighet att välja mellan ett antal fördelningar med givna
parametrar eller kan du konstruera en egen diskret sannolikhetsfördelning. Välj ett av
alternativen genom att mata in tillhörande siffra.
• Du får en figur med täthetsfunktion eller sannolikhetsfunktion för din valda fördelning.
Välj nu hur många mätningar du ska göra från denna fördelning och mata in detta antal.
• I kommandofönstret visas resultatet av din mätningar (de 10 första om du valt ett
stort antal), d.v.s. R har hämtat slumptal från din valda fördelning. Summan av alla
mätningarna skrivs ut. I din figur markeras mätningarna med kryss.
• Antag nu att du gör upprepade serier med det antal mätningar, n, som du valt. För
varje serie beräknas summan av dina mätningar. Hur varierar då summan? Mera matematiskt beskrivet: Om X1 , X2 , . . . , Xn är oberoende med den fördelning du valt, vad är
då fördelningen för summan X1 + X2 + . . . + Xn ?
• Undersök detta genom att simulera N serier med det antal mätningar (n) du valt. Rutinen ritar sedan upp ett histogram för summan. Ange alltså ett värde på N, tänk på att
välja N tillräckligt stort så att du kan få en uppfattning av fördelningen i histogrammet.
2
• Centrala gränsvärdessatsen säger att om du valt ett tillräckligt stort antal mätningar
kommer fördelningen för summan att bli ungefär normalfördelning. Rutinen ger dig
möjlighet att anpassa en normalfördelning till data. Du kan låta R sköta om det och
din uppgift blir då att undersöka grafiskt om du tycker att approximationen verkar bra.
Till din hjälp har du också en Q-Q plot där du kan se om summan verkar passa till en
normalfördelning.
• Du kan också anpassa normalfördelningen själv och måste då fundera på vilka värden på
väntevärde och standardavvikelse som gäller (prova gärna detta själv som en utmärkt
övning!).
• Om du vill köra rutinen igen kan du undvika den interaktiva fasen genom att direkt
skriva in dina val i anropet. Exempelvis ger cgs(2,10,1000,1) att 1000 serier med
vardera 10 mätningar slumpas från en likformig fördelning, R(0, 4). Histogramet för de
1000 summorna plottas, normalfördelning anpassas och en Q-Q plot ritas.
Uppgift 6
1. Välj rektangelfördelning, antal=2 i R-funktionen cgs(). Vilka värden kan summan av
två mätningar ligga mellan? Verkar histogrammet rimligt?
2. Öka antalet mätningar i rektangelfördelningen. Vad händer om du tar antalet mätningar
till 5? Eller ökar till 10?
3. Försök anpassa ”rätt” normalfördelning till histogrammet, d.v.s. tänk ut värdena på
väntevärde och standardavvikelse.
(b−a)2
Ledning: Om s.v. X är uniform(a,b): ⇒ E[X] = a+b
2 och att V [X] =
12 .
4. Exponentialfördelning: Gör nu motsvarande för exponentialfördelningen. Hur många
mätningar behöver ni ta innan ni tycker att summan är ungefär normalfördelat? Verkar
fördelningen gå snabbare eller långsammare mot en normalfördelning än det gjorde för
den likformiga fördelningen. Vad beror detta på?
5. Normalfördelning: Vad händer om ni tar antal=2? Kan du förklara detta?
6. Undersök gärna på motsvarande sätt vad som händer då man bildar summor från
binomial- eller poissonfördelningen.
7. Testa gärna med en egen diskret fördelning, tex binomialfördelningen.
8. Du har tittat på vad som händer med summor av variabler. Vad händer om man i stället
tar medelvärdet av variablerna (mätningarna)?
Uppgift 7 (Har ni gjort den under labpass 2 behöver ni inte göra om den)
Centrala gränsvärdessatsen i praktiken: På 35 patienter med Hodgkins sjukdom mätte man
antalet T4 celler i blodet (antal/mm3 ). Samtidigt mätte man motsvarande antal hos 35 patienter som hade andra sjukdomar (Non-Hodgkins). Data ligger i filen Hodgkindata.RData
som du hittar på kursens hemsida. Läs in data via Workspace-fönstrets öppna-ikon. Du har
nu fått två nya variabler Hodgkin och NonHodgkin. Undersök om antalet celler i blodet
är normalfördelat för de båda grupperna.
3
Det är möjligt jämföra grupperna genom att bilda differensen mellan de två gruppmedelvärdena.
Kan du använda dig av centrala gränsvärdessatsen i detta fall? Kan du säga något om vilken
fördelning differensen i medelvärden har? är det ett stort problem att variabeln inte är normalfördelad i de båda grupperna från början? Kan man åtgärda detta på något sätt? Prova
och se vad som händer med
√ fördelningen för data om man istället använder en transformation
av värdena, till exempel X eller log X.
2
Fördelningsanpassning
2.1
Empiriska fördelningsfunktionen F ∗ (x) – normalfördelningspapper
Grafiska metoder används främst för tre ändamäl: skattning av parametrar, validering av
fördelning samt skattning av kvantiler. Den grafiska tekniken bygger kort på att, vid givet
slumpmässigt stickprov x1 , x2 ,...,xn :
1. Först ordnas stickprovet, betecknas x(1) , x(2) ,...,x(n) .
2. Man skattar fördelningsfunktionen F (x) med det vi kallar för den empiriska fördelningsfunktionen F ∗ (x). Den definieras som:
F ∗ (x)


 0
; x < x(1)
i/n ; x(i) ≤ x < x(i+1)
=

 1
; x(n) ≤ x
3. Därefter plottas de n stycken talparen (x(i) , ( ni )). Plottningspositionen ni som vi använder
för den empiriska fördelningsfunktionen har en del fördelar men också vissa nackdelar,
tex. att x(n) kommer att vara den position som svarar mot 1 hos fördelningsfunktionen.
i
Andra val av plottningspositioner förekommer: exempelvis, (x(i) , ( n+1
) eller (x(i) , ( i−1/2
n ),
se Holmquist B. Matematisk statistik för M och V, Kompletteringar och tillämpningar,
1996.
I R kan F ∗ (x) ritas med hjälp av funktionen plot(...,type="s",...). Nedansåtende kommandorader exemplifierar teckniken med hjälp av 100 observationer från en s.v. X ∈ N (2, 1).
>
>
>
>
>
X<-rnorm(100,2,1)
sortX<-sort(X)
Fn<-seq(1,length(sortX),1)/(length(sortX)+1) # plotposition foer
plot(sortX,Fn,type="s",col="blue")
grid()
F_n
• På y-axeln har vi F ∗ (x). Använd denna för att skatta medelvärdet, kvartilerna samt
medianen i fördelningen.
• Eftersom vi känner µ och σ i det här fallet kan vi komplettera figuren med den riktiga
fördelningsfunktionen, FX (x). Gör det, glöm inte att använda points(...) istället för
plot(...) innan du plottar ovanpå F ∗ (x).
4
2.2
Normplot i R-package nsRFA
Om vi vet eller misstänker att stickprovet kommer från en normalfördelning kan vi istället
plotta det ordnade stickprovet i ett normalfördelningspapper. Skalan på y-axeln i ett normalfördelningspapper är anpassad så att observationerna kommer att följa en rät linje om de
är normalfördelade. Jämför teknken med qqnorm. Om vi får någon kurvatur indikerar detta
alltså att observationerna inte är normalfördelade. Om man har installerat ett packages som
heter nsRFA i R kan man direkt plotta ett eller flera stickprov i ett normalfördelningspapper
med kommandot normplot() och normpoints(). Använd ?normplot för att komma underfund med funktionen. Plotta därefter stickprovet X i ett normalfördelningspapper.
> ?normplot
> normplot(X)
Uppgift
Skatta nu medelvärdet µ och standardavvikelsen σ i normalfördelningsplotten, skattningstekniken är en direkt tillämpning av moment 2 i Lektionsblad 4. Stämmer skattningarna med det
använda stickprovet?
Om R-packages nsRFA inte redan finns installerat i din dator kan man enkelt göra det i
RStudio genom att välja flicken packages i nedre högra fönstret och därefter kryssa för de
packages man vill installera. Det finns en hel del olika fördelnnigspapper man kan ha nytta av
i nsRFA.
2.3
Jordbävningar
Vi ska nu studera ett datamaterial där data insamlats under perioden den 16 december
1902 t.o.m. den 4 mars 1977. Det rör sig om tidsintervall, mätt i dagar, mellan kraftiga
jordbävningar världen runt. Jordbävningar med en magnitud på åtminstone 7,5 på Richterskalan finns representerade, alternativt jordbävningar med över 1000 dödsoffer.
Datamaterialet finns på kursens material-hemsida under namnet Quakes.RData. Läs in filen
genom kommandot load("Quakes.RData"). De numeriska värdena finns lagrade i en vektor
med namn quakeper. Använd length för att finna antalet tidsperioder. Som tidigare ritar vi
histogram och beräknar diverse läges- och spridningsmått:
>
>
>
>
>
>
hist(quakeper,freq=FALSE) # histogrammets totala area blir 1, taethetsfuktion#
m<-mean(quakeper)
med<- median(quakeper)
s<-sd(quakeper)
s2<- var(quakeper)
myrange<-range(quakeper)
Uppgift
Använd data och fundera: verkar det troligt att det kan gå längre period än 5 år mellan
kraftiga jordbävningar? Hur ofta inträffar de?
I R finns en del användbara villkorssatser som gör det enkelt att skapa nya vektorer och
matriser med hjälp av lämpliga bivillkor, Man kan alltså på så sätt i en given vektor eller
5
matris finna element som uppfyller ett aller annat intressant villkor. För att exmpelvis finna
de perioder mellan jordbävningar med längd kortare än 1000 dagar (c:a 3 år) kan man skriva:
> less1000 <- quakeper[quakeper < 1000];
> length(less1000)
Första kommandot skapar en vektor som vi kan kalla vad som helst, tex. less1000. Den
innehåller de element i ursprungsvektorn quakeper vilka uppfyller villkoret. För att få reda
på hur många element som uppfyller villkoret använder vi helt enkelt length(less1000).
Uppgift
Vi vill uppskatta sannolikheten för att en period mellan jordbävningar är kortare än 1000
dagar genom att beräkna motsvarande andel i datamaterialet Vi har i själva verket beräknat
täljaren fall i kommandoserien ovan, och nämnaren ges helt enkelt av length(quakeper).
Beräkna nu den intressanta kvoten och notera ditt svar. Hur stor är sannolikheten att det
dröer mer än 200 dagar mellan två stora jordbäningar?
Anmärkning. Den storhet som beräknades som mean(quakeper) benämnes iblandåterkomsttid (engelska: return period), beteckna den tex. med Tr , detta är egentligen en skattning av
vänteärdet. Intensiteten av de händelser som studeras kan beräknas som 1/T och studeras
ofta i statistisk riskanalys, brukar betecknas med λ.
2.4
Anpassning till standardfördelning
Enligt gängse statistisk teori är tidsavståndet för två händelser som uppträder slumpmässigt
i tiden efter varandra exponentialfördelade, s.v. T ∈ Exp(λ). (Se sid 98 i kurboken). Dess
Fördelningsfunktion kan då skrivas som FT (t) = 1 − e−λ·t , t > 0.
Uppgift
Avgör nu om tidsavståndet mellan två efterföljande jordbävningar är exponentialfördelad
genom att jämföra den empiriska fördelningsfunktionen, F ∗ (t) , med den teoretiska, FT (t)
R-tips:
>
>
>
>
>
>
>
sortX<-sort(quakeper)
Fn<-seq(1,length(sortX),1)/(length(sortX)+1) # plotposition foer
plot(sortX,Fn,type="s",col="blue")
grid()
taxis<-seg(0.01,max(quakeper,1))
FT<- 1-exp(-\lambda*taxis)
points(taxis,FT,type="l",col="red")
F_n
Uppgift Hodgkin-NonHodgkin
För att avgöra om två datamaterial, stickprov, är lika ska de ha samma fördelning de ska
ju vara stickprov från samma population. Avgör nu om de två stickproven i datamaterialet
Hodgkindata.RData har samma fördelning med hjälp av dess empiriska fördelningsfunktioner
och normalfördelningspapperet, normplot(). Tips små storheters.v., X, som är positiva men
med värden nära noll kan ofta modelleras med hjälp av en lognormalfördelning, s.v. log X ∈
N (µ, σ 2 ). Stämmer detta?
6