Syftet med den här laborationen är att du skall bli mer förtrogen med följande viktiga områden inom matematisk statistik

FMS012/MASB03: M ATEMATISK STATISTIK , 9 HP , HT-16

Datorövning 3

Hypotesprövning och styrka

Syftet med den här laborationen är att du skall bli mer förtrogen med följande viktiga områden inom matematisk statistik

• Hypotesprövning

• Intervallskattning

Dessutom får du möjlighet att arbeta igenom ett något större verkligt problem. Vi kommer att ägna oss åt statistisk analys av radonmätningar i bostadshus och försöka bedöma om gällande gränsvärden kan anses vara över- eller underskridna.

Specialrutiner finns att hämta på kursens hemsida:

http://www.maths.lu.se/kurshemsida/fms012masb03/

1 Förberedelseuppgifter

1. Läs igenom denna handledning.

2. Förvissa dig om att du förstår hur hypotesprövning går till.

3. Redovisas vid laborationens start! Vi har ett stickprov x 1 , . . . , x 5 från X i ∈ N(μ, σ) där σ = 1 är känd och vi skattar μ med μ ^∗ = ¯ x. Vi vill testa H 0 : μ = 0 mot H 1 : μ 6= 0 på signifikansnivån α = 0.05 med hjälp av en teststorhet.

(a) Skriv upp hur teststorheten ser ut och ange ett villkor för när H ₀ ska förkastas.

(b) Räkna ut väntevärdet för teststorheten när det sanna μ-värdet är μ = 1.

4. Redovisas vid laborationens start! Vi har en observation x = 3 från X ∈ Po(μ) och vill testa H ₀ : μ = 8 mot H 1 : μ < 8 på signifikansnivån α = 0.05 med hjälp av direktmetoden.

Beräkna testets P-värde och avgör om H ₀ ska förkastas eller ej.

2 Styrkefunktion

2.1 Normalfördelning med känt σ

Vi ska undersöka styrkan h(μ) = P(förkasta H 0 om μ är det rätta värdet) hos ett test av nollhypote-

sen H ₀ : μ = μ 0 och se hur den beror på, bl.a., stickprovsstorleken n. Till vår hjälp har vi några speci-

alskrivna funktioner som kan laddas ned från kurshemsidan. Vi börjar med fallet där X i ∈ N(μ, σ)

med känd standardavvikelse σ och μ ^∗ = ¯ x. Ladda ner normstyrka.m från kurshemsidan och se

vad den gör:

>> help normstyrka

>> normstyrka

Den övre figuren visar täthetsfunktionen för teststorheten T = μ ^∗ − μ 0

D(μ ^∗ ) , där D(μ ^∗ ) = σ

√ n , dels när μ = μ 0 och dels när μ = μ 1 . Det röda området är det kritiska området med arean lika med signifikansnivån α. Den blå arean representerar sannolikheten att inte förkasta H 0 när det sanna μ-värdet är μ 1 .

I den undre figuren har vi styrkefunktionen h(μ) = P(förkasta H 0 om μ är det sanna värdet) med h(μ 1 ) markerad.

Uppgift: För vilka värden på teststorheten T ska vi förkasta H 0 i det här fallet?

(Tips: norminv(1-0.05/2) ger λ 0.05/2 ) Jämför med figuren och förberedelseuppgift 3(a).

Uppgift: Var ligger toppen på fördelningen för teststorheten när μ = μ 1 ? Kontrollera att det stäm- mer med förberedelsesuppgift 3(b).

Uppgift: Hur stor är sannolikheten att upptäcka att μ 6= 0 när det sanna μ-värdet är 1?

Uppgift: Hur ändrar sig styrkan när det sanna μ-värdet ökar?

Uppgift: Vad händer med styrkan när det sanna μ-värdet närmar sig μ 0 ?

Uppgift: Experimentera med olika värden på μ 1 (t.ex. normstyrka(’mu1’,2) ger styrkan när μ 1 = 2). Ungefär hur stort (eller litet) behöver det sanna μ-värdet vara för att vi ska ha minst 80 % sannolikhet att upptäcka att μ 6= 0?

Uppgift: Experimentera med olika värden på stickprovsstorleken n (t.ex. normstyrka(’n’,10) ger styrkan när n = 10). Ändrar sig det kritiska området när du ändrar n?

Uppgift: Ändrar sig styrkan h(μ 1 ) när du ändrar n?

Uppgift: Hur stort behöver n vara för att vi ska ha minst 80 % sannolikhet att upptäcka att μ 6= 0

när det sanna μ-värdet är 1? (Pröva dig fram).

Uppgift: Ändra signifikansnivån α till 1 % (normstyrka(’alfa’,0.01)). Vad hände med det kritiska området? Blev det större eller mindre?

Uppgift: Vad hände med styrkan när μ = μ 1 ? Blev den större eller mindre?

Uppgift: Hur stort behöver n nu vara för att vi ska ha minst 80 % sannolikhet att upptäcka att μ 6= 0 när det sanna μ-värdet är 1?

2.2 Normalfördelning med okänt σ

I de flesta praktiska situationer känner vi inte σ utan den måste skattas med σ ^∗ = s. Det gör att teststorheten T = μ ^∗ − μ 0

d (μ ^∗ ) , där d (μ ^∗ ) = s/ √

n blir t(n − 1)-fördelad när H 0 är sann. Den kommer då att variera mer än tidigare, eftersom s i nämnaren också varierar slumpmässigt. Det gör det besvärligare att beräkna styrkan men med hjälp av Matlabs funktioner går det. Ladda ner tstyrka.m från kurshemsidan och se vad den gör:

>> help tstyrka

>> tstyrka

Den övre figuren visar nu täthetsfunktionen för teststorheten T = μ ^∗ − μ 0

d (μ ^∗ ) , där d (μ ^∗ ) = s

√ n , dels när μ = μ 0 och dels när μ = μ 1 . När μ = μ 0 blir det en t(n − 1)-fördelning som är symmetrisk kring 0. När μ = μ 1 får vi istället en icke-central t-fördelning ¹ . Observera att den inte är symmet- risk. Denna besvärliga fördelning är en anledning till att man ofta föredrar att anse att σ är känd när man beräknar styrkan för ett test.

Uppgift: För vilka värden på teststorheten T ska vi förkasta H 0 i det här fallet?

Tips: tinv(1-0.05/2,5-1) ger t 0.05/2 (5 − 1). Jämför med figuren.

Uppgift: Hur stor är nu sannolikheten att upptäcka att μ 6= 0 när det sanna μ-värdet är 1? Är den större eller mindre än när σ var känt?

Uppgift: Experimentera med olika värden på μ 1 (t.ex. tstyrka(’mu1’,2) ger styrkan när μ 1 = 2). Ungefär hur stort (eller litet) behöver det sanna μ-värdet nu vara för att vi ska ha minst 80 % sannolikhet att upptäcka att μ 6= 0?

1

Om X ∈ N(0, 1) och Y ∈ χ

(f ) är oberoende så är X + Δ

pY /f icke-central t-fördelad med f frihetsgrader och centreringsparameter Δ. I vårt fall är f = n − 1 och Δ = μ

− μ

σ/

√ n . Den intresserade kan läsa mer på

https://en.wikipedia.org/wiki/Noncentral_t-distribution.

Uppgift: Experimentera med olika värden på stickprovsstorleken n (t.ex. tstyrka(’n’,10) ger styrkan när n = 10). Ändrar sig det kritiska området när du ändrar n nu? Varför?

Uppgift: Hur stort behöver n nu vara för att vi ska ha minst 80 % sannolikhet att upptäcka att μ 6= 0 när det sanna μ-värdet är 1?

Uppgift: Ändra signifikansnivån α till 1 % (tstyrka(’alfa’,0.01)). Vad hände med det kri- tiska området? Blev det större eller mindre?

Uppgift: Vad hände med styrkan när μ = μ 1 ? Blev den större eller mindre?

Uppgift: Hur stort behöver n nu vara för att vi ska ha minst 80 % sannolikhet att upptäcka att μ 6= 0 när det sanna μ-värdet är 1?

Uppgift: Experimentera vidare med att ändra n, μ 0 , μ 1 , σ och α och se vad som händer. Du kan också göra ensidiga test, se hjälpfunktionen för ett exempel.

3 Radon

3.1 Något om radonmätningar

Radon är en ädelgas som är radioaktiv. Den vanligast förekommande isotopen har en halverings- tid på 3.8 dygn. Radonisotopen sönderfaller till nya ämnen, s.k. radondöttrar, som i sin tur är radioaktiva med mycket kort halveringstid. Vid sönderfallen bildas alfa-partiklar, som, när de far fram, kan orsaka skada i sin allra närmaste omgivning. Om gasen eller någon av döttrarna har in- andats utgör lungvävnaden den närmaste omgivningen. Radon och dess döttrar är delar av en lång sönderfallskedja som startar med uran och slutar med bly.

Ett sätt att mäta radonkoncentrationen i inomhusluften är att hänga upp en alfa-känslig film. När den träffas av en alfa-partikel uppstår en skada i filmen i träffpunkten. Denna skada förstärks vid

”framkallning” av filmen så att det blir ett hål i filmen. Bilden nedan visar hur ett hål kan se ut efter framkallning då man tittar på filmen i mikroskop. Hålen har maximalt diametern 7 μm. Antalet hål på en yta är ett mått på radonkoncentrationen.

3.2 Statistisk modell

För att kunna göra en ordentlig statistisk analys av ett mätmaterial behöver vi mer statistisk kunskap

om radioaktivt sönderfall. Det visar sig att tidpunkterna och platserna (rumskoordinaterna) för

sönderfallen bildar en s.k. poisson-process (efter den franske matematikern Siméon Denis Poisson (1781–1840)). Poisson-processen behandlas utförligt i fortsättningskursen i stokastiska processer.

Enkelt kan man säga att sannolikheten för att en given radonatom skall sönderfalla i ett givet tidsintervall är fix, och oberoende av vad som har hänt tidigare. Bl.a. innebär detta att antalet hål på en given yta av en film är poisson-fördelat med ett väntevärde som är proportionellt mot radonkoncentrationen, exponeringstiden och ytans storlek. Vidare är antalet hål på olika disjunkta (ej överlappande) ytor på en film oberoende stokastiska variabler. Detta är vad som visar sig väsentligt i den fortsatta analysen.

Det datamaterial som vi skall arbeta med har uppmätts genom att ett antal rum i en bostad har försetts med var sin film. Dessa filmer har efter framkallning avlästs på tio olika icke överlappande ytor, med fix storlek, var. Vi inför följande beteckningar:

n = antalet upphängda filmer, dvs antalet rum, γ i = radonkoncentrationen i rum i, mätt i Bq/m ³ ,

X ij = antalet hål i film i på yta j, i = 1, . . . , n, j = 1, . . . , 10.

Enligt ovan gäller då X ij ∈ Po(K γ i ),

där proportionalitetskonstanten K beror på avläsningsytornas storlek och exponeringstiden, men också på bl.a. förstoringen vid avläsningen av filmerna.

4 Arbete med data

Datamaterialet är uppmätt i en nybyggd bostad den 24/3–25/4 1994. Detta skall tolkas så att filmerna hängdes upp vid en viss tidpunkt den första dagen och togs ned vid samma tidpunkt den sista dagen.

i Rum X ij

1 Vardagsrum 20 17 22 15 20 22 24 22 34 20

2 Sovrum 14 15 17 13 14 11 15 16 22 15

3 Mikaels rum 11 17 19 14 25 17 18 16 23 21 Datamaterialet finns i radon200.dat och läses in i Matlab:

>> r200=load(’radon200.dat’)

Kolonn 1 i variabeln r200 innehåller nu mätvärdena för vardagsrummet, kolonn 2 sovrummet och kolonn 3 Mikaels rum.

Konstanten K är 0.0962 för en yta vid 30 dagars exponering. Eftersom den aktuella expone- ringstiden är längre måste en kompensation för detta göras. Enligt resonemanget i förra styc- ket skall detta helt enkelt göras linjärt, eftersom väntevärdena för X -variablerna är proportio- nella mot exponeringstiden. Eftersom våra filmer exponerats 32 dagar bör vårt värde på K vara 0.0962 · 32/30 = 0.1026:

>> K=0.0962*32/30

Syftet med analysen av datamaterialet är att utreda om gränsvärdet på 200 Bq/m ³ överskrids eller

inte. Vi kommer att beräkna punktskattningar av radonkoncentrationen dels för rummen var för

sig, dels för hela huset. Punktskattningarna kommer att kompletteras med motsvarande intervall-

skattningar.

4.1 Punkt- och intervallskattningar

Vi startar med att studera de tre rummen var för sig. Tänk igenom att en väntevärdesriktig punkt- skattning γ ^∗ _i av γ i , i = 1, 2, 3, ges av γ ^∗ _i = 1

10K

10

X

j=1

X ij = X ¯ i

K där ¯ X i är medelvärdet i rum i.

Beräkna skattningarna för datamaterialet ovan:

>> g3rum = mean(r200)/K

Uppgift: Vad blev de tre γ-skattningarna?

För att gå vidare i vår statistiska analys och (så småningom) beräkna konfidensintervall behöver vi ta reda på de statistiska egenskaperna hos punktskattningarna.

Vi har att V(γ ^∗ _i ) = V



 1 10K

10

X

j=1

X ij



 = 1

(10K ) ² ·

10

X

j=1

V (X ij ) = 10K γ i

(10K ) ² = γ i

10K vilket ger medelfelet d (γ ^∗ _i ) för vart och ett av de tre rummen:

>> d3rum = sqrt(g3rum/(10*K)) Uppgift: Vad blev de tre medelfelen?

För att få en uppfattning om hur stor radonkoncentrationen kan tänkas vara i de olika rummen gör vi 95 % konfidensintervall för γ i . Det förutsätter att vi kan normalapproximera, dvs att Y i = P 10

j=1 X ij ∈ Po(10K γ i ) där 10K γ i > 15:

>> g3rum*10*K

Uppgift: Kan vi normalapproximera i alla tre rummen?

Eftersom γ ^∗ _i ∈

∼ N(γ i , p

γ i /(10K )) ges konfidensintervallet av I γ

= γ ^∗ _i ± λ α/ 2 · d(γ ^∗ _i ):

>> I3rum=[g3rum-norminv(1-0.05/2)*d3rum; g3rum+norminv(1-0.05/2)*d3rum]

Den första raden ger de lägre gränserna i de tre intervallen, den andra raden ger de övre gränserna.

Uppgift: Verkar det som om radonkoncentrationen ligger under eller över 200 Bq/m ³ i något av

rummen?

4.2 Hypotesprövning med direktmetoden

Man kan också välja att utföra analysen som ett hypotesprövningsproblem. Den som ska bo i ett rum vill testa

H ₀ :γ i = 200 Bq/m ³ , H ₁ :γ i < 200 Bq/m ³ .

Uppgift: Varför vill invånarna ha ett ensidigt test åt detta hållet?

Tidigare gjorde vi punkt- och intervallskattingar av γ i , vilket ger kvantitativ information om var de sanna värdena kan tänkas ligga. För att göra hypotestest är konfidensmetoden som vi använde i föregående avsnitt användbar. Problemet i det här fallet är dock dels, att det baseras på normal- approximation och testet blir därmed inte exakt, dels, och värre, att signifikansnivån blir fel. Den ba- seras ju på fördelningen när H 0 är sann och vi ska alltså använda γ i = 200 i D 0 (γ ^∗ _i ) = p200/(10K ) istället för d (γ ^∗ _i ). För att komma runt båda problemen kan vi i det här fallet i stället använda di- rektmetoden, dvs vi räknar ut ett P-värde som

P = P(Få det vi fått eller värre om H 0 är sann) och förkastar H 0 om P < α.

För att räkna utan normalapproximation kan vi räkna direkt med observationerna X ij ∈ Po(K γ _i ) och framförallt utnyttja att summan av observationerna i ett rum också är Poissonfördelad, Y i = P 10

j=1 X ij ∈ Po(10K γ _i ).

>> gamma0 = 200;

>> mu03rum = 10*K*gamma0 % väntevärdet för summan när H0 är sann

>> y3rum=sum(r200) % summorna i de tre rummen

>> P3rum=poisscdf(y3rum,mu03rum) % P(Y_i <= y_i) Uppgift: Ska nollhypotesen förkastas i något av rummen?

Summan av samtliga observationer i alla rummen är också poissonfördelad. Vi räknar ut ett P- värde som gäller för hela huset och avgör med direktmetoden om H 0 : γ = 0 skall förkastas, där γ är medelradonkoncentrationen i hela huset:

>> mu0hus = 3*mu03rum % väntevärdet för hela huset när H0 är sann

>> yhus=sum(y3rum) % summan för hela huset

>> Phus=poisscdf(yhus,mu0hus)

Uppgift: Kan vi påstå att medelradonkoncentrationen i huset understiger gränsvärdet?

4.3 Data från äldre hus (för den som har tid)

Nu kan du jobba mer självständigt med ett annat datamaterial av liknande typ. Data är uppmätt i ett äldre hus den 6/12 1993–4/3 1994.

Rum X ij

Sovrum 1 13 9 11 12 10 12 12 14 9 12

Sovrum 2 10 12 8 10 11 10 15 12 12 13

Gillestuga 10 10 15 6 7 10 14 16 12 10

Datamaterialet finns i radon400.dat. Konstanten K är nu 0.00663 (för en yta) vid 30 dagars exponering, dvs. med tanke på den aktuella perioden skall vi räkna med korrigerat K = 0.00663 · 88/30 = 0.0194.

Gränsvärdet för den här typen av bostäder är 400 Bq/m ³ . Om detta överskrids kan fastighetsägaren åläggas att vidta (dyra) åtgärder.

Uppgift: Hur ser analysen ut om ni gör den från ”fastighetsägarens perspektiv”? Hur skiljer det sig från de ”inneboendes perspektiv”?

Uppgift: Skatta γ i för de tre rummen var för sig och beräkna motsvarande konfidensintervall.

Uppgift: Testa H ₀ : γ i = 400 mot H ₁ : γ i > 400 med direktmetoden. Är det något av rummen som behöver åtgärdas?

Uppgift: Testa H ₀ : γ = 400 mot H 1 : γ > 400 för hela huset med direktmetoden. Verkar det bli

dyrt för fastighetsägaren?