• No results found

2Diskretvariabel:Grobarhethossolrosfrön 1Förberedelseuppgifter Datorövning2BetingadfördelningochCentralagränsvärdessatsen

N/A
N/A
Protected

Academic year: 2021

Share "2Diskretvariabel:Grobarhethossolrosfrön 1Förberedelseuppgifter Datorövning2BetingadfördelningochCentralagränsvärdessatsen"

Copied!
8
0
0

Loading.... (view fulltext now)

Full text

(1)

FMS012/MASB03: M

ATEMATISK STATISTIK

, 9

HP

, HT-16

Datorövning 2

Betingad fördelning och Centrala

gränsvärdessatsen

Syftet med den här laborationen är att du skall bli mer förtrogen med några viktiga områden inom kursen nämligen

• Diskreta fördelningar

• Betingade fördelningar

• Centrala gränsvärdessatsen

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 vad en sannolikhetsfunktion är och vad Centrala gränsvärdes- satsen innebär.

3. Redovisas vid laborationens start! Vi köper en påse med 7 solrosfrön. På baksidan i den finstilta texten står det att grobarheten är 75 %. Ange fördelningen för antalet frön som kommer att gro (om de gror oberoende av varandra) samt fördelningens väntevärde och varians.

4. Redovisas vid laborationens start! Vi beräknar medelvärdet ¯X av de oberoende s.v. Xi ∈ Po(3),i = 1, . . . , n. Ange väntevärde och varians för ¯X . Vilken fördelning får ¯X (approxima- tivt) närn är stort? Ungefär hur stort måste n vara för att approximationen ska bli bra?

2 Diskret variabel: Grobarhet hos solrosfrön

Vi vill simulera antalet frön som kommer att gro bland de sju fröna i en fröpåse. Det kan vi göra på två sätt. Det mest rättframma är att simulera 7 frön och räkna antalet som gror. Funktionen rand(1,n)ger en radvektor medn rektangelfördelade slumptal, U , mellan 0 och 1. För att sanno- likheten att ett frö kommer att gro skall blip kan vi helt enkelt se efter om U ≤ p. I så fall kommer fröet att gro. OmU > p så kommer det inte att gro. För att få reda på antalet frön som kommer att gro bland de 7 summerar vi den resulterande 0/1-variabeln:

(2)

>> n = 7;

>> p = 0.75;

>> U = rand(1,n) % 1 rad och n kolumner med observation från R(0,1)

>> U<=p % 0 = gror inte, 1 = gror

>> X = sum(U<=p) % antal frön som gror

Uppgift: Jämför resultatet av U=rand(1,n) och U<=p och förvissa dig om att du förstår vad som hände.

Uppgift: Hur många frön grodde?

Ett smidigare sätt är att utnyttja att vi vet att antalet frön som kommer att gro är Bin(7, 0.75)- fördelat. Då kan vi simuleraX direkt med hjälp av MATLABs färdiga rutiner:

>> help binornd

>> X = binornd(n,p)

Uppgift: Gör om simuleringen några gånger. Hur många frön brukar gro?

Antalet frön som kommer att gro varierar uppenbarligen från gång till gång, dvs från påse till påse.

För att se hur vanligt det är med olika antal frön som kommer att gro simulerar viN = 100 påsar och ritar ett stolpdiagram (vi har ju en diskret variabel).

>> N = 100; % antal fröpåsar

>> X = binornd(n,p,N,1) % X = antal groende frön i var och en av Nx1 påsar

>> antal = hist(X,0:n) % antal = antal av påsarna som ger 0,...,n frön

>> bar(0:n,antal) % stolpdiagram

>> xlabel(’antal frön som gror’)

>> ylabel(’antal påsar’)

Uppgift: Var det någon av påsarna som inte hade några groende frön alls?

Uppgift: Hur många av påsarna gav 5 groende frön?

Uppgift: Hur många av påsarna gav högst 2 groende frön?

Vi vill nu jämföra våra 100 påsar med den teoretiska sannolikhetsfunktionen. För att göra det måste vi skala om y-axeln.

(3)

>> bar(0:n,[antal/N; binopdf(0:n,n,p)]’) % rita två uppsättningar staplar

>> xlabel(’antal frön som gror’)

>> ylabel(’andel påsar’)

>> legend(’simulering’,’exakt’,’Location’,’NorthWest’)

De blå stolparna är vårt simulerade resultat och de röda (gula?) är den teore tiska sannolikhetsfunk- tionen.

Uppgift: Hur stor andel av de simulerade 100 påsarna hade precis 5 groende frön?

Uppgift: Hur stor andel av de simulerade 100 påsarna hade högst 2 groende frön?

Uppgift: Hur stor är sannolikheten att en påse ska ha precis 5 groende frön? (Kolla att det stämmer med binopdf(5,n,p))

Uppgift: Hur stor är sannolikheten att en påse ska ha högst 2 groende frön? (Kolla att det stämmer med binocdf(2,n,p))

Uppgift: Experimentera med att ändra grobarheten från 0.75 till nått annat. Hur ändrar sig för- delningen? Hur ändrar sig sannolikheten att få högst 2 frön som gror?

Uppgift: Ändra också antalet frön i påsarna. Hur ändrar sig fördelningen närn minskar eller ökar?

Hur ändrar sig sannolikheten att få högst 2 frön som gror?

Uppgift: Ändra också pån och p samtidigt.

Om np(1 − p) > 10 kan binomialfördelningen approximeras med en normalfördelning. Vi kan jämföra fördelningsfunktionerna och se hur bra det blir:

>> n = 7;

>> p = 0.75;

>> np = n*p % väntevärde

>> npq = np*(1-p) % varians

>> x = linspace(np-4*sqrt(npq),np+4*sqrt(npq)); % mu +/- 4 sigma

>> stairs(0:n,binocdf(0:n,n,p),’b-’) % Stegfunktion p.g.a diskret s.v.

>> hold on

>> plot(x,normcdf(x,np,sqrt(npq)),’r-.’) % men den här är kontinuerlig

>> hold off

(4)

Uppgift: Pröva med lite olika värden på n och p. Testa både när det går bra att normalapproxi- mera och när det inte går. Beräkna sannolikheten att högst 2 frön gror, både exakt och med normalapproximation (gärna med halvkorrektion).

3 Simulering med hjälp av betingad fördelning: Skördeutfall

Vi tänker oss nu att varje solrosfrö som gror ger upphov till ett Poissonfördelat antal nya frön, i medeltal 50 frön per groende solros. Frön som inte gror ger naturligtvis inga nya frön. Vi är intresserade av fördelningen för det totala antalet nya frön som en fröpåse med 7 frön och 75 % grobarhet kan ge upphov till.

Vi har, som tidigare, X = ”antal frön som gror” ∈ Bin(7, 0.75). Då kommer vi att få att den betingade fördelningen förY = ”antal nya frön”, givet att vi fick X = x frön som grodde, blir Y | X = x ∈ Po(50·x) där x = 0, . . . , 7. Det är alltså summan av x stycken oberoende Po(50)-fördelade variabler, en för varje groende frö. Fördelningen förY ges då av (Satsen om Total Sannolikhet)

pY(y) = P(Y = y) = X

k

P(Y = y | X = x) · P(X = x) =

7

X

x=0

pY |X =x(y) · pX(x)

=

7

X

x=0

e−50x·(50x)y y! ·7

x



· 0.75x · 0.257−x =Nått gräsligt!

För att ta reda på hur denna hiskliga fördelning ser ut börjar vi med att rita upp var och en av de 8 olika möjliga Poissonfördelningarna. Detta är de 8 olika varianterna av betingade fördelning- ar vi har: Po(0), Po(50), Po(100), . . . , Po(350). Vi ritar de 8 sannolikhetsfunktionerna i samma figurfönster men i varsin delfigur för att få lite överblick.

>> close all

>> figure(1)

>> bar(0:7,binopdf(0:7,7,0.75)) % antalet frön som gror

>> xlabel(’antal frön som gror’)

>> figure(2)

>> y = 0:450; % ett lagom intervall för x-axeln

>> for k = 0:7

subplot(4,2,k+1) % rita i delfönster nr 1..8 i en 4x2-plan bar(y,poisspdf(y,50*k))

xlabel(’antal nya frön’) end

Uppgift: Hur ändrar sig fördelningen när antalet groende frön ändrar sig?

Uppgift: Tänk efter hur fördelningen för Y borde se ut, när vi viktat ihop dessa 8 fördelningar med vikter enligt binomialfördelningen för antalet groende frön.

(5)

Även om vi inte vill räkna ut sannolikhetsfunktionen förY så är det ganska enkelt att låta Matlab göra det:

>> pY=zeros(size(y)); % Fyll först p_Y(y) med nollor.

>> for k=0:7 % Uppdatera p_Y(y) för varje k pY=pY+poisspdf(y,50*k)*binopdf(k,7,0.75);

end

>> figure(3)

>> bar(y,pY)

>> xlabel(’antal nya frön’)

Uppgift: Ser fördelningen ut som du hade väntat dig?

Den specialskriva funktionen solrosor.m, som finns på kurshemsidan, ritar upp sannolikhets- funktionen för Y där Y | X = x ∈ Po(μ · k) och X ∈ Bin(n, p) för valfria värden på n, p och μ.

Spara ner funktionen på din dator och använd den för att rita om figuren:

>> help solrosor

>> solrosor(7, 0.75, 50)

Uppgift: Experimentera med olika värden på n, p och μ. Vad händer om antalet frön i påsen, n, minskar eller ökar? Om grobarheten,p, minskar eller ökar? Om medelantalet nya frön per frö som gror, μ, minskar eller ökar?

Uppgift: Vad händer om grobarheten är 100 %?

Uppgift: Kan du få fördelningen att se normalfördelad ut?

4 Centrala gränsvärdessatsen

Vi skall nu titta lite närmare på Centrala Gränsvärdessatsen (CGS). Vi börjar med en liten simule- ring från en känd fördelning, två slumpmässiga obervationerx1,x2 frånX ∈ Po(μ) där μ = 3. Vi ska sedan beräkna medelvärdet ¯x och se hur nära väntevärdet μ det hamnar.

>> my = 3; % det sanna my-värdet

>> x = poissrnd(my,2,1) % en 2x1-matris med Po(my)-slumptal

>> xmedel = mean(x) % medelvärdet

Uppgift: Gör om simuleringen och medelvärdesberäkningen några gånger. Verkar medelvärdet variera mindre än de enskilda observationerna? Borde den det? I så fall, hur mycket mindre?

(6)

Låt oss göra om simuleringarna ett stort antal gånger så vi får bättre uppfattning om hur medelvärdet beter sig:

>> my = 3;

>> n = 2; % antal termer i medelvärdet

>> M = 1000; % antal simuleringar

>> x = poissrnd(my,n,M) % n x M-matris. x1 i första raden, xn i sista.

>> xmedel = mean(x) % M st medelvärden

>> subplot(2,1,1)

>> hist(x(1,:),0:15) % histogram över de Mst x1-värdena

>> subplot(2,1,2)

>> hist(xmedel,0:0.5:15) % histogram över de Mst x-medelvärdena

Uppgift: Experimentera med lite olika värden pån och se vad som händer med medelvärdet. Du kan behöva ändra klassbredden i det undre histogrammet för att se något, t.ex. 0:0.1:15, som ger klasser från 0 till 15 med bredd 0.1.

Uppgift: Jämför variationen hos de enskilda observationerna i övre figuren och variationen för skattningarna i undre figuren. Ändrar sig variationen hos observationerna när vi ändrarn?

Det är uppenbart att medelvärdet ¯X varierar när vi får nya observationer. Det varierar också mindre när den baseras på fler observationer. Dess fördelning ska också bli mer och mer lik en normalför- delning ju fler observationer vi har. Det ska vi undersöka närmare nu.

Eftersom ¯Xn = 1 n

n

X

i=1

Xi där

n

X

i=1

Xi ∈ Po(nμ) ≈ N(nμ,

nμ) om nμ > 15, så gäller också att X¯nN(μ, p

μ/n) om nμ > 15, dvs om n > 15/μ.

Vi hade tidigare μ = 3, dvs n = 5 var tillräckligt stort för att medelvärdet skulle bli ungefär normalfördelat. Vi simulerar nu från en Poissonfördelning med ett mindre μ så att det är lättare att se hur ”onormalfördelat” det blir för smån.

För att ni ska slippa skriva så mycket kod finns den specialskrivna Matlab-funktionen poisson_cgs.m att ladda ner från kurshemsidan.

>> help poisson_cgs

>> poisson_cgs

Uppgift: Om ¯Xn är ungefär normalfördelad ska den oftast (i 95 % av fallen) hålla sig inom de trumpetformade gränserna. Vad händer när n är för litet? Vad händer när n växer? Jämför med de utritade fördelningarna för de tren-värdena.

Uppgift: För smån kan normalfördelningen uppenbarligen bli negativ. Blir ¯Xnnågonsin negativ?

(7)

Uppgift: Hur är det med den övre gränsen? Hamnar ¯Xnöver den så sällan som den borde närn är litet? Hur blir det närn växer?

Experimentera vidare med olika värden på μ, t.ex.:

>> poisson_cgs(’my’,3) % my=3

>> poisson_cgs(’my’,3,’n’,10,’logx’,false) % my=3, n=1..10, ej logskala

>> poisson_cgs(’my’,0.01,’n’,2000) % my=0.01, n=1..2000

4.1 Allmän diskret fördelning (om du hinner)

Vi ska nu experimentera med en diskret fördelning där vi inte på förhand vet när n är tillräckligt stort. Vi börjar med den likformiga fördelningen över 1 t.o.m. 6, dvs ett tärningskast. Mata in denna sannolikhetsfunktion i form av en vektor och rita upp den:

>> pX=[1 1 1 1 1 1]/6

>> bar(1:6,pX)

Uppgift: Tänk efter vilka värden summan av två tärningskast kan anta och hur fördelningen bör se ut.

Den specialskrivna Matlab-funktionen p_summaplot.m (finns på kurshemsidan) beräknar och ri- tar upp fördelningen för en summa avn observationer från fördelningen pX(k) tillsammans med en normalfördelning med samma väntevärde och standaradavvikelse som summan.

>> help p_summaplot

>> p_summaplot(pX,1) % samma fördelning som innan.

Uppgift: Fundera på hur man har räknat ut väntevärdet och standardavvikelsen!

Nu ritar vi upp fördelningen för summan av tvåX från denna fördelning:

>> p_summaplot(pX,2)

Uppgift: Stämmer det med det du trodde? Ser det normalfördelat ut?

Uppgift: Vad hände med väntevärdet och standardavvikelsen?

(8)

Det behövs tydligen lite fler termer i summan innan det blir normalfördelat.

Uppgift: Ökan tills fördelningen börjar se normalfördelad ut.

Uppgift: Hur ändrar sig väntevärdet och standardavvikelsen när du ökarn?

Eftersom fördelningen var symmetrisk och ”snäll” blev den ganska snabbt normalfördelad. Använd en sned fördelning:

>> pX = [1 1 2 8]/12

>> p_summaplot(pX,1) ...

Uppgift: Hur många komponenter behövs det nu i summan för att fördelningen ska kunna ap- proximeras med en normalfördelning?

Uppgift: Hitta på fler fördelningar och experimentera. Några varianter att testa:

>> pX=[3 2 1 0 1 2 3]/12 % U-formad fördelning

>> pX=[10 5 2 0 0 0 0 1]/18 % Skev med ett extremt värde

References

Related documents

En större andel patienter som fått stomi på grund av rektalcancer, jämfört med de patienter med stomi som inte fått denna sjukdomsdiagnos, uppgav att de var missnöjda med både

Detta för att ge en bakgrund till hur socialtjänsten arbetar med våld i nära relation idag och därmed även till socialtjänstens förutsättningar för att arbeta med

Det som kan uppstå är att de hamnar i samma rockring då kör jag med att en av eleverna får backa om de inte kommer överens om vem som ska backa så får de köra sten sax och påse

Vissa av deltagarna upplevde detta olika beroende på situationen och uttryckte att det blev mer pinsamt i offentligheten där personerna runt omkring inte alltid hade någon kunskap

Förekomsten av mycket hygroskopiska föreningar i aerosoler kan påskynda processen för bildandet molndroppar, medan närvaron av mindre hygroskopiska ämnen kan förlänga den tid som

I Kardborreprojektet, som genomförts i samarbete mellan Bioresurs och Nationellt centrum för matematikutbildning (NCM), finns för- slag på hur man kan använda frön

När fröet tar upp vatten, ökar enzymaktiviteten och upp- lagsnäringen bryts ner till enklare föreningar som transporteras till fröets tillväxtpunkter. Först börjar roten

Förändringar av kottars form vid olika fuktighet gör att djur och andra figurer som tillverkas av kottar också kommer att ändra form.. Från kotte