2Matlab–deförstastegen 1Förberedelseuppgifter Datorövning1:Fördelningar

FMS012/MASB03: M

ATEMATISK STATISTIK

, 9

HP

, HT-16

Datorövning 1: Fördelningar

I denna datorövning ska du utforska begreppen sannolikhet och fördelningar genom numeriska exempel i Matlab.

Du behöver en Matlab-installation som inkluderar Statistics Toolbox. De extra filer du behöver finns att ladda ner från kursens hemsida www.maths.lu.se/kurshemsidor/fms012masb03

1 Förberedelseuppgifter

1. Läs igenom denna handledning.

2. Förvissa dig om att du förstår vad täthetsfunktion och fördelningsfunktion är och hur de förhåller sig till varandra.

3. Redovisas vid laborationens start! Skriv upp täthetsfunktionen förX ∈ N (μ, σ)-fördelad s.v. och skissa upp den. Ange väntevärde och standardavvikelse förX .

4. Redovisas vid laborationens start! OmX har en ”standardnormalfördelning”, vad är då μ och σ?

2 Matlab – de första stegen

(Hoppa över detta om du redan är van vid Matlab) Matlab¹ tillåter användaren att kombinera numeriska beräkningar med avancerad grafik. Kortare kommandon kan köras interaktivt, men för mer komplicerade problem är det möjligt att skriva program och definiera egna funktioner.

I tillägg till Matlab finns flera s.k. Toolboxes (verktygslådor) för olika tillämpningar, t.ex. signalbe- handling, reglerteori och finita elementmetoder.

Under datorövningarna kommer vi att använda, bland annat, Statistics Toolbox.

Matriser och vektorer i Matlab

Matlab kan användas som en avancerad miniräknare: de vanligaste funktionerna är fördefinierade.

Vid Matlab-prompten (>>), kan du t.ex. beräkna√

1.19²− 1 + sin(π/2) + e²genom att skriva

>> sqrt(1.19^2-1)+sin(pi/2)+exp(2) och resultatet dyker upp.

1Mer information om Matlab finns på http://www.mathworks.com/

Uppgift: Beräkna ovanstående uttryck i Matlab.

När du vill veta mer om de fördefinierade funktionerna i Matlab är help-kommandot använd- bart. Det är en god idé att använda det under övningarna, även om det inte står uttryckligen i handledningen. Börja med

>> help help

Sedan kan du, t.ex. skriva

>> help log

för att få reda på vilken bas Matlab använder för logaritmfunktionen.

Matlab är en förkortning av Matrix laboratory, och matriser och vektorer är karakteristiskt för Matlab. Alla data är lagrade i vektorer och matriser. (Med vektor menar vi en rad- eller kolonn- matris.) Matrisen

A = 2 0

3 1



skrivs in i Matlab på följande sätt:

>> A = [2 0; 3 1]

En radvektor kan t.ex. skrivas på följande sätt:

>> v=[0 0.1 0.2 0.3]

och en kolonnvektor:

>> vv=[0; 0.1; 0.2; 0.3]

Kommandot length (eller size) ger storleken på vektorn eller matrisen:

>> vLength=length(v)

>> vvLength=length(vv)

>> ASize=size(A)

Vi kan plocka ut enskilda element ur vektorer och matriser på följande vis. Säg att vi vill komma åt värdet av det fjärde elementet i vektorn v och dessutom värdena på de tre första elementen. Det gör vi på följande vis:

>> v(4), v(1:3) eller, tillsammans:

>> v([4 1:3])

Elementen i en vektor kan sorteras i stigande ordning:

>> u=[8 -3 2.7]

>> uSorted=sort(u)

Att hantera variabler och data

Vi har nu definierat ett antal variabler och en lista på de variabler som finns i Matlabs minne fås med kommandot who. Kommandot whos ger samma lista men utökad med storleken på variablerna.

Uppgift: Kör båda kommandona. Känner du igen variablerna i listan?

Vi avslutar sessionen med att ta bort variablerna.Alla variablerna tas bort med kommandot clear.

Använd help clear för att ta reda på hur du tar bort bara vissa variabler.

Grafik

I den här delen ska du göra några enkla figurer i Matlab. Efteråt kommer du att kunna rita en figur över en funktionx 7→ f (x). Som exempel, låt oss välja f (x) = sin x för 0 < x < 4π. Skapa först variablerna x respektive y:

>> x=0:0.05:4*pi % Detta efter procenttecknet är en kommentar, skriv

% inte av den.

>> y=sin(x) % x=0, 0.05, 0.1, 0.15, ..., 4pi Två vektorer av samma längd kan ritas mot varandra så här:

>> plot(x,y)

Ett grafikfönster dyker nu upp (om det inte redan fanns ett) som figuren ritas i. Figuren heter figure(1). Man kan ha flera grafikfönster. Vill du att nästa figur ska komma i ett nytt fönster, istället för att rita över den första figuren, så kan du skapa ett nytt grafikfönster med kommandot figure(2).

Man kan ge flera optioner till plot-kommandot, t.ex., färg:

>> plot(x,y,’r’) % Röda streck

Vi kan också välja att rita ut de enskilda punkterna som stjärnor istället för att dra linjer mellan dem:

>> plot(x,y,’*’)

Man kan också kombinera optionerna. Se help plot för att ta reda på vad följande kommando borde göra. Kolla sedan att det blev så också:

>> plot(x,y,’md-’)

Man kan använda kommandot axis för att titta på en bestämd det av figuren. Pröva med

>> axis([0 10 -1.5 1.5]) Vad hände?

Det är ofta lättare at tolka en figur om man lägger till en grid: ta reda på hur kommandot grid används och lägg till en grid till din figur.

Det aktuella figuren töms med kommandot clf. Det tomma fönstret blir kvar. Vill du ta bort det också ska du använda close istället.

3 Relativa frekvenser och fördelningar

(Här börjar den egentliga laborationen)

I denna del ska vi använda numeriska exempel i Matlab för att studera koncepten ” sannolikhet”

och ”fördelning”. Målet är att du ska få en intuitiv känsla för sannolikhetsresonemang, snarare än att konfronteras med teori.

Data-undersökning

För att illustrera syftet använder vi artificiella data som är simulerade från en statistisk fördelning.

Detta i motsats till verkliga data där det inte finns några etiketter som säger vilken fördelning det är. Trots att vi vet hur data genererades är det ändå användbart och man använder ofta simulerade data i skattningar och test i mer komplicerade situationer.

För att skaffa dig ett slumpmässigt dataset med 50 värden, skriv

>> data=randn(1,50)

Uppgift: Vilken fördelning kommer ditt slumpmässiga stickprov från (använd help randn)? Vil- ka värden har parametrarna i den? Skriv ner täthetsfunktionen.

En god regel, när man står inför ett nytt datamaterial, är att rita upp det på några olika sätt.

Vi börjar med att göra ett histogram:

>> hist(data)

Uppgift: Ser det ut som du väntade dig? Jämför med täthetsfunktionen.

Använd nu kommandot

>> figure(2) % Ritar i ett nytt fönster

>> plot(data,’.’)

och relatera det till histogrammet.

Uppgift: Jämför histogrammet med ploten. Hur syns egenskaperna hos data i histogrammet, och tvärtom?

Ett annat sätt är att rita de sorterade data, med ordningsnumret på y-axeln:

>> figure(3)

>> plot(sort(data),1:length(data),’.’)

Uppgift: Jämför denna plot med figure(1) och figure(2). Hur hänger de ihop med varandra?

Uppgift: Välj ut några datapunkter i figure(2) och försök hitta dem i de andra två figurerna.

I figure(3) kan vi t.ex. avläsa hur många av observationerna som är mindre än eller lika med ett visst tal.

Uppgift: Väljx = 1.1 och försök avgöra i figuren (det går att zooma) hur många av värdena som är mindre än eller lika med 1.1.

När antalet observationer i stickprovet ökar kan vi tolka kvoten som sannolikheten att få ett värde mindre än eller lika medx. Kvoten kan beräknas så här:

>> ratio = sum(data<=1.1) / length(data) Uppgift: Stämmer det med din uppskattning från figuren?

För att förstå hur data<=1.1 fungerar så jämför vi det med ursprungsdata:

>> data

>> data<=1.1

Vad är det som händer?

Uppgift: Pröva med några andra värden påx. Hur borde andelen ändra sig när x ökar respektive minskar? Jämför med figuren.

Den omvända proceduren, hitta det värde x som motsvarar en given sannolikhet, dvs en given kvantil, är ofta viktigare. Vi återkommer till det lite senare.

Vi kan naturligtvis låta datorn välja ett stort antal värden att undersöka och sedan försöka få en överblick. Eftersom vi har ett ändligt antal observationer så blir antalet, eller andelen, observationer som än mindre än eller lika med ett visstx-värde en stegfunktion som vi kan rita upp:

>> figure(4)

>> stairs(sort(data),(1:length(data))/length(data),’-’)

>> grid on

−3 −2 −1 0 1 2 3 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figur 1: Empirisk fördelningsfunktion, ett exempel

Figuren bör likna Figur 1 i handledningen och din egen figure(3), bortsett från y-skalan. Den visar hur värdena är fördelade och denna typ av figur kallasempirisk fördelningsfunktion (empirical distribution function². För ett värde på x-axeln, t.ex. 1.1, hittar vi, på y-axeln, andelen värden som är mindre än eller lika värdet på x-axeln.

Uppgift: Kolla att andelen värden som är mindre än eller lika med 1.1 stämmer med det du fick fram tidigare.

Större stickprov. Fördelningsfunktionen för en slumpvariabel

Låt oss nu studera ett större datamaterial, t.ex. 2000 observationer från samma fördelning som tidigare. Vi simulerar data och ritar dem i en ny figur:

>> data=randn(1,2000);

>> figure(5)

>> hist(data)

>> figure(6)

>> stairs(sort(data),(1:length(data))/length(data),’.-’)

>> grid on

Uppgift: Jämför histogrammet med det i figure(1). Hur förändrades det när du fick fler obser- vationer?

Uppgift: Jämför den empiriska fördelningsfunktionen med den i figure(4). Hur förändrades den?

2Fördelningsfunktioner kallas oftacumulative distribution functions.

Uppgift: Vad blir nu andelen värden som är mindre än eller lika med 1.1?

Med många observationer närmar sig resultatetfördelningsfunktionen, dvs, för en slumpvariabel X , funktionen

FX(x) = P(X ≤ x).

I vårt fall valdes X från en normalfördelning; vi hade X ∈ N(0, 1). Vi ritar in den teoretiska fördelningsfunktionen, normcdf, i samma figur som de två empiriska:

>> x=linspace(-4,4,500); % 500 värden jämnt fördelade mellan -4 och 4

>> figure(4)

>> hold on % Fortsätt rita fler saker i samma figur.

>> plot(x,normcdf(x),’r’)

>> hold off % Sluta rita i samma figur.

>> figure(6)

>> hold on % Fortsätt rita fler saker i samma figur.

>> plot(x,normcdf(x),’r’)

>> hold off % Sluta rita i samma figur.

För alla fördelningsfunktioner FX, har vi att FX(x) → 1 när x → ∞ och att FX(x) → 0 när x → −∞.

Uppgift: Tolka figuren. Vad är det på x- och y-axlarna?

Uppgift: Jämför hur bra de empiriska fördelningsfunktionerna följer den teoretiska i de två figu- rerna. Vad hände när antalet observationer ökade?

Uppgift: Läs av P(X ≤ 1.1) ur den teoretiska fördelningsfunktionen i figuren och jämför med dina tidigare skattningar. Jämför också med det exakta värdet som kan fås med normcdf(1.1).

Kvantiler

Begreppetkvantil är viktigt. Kvantilen kan definieras på olika sätt men vi (och många andra) an- vänder följande definition: kvantilen är det talxαsom uppfyller

P(X ≤ xα) = 1 − α (1)

där α är ett tal mellan 0 och 1 (vanliga val är: 0.05, 0.01, 0.001).

Uppgift: Läs av kvantilenx0.05 där α = 0.05 ur dina figurer, med hjälp av definitionen (1). Både som skattningar i de två empiriska fördelningsfunktionerna och exakt i den teoretiska. Jämför med det exakta värdet, som kan fås med norminv(1-0.05).

Uppgift: Experimentera med att ändra antalet observationer. Simulera nya slumptal, rita nya histo- gram och empiriska fördelningsfunktioner, samt skatta P(X ≤ 1.1) och x0.05.

Uppgift: Använd ett mycket litet dataset, t.ex. 5 observationer och gör om simuleringarna och skattningarna några gånger. Verkar de tillförlitliga?

Uppgift: Använd ett större dataset, t.ex. 100 observationer och gör om simuleringarna och skatt- ningarna några gånger. Verkar de mer tillförlitliga nu?

Hur datasetets storlek påverkar osäkerheten i uppskattningarna kommer vi tillbaka till under hela resten av kursen.

Andra fördelningar

Vi ska nu rita upp några normalfördelningar, N (μ, σ), och se hur de ändrar sig när vi ändrar på parametrarna μ och σ.

>> close all % stäng alla gamla figurer

>> x = linspace(0,10,1000); % Genererar 1000 tal jämnt utspridda

% mellan 0 och 10.

>> figure(1)

>> plot(x,normpdf(x,2,0.5)) % N(2, 0.5)

>> hold on % Lås plotten, övriga ritas i samma bild.

>> plot(x,normpdf(x,7,0.5),’r’) % N(7, 0.5) i rött

>> plot(x,normpdf(x,5,2),’g’) % N(5, 2) i grönt

>> plot(x,normpdf(x,5,0.2),’y’) % N(5, 0.2) i gult

>> hold off % Lås upp plotten

>> xlabel(’x’)

>> title(’Täthetsfunktioner, f(x)’)

>> figure(2)

>> plot(x,normcdf(x,2,0.5))

>> hold on

>> plot(x,normcdf(x,7,0.5),’r’)

>> plot(x,normcdf(x,5,2),’g’)

>> plot(x,normcdf(x,5,0.2),’y’)

>> hold off

>> xlabel(’x’)

>> title(’Fördelningsfunktioner, F(x)’)

Uppgift: Vad händer med fördelningen när μ ändras? Vad representerar μ i fördelningen?

Uppgift: Vad händer med fördelningen när σ ändras? Vad respresenterar σ i fördelningen?

Uppgift: Fördelningsfunktionen är ju integralen av täthetsfunktionen. Relatera dem till varandra i figuren. Hur ändrar sig, t.ex. fördelningsfunktionen närx ligger nära μ jämfört med när x ligger långt från μ? Hur ser täthetsfuktionen ut då (stor eller liten?)

Uppgift: Experimentera med andra värden på μ och σ och se vad som händer. Du kan behöva ändra x för att för att få plats i figuren (tips: det allra mesta av en normalfördelning ryms inom μ ± 4σ).

Jfr. Uppgift 6.7 : Elförbrukningen (kWh) vid en kemisk tillverkningsprocess varierar från dag till dag som en s.v.X ∈ N (180, 5).

Uppgift: Rita upp fördelningsfunktionen förX och avläs sannolikheten att elförbrukningen en viss dag är minst 170 kWh. Jämför med det exakta värdet 1-normcdf(170,180,5).

Uppgift: Utnyttja figuren för att bestämma P(170 ≤ X ≤ 195). Jämför med exakta värdet normcdf(195,180,5)-normcdf(170,180,5).

Uppgift: Läs av 1 %-kvantilen för elförbrukningen i figuren.

Jämför med det exakta värdet norminv(1-0.01,180,5).