Download Report

Lunds univrsitet
Matematikcentrum
Matematisk statistik
Biostatistisk grundkurs, MASB11
Laboration 2 VT-2015, 150205
Felrisker – Fördelningar och Simulering
Introduktion
Syftet med laborationen är dels att vi skall bekanta oss med lite av de olika funktioner som
finns i R vad det gäller simulering och dels att öka förståelsen för vissa grundläggande områden
inom sannolikhetsteorin t.ex. frekvenstolkning av sannolikhetsbegreppet, slumpmässiga urval,
slumpvariabel, sannolikhetsfördelning.
• Laborationen redovisas med en kort skriftlig rapport som ska vara inlämnad senast innan
nästa laborationstillfälle tisdag 150217. Rapporten ska omfatta avsnitt 1, 2 och 3 och
vara en syntes av vad ni kommit fram till. Förberedelsuppgifterna ska också lämnas in
– handskrivna från var och en.
Förberedelseuppgifter
1. Räkna uppgift 3.14 i kursboken. Tilläggsfrågor: Beräkna den intressanta sannolikheten
att en patient har sjukdomen om testet är negativt. Vilken egenskap hos testet skall
man försöka ändra för att denna sannolikhet skall bli mindre? Skall man försöka få
sannolikheten för positivt falskt svar att bli 0 eller sannolikheten för sant negativt svar
att bli 1?
2. Räkna uppgifterna 4.8 och 4.9 i kursboken.
3. I en befolkning är 20% rökare och vi väljer slumpmässigt ut 5 personer.
(a) Vad är sannolikheten att ingen av de fem är cigarrettrökare?
(b) Vad är sannolikheten att alla fem är cigarrettrökare?
(c) Vad är sannolikheten att minst 2 är cigarrettrökare?
(d) Vad är väntevärde och varians för antalet cigarrettrökare?
1
Diagnostik, felrisker och hypotesprövning
Vi skall börja med att studera några aspekter av det sjukdomstest som presenterades i
övningsuppgift 3.14, och som du bör ha arbetat dig igenom för att få fullt utbyte av detta avsnitt.
Det finns ett färdigt R-script, Uppg314.R på kurshemsidan, som du kan använda dig av. Undersök vad funktionen gör och vilka inparametrar som behövs. Anropa därefter först funktionen med de värden som anges som standard i scriptet. ändra sedan sannolikheten för positivt
utslag om personen har sjukdomen till 0.999 och låt funktionen rita nya kurvor. ändra slutligen
sannolikheten för positivt utslag om personen inte har sjukdomen till 0.01 och låt funktionen
rita nya kurvor. Studera de tre figurerna och besvara följande frågor:
• Identifiera de fyra kurvorna i varje figur, vilka sannolikheter representerar de?
• Kurvorna hör ihop parvis, på så sätt att den ena enkelt kan rekonstrueras ur den andra.
Förklara hur.
• Diskutera vilka egenskaper du som patient skulle värdera högst vid ett test av detta slag
och relatera det till kurvorna i figurerna.
Exemplet med sjukdomstestet illustrerar väl de felrisker man måste ta i beaktande, när man
skall konstruera ett hypotestest. Om vi, rent allmänt, ställer upp en hypotes H0 som vi vill
pröva, kan vi hamna i någon av de situationer som beskrivs i figuren nedan.
H0 förkastas
H0 förkastas ej
H0 sann
fel typ I
OK
H0 falsk
OK
fel typ II
Hur teststorheten kommer att se ut och hur felriskerna skall beräknas beror på hur nollhypotesen H0 formulerats. Nollhypotesen skall först och främst väljas på ett sådant sätt att ett
förkastande av nollhypotesen ger ett tydligt svar på den fråga, man söker svar på.
• Utgå från din diskussion i sista frågepunkten ovan och formulera (i ord) en nollhypotes,
som du finner adekvat. Identifiera sedan de två felriskerna P [fel typ I] och P [fel typ II]
och tala om hur stora de blir i ditt fall. (Använd de siffror som ges i uppgift 3.14.)
2
Sannolikhet enligt frekvenstolkningen — Kast med tärning
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
2.1
Uppvärmning
Om vi kastar en symmetrisk tärning förväntar vi oss, att i det långa loppet skall alla sex
sidorna ha kommit upp ungefär lika många gånger. Detta betyder till exempel att om vi
räknar antalet treor som kommit upp, så bör detta antal utgöra ungefär en sjättedel av det
totala antalet kast.
Eftersom frekvenstolkningen handlar om vad som händer i det långa loppet — i vårt exempel
efter många tärningskast — kan det ju lätt hända att man tröttnar och avbryter sitt försök
2
innan man hunnit skaffa tillräckligt mycket data (att kasta en tärning tiotusen gånger kan ju
bli lite jobbigt). Det finns dock pionjärer som offrat sig och verkligen gjort detta. För att visa
vad som händer, utan att bevisa någonting, räcker det att utföra en datorsimulering, det vill
säga, låta en dator utföra försöket i stället.
Vi skall nu med datorns hjälp simulera hundra tärningskast och studera den relativa frekvensen
av treor.
Simulera nu de hundra tärningskasten. Glöm inte att avsluta kommandot med semikolon,
annars kommer skärmen att fyllas med hundra stycken slumptal.
>> X <- floor(6*runif(100,0,1)+1);
Funktionen floor avrundar nedåt. Vi vill nu kontrollera att elementen i X verkligen har en
fördelning som ett tärningskast.
Första steget blir att räkna antalet treor. I R finns ett antal relationsoperatorer, som tillåter
jämförelser av matriser. Med kommandot
>> Y <- X==3;
får vi en vektor av samma dimension som X och som enbart innehåller ettor och nollor. På
varje plats där X har en trea, har Y en etta, och på varje plats där X har ett element som
inte är en trea, har Y en nolla. Genom att räkna antalet ettor i vektorn Y, får vi alltså reda
på hur många treor, som finns i vektorn X.
De successiva relativa frekvenserna av treor kan vi nu beräkna med följande kommando:
>> relf <- cumsum(Y)/seq(1,100)
Funktionen cumsum ger en vektor där element nummer i är summan av de i första elementen
i inparametern, i vårt fall Y. Notationen seq(1,100) är en vektor med talen 1 till och med 100.
övertyga dig om att relf innehåller de successiva relativa frekvenserna. Nu kan vi plotta de
relativa frekvenserna:
>> plot(relf)
2.2
Differensen mellan antalet treor och fyror
Nu ska analysen kompletteras genom att undersöka hur den successiva skillnaden av antalet
treor och fyror ser ut som funktion av antalet kast. Börja med att skapa en ny variabel Z
som innehåller den kumulerade summan av antalet fyror och beräkna därefter den successiva
skillnaden mellan antalet treor och fyror. öppna därefter ett nytt grafikfönster och plotta
skillnaden mellan antalet treor och antalet fyror som funktion av antalet kast:
>>
>>
>>
>>
win.graph()
Z <- X==4;
diff <- cumsum(Y)-cumsum(Z);
plot(diff)
• Ta en liten paus och fundera över följande fråga: Törs du dra några slutsatser eller
eventuellt ställa upp något antagande utifrån kurvorna över de relativa och absoluta
frekvenserna? Vilka resultat förväntar du dig rent teoretiskt vid många kast med en
symmetrisk tärning?
Hundra kast är kanske lite för lite. Simulera istället tiotusen tärningskast och beräkna relativa
frekvensen treor som funktion av antalet kast på samma sätt som ovan (spara gärna resultatet
3
i en ny variabel, till exempel relf1). Kan du säga något om den relativa frekvensen treor, när
antalet kast är stort?
Beräkna också skillnaden mellan antalet treor och antalet fyror och spara den tex. i diff1.
Kan du säga något om skillnaden när antalet kast är stort?
För att bättre kunna se vad som händer i försökets början, kan man plotta relf1 mot log(antal kast):
>> plot(log10(seq(1,10000)),relf1)
Om man använder kommandot points(...) efter att man skapat en figur med plot(...) kan
man låsa figuren i bildfönstret, så att man kan rita nya kurvor ovanpå.
Gör ytterligare en försöksserie med tiotusen kast. Analysera och plotta relativa frekvensen
treor (som du sparat i en ny variabel, till exempel relf2) gentemot log(antal kast) på samma
sätt som ovan. För att kunna skilja de två kurvorna åt, kan man plotta den andra kurvan med
en annan signatur eller färg, till exempel
>> plot(log10(seq(1,10000)),relf2,col ="red"’)
Notera såväl skillnader som likheter mellan de två försöksserierna.
Rita också upp differenserna mellan treor och fyror för de två försöksserierna och notera
skillnader och likheter.
2.3
Analys av försöket kasta en tärning
För att belysa vad som sker vid det här experimentet (kast med tärning) skall du utföra
försöket flera gånger under likartade förhållanden. Det vi kallar olika realiseringar. Försöket
skall utföras cirka 7 gånger.
För att utröna vad som sker i långa loppet behövs det ett stort antal kast, helst 100 000, i
varje realisering. Bilden blir klarare ju fler kast du använder; numera bör internminnet hos
alla datorerna räcka för så stora simuleringar, vid varje realisering skall du i varje fall använda
lika många kast, dock minst 10 000. Följande frågor skall besvaras, varje fråga skall åtföljas av
en lämplig figur där resultatet av de olika realiseringarna framgår. Man kan med fördel utföra
beräkningarna med en så kallad script-fil i R. Se Introduktion till R.
• Hur många kast gissar du behövas för att vi säkert skall kunna uppskatta sannolikheten
för en trea respektive en fyra? Ett motiverat svar kan inte ger förrän längre fram i kursen.
• Vilka slutsatser kan du dra angående skillnaden av de relativa frekvenserna av treor
och fyrar i en realisering? är resultatet i de olika realiseringarna samstämmigt? Vilket
resultat förväntar du dig rent teoretiskt vid många kast med en symmetrisk tärning?
• Är kurvan över skillnaden mellan antalet treor och antalet fyror förenlig med din slutsats
från föregående punkt? Förklara också varför resultatet blev som det blev.
3
Binomialfördelningen
Om man till exempel vill skapa ett antal slumptal (25 stycken) från en binomialfördelning där
antale försök n=10 och sannolikheten för den lyckade händelsen a är 0.2 skriver man följande
kommando:
xbino1 <- rbinom(25,10,0.2)
4
I variabeln xbino1 lagras de simulerade värdena. För att kolla på resultatet kan man skriva
bino1 eller print(bino1) i kommandofönstret. Det som bland annat är intressant att ta reda
på i stickprovet är om sannolikhetsfördelningen, den relativa frekvensen av de olika utfallen
stämmer överens med den valda binomialfördelningen. Vi måste alltså räkna efter hur många
händelser av respekive utfall som finns i stickprovet. I R kan man göra detta genom att först
kategorisera resultatet i variabeln bino1 och därefter beräkna den relativa frekvensen:
xomega <- seq(0,10,1) #utfallsrummet
xkategori <- factor(xbino1) # kategorisering av resultatet i xbino1.
xfreq <- table(xkategori) # r\"{a} absoluta frekvensen.
relfreq <- prop.table(xfreq) # skattningen av sannolikhetsfunktionen.
Kontrollera nu vilka av de möjliga utfallen som kom med stickprovet. Kom alla med? För att
kunna plotta den skattade sannolikhetsfunktionen mot respektive utfall måste man skapa en
vektor – variabel som innehåller dessa värden. Gör nu detta:
xutfall<-c(?,?,?,...,?,?)
När man gjort detta kan man plotta den skattade sannolikhetsfunktionen i ett stolpdiagram
i R med hjälp av plot-funktionen:
plot(xutfall, relfreq, type="h",col="blue")
# alternativt med barplot
Den teoretiska sannolikhetsfunktionen pX (x) för motsvarande binomialfördelade slumpvariabel, s.v. X, kan man enkelt plocka fram i R med hjälp av R-kommandot dbinom(x,n,p). Vill
man ha tag på fördelningsfunktionen FX (x) använder man R-kommandot x,n,p. I vårt fall
skriver man tex:
px <- dbinom(x,10,0.2) # x={0,1,...,10} utfallsrummet.
points(xomega, px, type="h", col="red") # plottar i samma figur.
• Undersök nu hur många stickprov frå en s.v.X ∈ Bin(10, 0.2) som behövs för att hitta en
rimlig överensstämmelse mellan den skattade sannolikhetsfunktionen och den teoretiska
sannolikhetsfunktionen pX (x). Gäller samma sak för den s.v. X ∈ Bin(10, 0.5) eller för
den s.v. X ∈ Bin(25, 0.7)?
4
R-script Uppg314.R
ps<-0.9
pf<-0.05
p<-seq(0,1,0.001)
ppos<-ps*p+pf*(1-p)
p7kpos<-(ps*p)/ppos
p7kneg<-(1-ps)*p/(1-ppos)
win.graph()
plot(p,p7kpos,type="l",col="red",xlab="p, relativsjukdomsfrekvens",ylab="P, Sannolikhet")
points(p,1-p7kpos,type="l",col="blue")
points(p,p7kneg,type="l",col="dark red")
points(p,1-p7kneg,type="l",col="green")
grid()
title(main=\textrm{"Konfidens f\"{o}r sjukdomstest"})
5