• No results found

Laboration 1: Grundläggande sannolikhetsteori, simulering och dataanalys

N/A
N/A
Protected

Academic year: 2022

Share "Laboration 1: Grundläggande sannolikhetsteori, simulering och dataanalys"

Copied!
9
0
0

Loading.... (view fulltext now)

Full text

(1)

Laboration 1

Matematisk statistik AK för Π och E, FMS012, HT14/VT15

Laboration 1: Grundläggande sannolikhetsteori, simulering och dataanalys

Syftet med den här laborationen är att du skall bli mer förtrogen med några grundläg- gande områden inom matematisk statistik nämligen

• Betingade sannolikheter

• Stokastiska variabler och deras fördelningar.

• Simulering

• Dataanalys

Vidare skall du lära dig lite om hur man kan använda Statistics toolbox i matlab.

1 Förberedelseuppgifter

Som förberedelse till laborationen bör du läsa igenom Kapitel 3 och 8 samt hela labo- rationshandledningen. Repetera dessutom det som sades om normalfördelningsplot på föreläsning 3.

Till laborationens start har du med dig lösningar till uppgifterna a) – g), dessa kan komma att samlas in av handledaren. Observera att godkända uppgifter är ett krav för att bli godkänd på laborationen.

a) Definiera den betingade sannolikheten för A givet B. (Behövs till avsnitt 5.2).

b) Definiera följande funktioner: fördelningsfunktion, täthetsfunktion och sannolik- hetsfunktion. Redogör för sambanden mellan dem. (Behövs till avsnitt 4).

c) Ange täthetsfunktionerna för normalfördelning, Rayleighfördelning. Hur ser san- nolikhetsfunktionen ut för en binomialfördelning? (Behövs till avsnitt 3).

d) Definiera α-kvantilen för en kontinuerlig s.v. X. (Behövs till avsnitt 3.5).

e) Den stokastiska variabeln X har täthetsfunktionen fX(x) = x/2 för 0 ≤ x ≤ 2.

Bestäm P(1 ≤ X ≤ 2). (Behövs till avsnitt 3).

f) Låt X vara en diskret stokastisk variabel med sannolikhetsfunktion pX(k) = 0.6 · 0.4k, k = 0, 1, 2 . . . .

Bestäm P(1 ≤ X < 4). (Uppvärmning till avsnitt 3).

g) Antag att du har tillgång till rektangelfördelade slumptal på intervallet (0, 1). Hur kan du med hjälp av inversmetoden generera slumptal från en exponentialfördelning med parameter λ? (Behövs till avsnitt 4).

(2)

2 Introduktion till Matlab (kan hoppas över om ni använt Matlab innan

I alla laborationer i den här kursen kommer matlab, som är ett interaktivt program för numeriska beräkningar, att användas. Matlab finns på de flesta datorer på LTH, och till fördelarna med programmet hör att det ser i stort sett likadant ut oberoende av på vilken sorts dator man kör det. Om du vill ha mer att läsa om Matlab så finns det olika handledningar, tex:

http://www.maths.lth.se/matstat/kurser/fms012/cdi/lab/fms012_lab0_13.pdf Börja med att logga in på ditt vanliga konto och starta Matlab genom att helt enkelt skriva matlab. Programmet avlutas med kommandot exit. Datamaterial och specialruti- ner finns att hämta på kursens hemsida:

http://www.maths.lth.se/matstat/kurser/fms012/pie/

Till att börja med kan man tänka på Matlab som en avancerad räknedosa som beräknar uttryck. Man skriver in vad man vill ha gjort och Matlab svarar.

>> 3*11.5+2.3^2/4 ans =

35.8225

Variabler tilldelas värden med tecknet = och finns sedan kvar i minnet. Pröva att tilldela några variabler värden.

>> a=1;

>> b=sqrt(36);

>> width=3.89;

>> who

Your variables are:

a b width

Kommandot who visar alltså vilka variabler som finns i minnet. En variabel, t ex b, kan raderas ur minnet med clear b. Läggs ett semikolon, ;, till efter en kommandorad, skrivs resultatet inte ut på skärmen. Matlab kan också hantera vektorer och matriser, och de hanteras precis lika enkelt som ovan.

>> x=[1 3 7]

x =

1 3 7

>> y=[2 1 8]’

y = 2 1 8

>> z=[1 2 ; 3 4]

z =

1 2

3 4

>> w=rand(1,4) w =

0.2190 0.0470 0.6789 0.6793

(3)

Tecknet ’ betyder som synes transponat och semikolon används för att skilja rader åt i matriser. Funktionen rand(m,n) ger en m × n-matris med slumptal som är rektangel- fördelade mellan 0 och 1. Kommandot whos visar vilka variabler som finns i minnet, och deras storlek.

>> whos

Name Size Bytes Class Attributes

a 1x1 8 double

w 1x4 32 double

width 1x1 8 double

x 1x3 24 double

y 3x1 24 double

z 2x2 32 double

Notationen för algebraiska uttryck är den vanliga, men kom ihåg att multiplikationsteck- net * tolkas som matrismultiplikation. Elementvis multiplikation mellan två matriser A och B av samma dimensioner skrivs A.*B.

I Matlab finns alla vanliga funktioner inbyggda, t ex exp log log10 sin asin cos acos tan atan.

Observera att log är den naturliga logaritmen. Pröva att plotta en funktion, t ex genom följande kommandon.

>> x=0.5:0.1:2

>> help log

>> y=log(x)

>> plot(x,y)

Precis som i C och andra programmeringsspråk kan man i matlabskriva egna funktioner, vilket kan vara mycket praktiskt om man skall utföra samma typ av beräkningar flera gånger. Hur man gör är i grunden rätt enkelt. Man skriver i stort sett bara de vanliga matlab kommandona i en särskild s.k. m-fil. Låt oss till exempel skriva en funktion som plottar sinus- och cosinus- och sparar värdena i en matris. Börja med att starta din favoriteditor, och börja skriv programkoden

function T=plotta(a,b)

% T=plotta(a,b)

% Plottar sinus och cosinus i intervallet [a,b] samt

% lagrar funktionsvärdena i matrisen T.

x=linspace(a, b, 100);

y=sin(x);

z=cos(x);

plot(x,y,’y’) hold on

plot(x,z,’g’) hold off T=[y’ z’];

Spara sedan filen under filnamnet plotta.m (suffixet .m är nödvändigt för att matlabskall kunna hitta filen). Nu kan du använda funktionen plotta på samma sätt som alla andra funktioner (om du t.ex. ger kommandot help plotta, så kommer de kommentarer du själv skrivit efter %-tecken fram).

(4)

3 Introduktion till Statistics toolbox

Här skall vi titta lite på vad man kan använda Statistics toolbox i matlabtill. Statistics toolbox är ett programpaket till matlabsom innehåller funktioner som är användbara inom statistikområdet. Om du skriver

>> help stats

så skrivs en lista ut på skärmen innehållande alla de funktioner som programpaketet innehåller. I den här laborationen kommer vi främst att vara intresserade av funktionerna ...pdf, ...cdf, ...inv och ...rnd.

För några fördelningarna har vi alltså följande funktioner

Fördelning fX(x) eller px(k) FX(x) F−1X (x) slumptal

Normalförd. X ∈ N(μ, σ) normpdf(x,μ,σ) normcdf(x,μ,σ) norminv(x,μ,σ) normrnd(μ,σ,r,c) Exponentialförd. X ∈ Exp(λ) exppdf(x,1/λ) expcdf(x,1/λ,σ) expinv(x,1/λ) exprnd(1/λ,r,c) Binomialförd. X ∈ Bin(n, p) binopdf(k,n,p) binocdf(x,n,p) binoinv(x,n,p) binornd(n,p,r,c) Om man sedan är intresserad av att se hur någon av funktionerna fungerar skriver du

bara

>> help funktionsnamn

T.ex., du är intresserad av vad funktionen normcdf gör. Skriv

>> help normcdf

Här nedan följer ett antal enkla uppgifter som ni skall lösa med hjälp av matlab. Kom ihåg att använda funktionen help flitigt!

1. Beräkna P(X ≤ 5) då X är a) Normalfördelad med μ = 2, σ = 3, b) Exponential- fördelad med λ = 1/3

Svar: . . .

2. Beräkna P(X > 57) då X ∈ Bin(200, 0.3).

Svar: . . .

3. Beräkna P(50 < X ≤ 90) då a) X ∈ N(100, 40), b) X ∈ Bin(200, 0.3).

Svar: . . .

4. Plotta täthetsfunktionen och fördelningsfunktionen för X ∈ N(3, 2).

5. Använd kommandot norminv för att finna 5%, 10%, 25%, 50%, 75% och 99%

kvantilerna för X ∈ N(50, 10).

Svar: . . .

(5)

4 Simulering

Det vanligaste sättet att simulera stokastiska variabler med en given fördelningsfunktion är att utgå från slumptal som är rektangelfördelade på (0, 1) och sedan transformera dem på lämpligt sätt. Vi skall i det följande koncentrera oss på simulering av kontinuerliga stokastiska variabler.

Antag att vi vill simulera variabler med fördelningsfunktionen F(x) och att F(x) är in- verterbar. Kalla inversen F−1(x) och låt U vara en Rekt(0,1)-fördelad stokastisk variabel.

Om vi sätter X = F−1(U) visar räkningen

P{X ≤ x} = P{F−1(U) ≤ x} = P{F(F−1(U)) ≤ F(x)} = P{U ≤ F(x)} = F(x) att X har fördelningsfunktionen F. Observera att det är väsentligt att F(x) är växan- de för att räkningen ovan skall vara korrekt, och det är ju också fallet ty F(x) är en fördelningsfunktion.

Nu skall vi generera exponentialfördelade slumptal med parametern λ, X ∈ Exp(λ). Hur gör vi detta om vi utgår från U ∈ R(0, 1)?

Svar: X = . . .

Skapa sedan 1000 exponentialfördelade slumptal med parametern λ = 1/3 med hjälp av Matlabs funktion rand och lagra dem i vektorn x. Det kan vara bra att veta att om man ger funktioner som exp, log etc en vektor som inparameter så returnerar Matlab en vektor av samma storlek med funktionen utförd elementvis.

Ett histogram över vektorn x med 30 klasser kan erhållas med

>> help hist

>> mid=.5:1:30;

>> hist(x,mid)

Vi vill dock också kunna jämföra histogrammet med den teoretiska täthetsfunktionen, och då blir det lite krångligare eftersom histogrammet måste skalas om för att vara jämförbart med en täthetsfunktion vars area under funktionen är 1.

>> [n,mid]=hist(x,mid);

>> width=1;

>> bar(mid,n/1000/width)

>> hold on

>> xx=linspace(0,30)

>> plot(xx,1/3*exp(-xx/3))

>> hold off

Det första kommandot beräknar histogrammet, men det ritas inte upp utan istället lagras mittpunkterna för klasserna i vektorn mid och antalet observationer i respektive klass lagras i n.

Med komandot bar ritas histogrammet med klassbredden 1 upp. Den här gången norme- rar vi antalet observationer dels med det totala antalet, vilket ger den relativa andelen observationer i klassen, dels med klassbredden, vilket ger en uppskattning av täthets- funktionens värde i mittpunkten, ty om Δ är klassbredden så gäller

fX(x)Δ ≈

Z x+Δ/2 x−Δ/2

fX(x) dx = FX(x+Δ/2)−FX(x−Δ/2) = P(x−Δ/2 < X < x+Δ/2).

(6)

Slutligen ritas den exakta tätheten upp. Ser det ut som om X är exponentialfördelad?

Svar: . . .

På samma sätt som att histogrammet kan ses som en uppskattning av täthetsfunktionen för X så kan den s.k. empiriska fördelningsfunktionen

FX(x) = {antal xi≤ x; i = 1, 2, . . . , N}

N

ses som en uppskattning av fördelningsfunktionen för X, FX(x) = P(X ≤ x).

Nu använder vi funktionen stairs för att jämföra den empiriska fördelningsfunktionen med den teoretiska

>> stairs(sort(x),(1:1000)/1000)

>> hold on

>> plot(xx,1-exp(-xx/3))

>> hold off

Ser det ut som om X är exponentialfördelad?

Svar: . . .

5 Dataanalys - kvalitetskontroll på serieproducerade mobiltelefoner

En GSM mobiltelefon måste uppfylla vissa krav (t.ex. mottagarkänslighet och störtålig- het) som telebolagen har ställt upp. Därför måste alla serieproducerade telefoner testas för att se om de uppfyller kravspecifikationen. Det kan ju vara så att p.g.a. felakti- ga/dåliga komponenter eller dåliga lödningar så klarar en serieproducerad mobiltelefon inte de givna kraven och måste kasseras, eller rättas till manuellt.

I filerna phase_error.mat och sensitivity.mat finns uppmätta värden på två av de krav som ställs på radiomottagaren på ett antal serieproducerade mobiltelefoner från ett icke namngivet stort svenskt telekommunikationsföretag. Mätningarna gjordes under ett antal dagar i oktober 97.

5.1 Beskrivning av fasfel- och känslighetskrav

I phase_error.mat har man mätt upp fasfelet i radiomottagaren i telefonen. GSM har ett digitalt modulationssätt kallat Gaussian Minimum Shift Keying, GMSK. Mycket enkelt förklarat (och inte helt korrekt) bygger idén på att nollorna och ettorna kodas som sinusvågor med olika fas. Man kan tänka sig att en nolla kodas som A sin(T

bt), 0 ≤ t ≤ Tb, medan en etta kodas som A sin(T

bt + π2), 0 ≤ t ≤ Tb. (För en mer korrekt förklaring av GMSK och andra digitala modulationssätt hänvisas du till kurser i digital transmissionsteori.) På grund av olineariteter och spridning i komponenter i telefonen kommer mottagna nollor och ettor att bli fasförskjutna, A sin(ωt + φ0), 0 ≤ t ≤ Tb, respektive A sin(ωt + π2 + φ1), 0 ≤ t ≤ Tb, där φ:na kan ses som slumpmässiga, dvs fasförskjutningen blir olika för varje nolla och etta.

Om fasfelet råkar bli för stort så kommer telefonen att kunna göra fel när den beslutar om det var en nolla eller en etta som sänts. Därför finns det ett test som kontrollerar hur

(7)

stort fasfelet är i varje telefon. Testet går till på följande sätt. En lång känd testsekvens av nollor och ettor sänds. För varje nolla och etta som tas emot av telefonen mäts fasfelet. Kravet är att det maximalt uppmätta fasfelet för en telefon måste vara mindre än 20 grader samtidigt som det genomsnittliga absolutfasfelet (standardavvikelse eller RMS-värdet) får vara maximalt 5 grader.

I sensitivity.mat har man mätt upp mottagarkänsligheten i telefonen. Det finns näm- ligen ett krav på att vid en mottagen signalstyrka (här kallad Pin) på −102 dBm in på antennen så skall avkodaren högst ge 2 % bitfel. Enheten dBm används ofta inom radiotekniken och definieras som 10 gånger 10-logaritmen av kvoten mellan en effekt (i Watt) och 1 mW. −102 dBm motsvarar alltså

10 log10

 Pin

10−3



= −102 → Pin≈ 63 fW.

Normalt testar man just vid nivån -102 dBm och kontrollerar att bitfelen understiger 2 %, men i sensitivity.mat har man testat (för ett mindre antal telefoner) för vilken insignalnivå 2 % bitfel erhålles. Detta har gjorts både för en radiokanal i mitten på GSMs 900MHz-band, (kring 947.5 MHz) samt längst ned på 900MHz-bandet (935 MHz).

5.2 Fasfel

Ladda in filen phase_error.mat. Detta görs med kommandot load phase_error. Mät- värdena för maximalt fasfel hittas i vektorn peak_phase_error och standardavvikelsen för fasfelet finns i rms_phase_error. Börja titta på det maximala fasfelet.

>> load phase_error

>> rpe=rms_phase_error;

>> ppe=peak_phase_error;

>> x=.25:.5:30;

>> hist(ppe,x)

>> grid

Vilken fördelning kan det maximala fasfelet tänkas ha? Ser det normalfördelat ut? Här kan du t.ex. använda kommandot normplot.

>> help normplot

>> normplot(ppe) Svar: . . .

Hur många telefoner är testade totalt samt hur många klarar inte specifikationskraven?

Tag reda på detta genom att använda kommandot length och find.

>> help length

>> help find

Gör en skattning av P(En telefons maximala fasfel > 20 grader) genom att utnyttja kommandona ovan.

Svar: . . .

(8)

I fabriken kontrolleras först det maximala fasfelet. De telefoner som inte klarade kraven sorteras ut och för resten kontrollerar man standardavvikelsen på fasfelet.

Titta nu på rms_phase_error.

>> x=.1:.2:8;

>> hist(rpe,x)

Vilken fördelning kan standardavvikelsen för fasfelet tänkas ha?

Svar: . . .

Utnyttja tidigare kommandon för att uppskatta P( standardavvikelsen för fasfelet > 5 grader | maximala fasfelet < 20 grader) samt P(en telefon klarar båda fasfelstesterna).

Svar: . . .

Kan man vara nöjd med ovanstående kvalitet?

Svar: . . .

5.3 Mottagarkänslighet

Ladda nu in sensitivity.mat

>> load sensitivity

I sensitivity finns mätningar av mottagarkänsligheten på 76 telefoner för en radiokanal kring 935 MHz, kolonn 1, och för en radiokanal kring 947.5 MHz, kolonn 2. Titta på dem var för sig.

>> slc=sensitivity(:,1);

>> smc=sensitivity(:,2);

>> x=-109:.5:-104;

>> figure(1)

>> hist(slc,x)

>> figure(2)

>> hist(smc,x)

Klarar alla telefoner kravspecifikationen −102 dBm?

Svar: . . .

I teorin bör det vara så att känsligheten är sämre för kanaler nära frekvensbandets änd- punkter än kanaler i frekvensbandets mitt på grund av det bandpassfilter som filtrerar ut mottagarfrekvensbandet (935-960 MHz). Kontrollera det genom att titta på histogram- met och ”se” om det verkar som om slc är ”sämre” än smc.

Svar: . . .

(9)

Här uppkommer nu en mycket intressant fråga.: Hur ”ser” man om en datamängd skiljer sig från en annan datamängd? Man kan t.ex. testa med ett s.k. hypotestest om medel- värdena för olika datamängder skiljer sig åt. Hypotestest kommer du att lära dig lite senare i kursen.

Svar till de numrerade uppgifterna 1. a) 0.84, b) 0.81

2. 0.65

3. a) 0.30, b) 0.93

4. x = linspace(-5,9,200); % (200 värden mellan -5 och 9) y = normpdf(x,3,2);

plot(x,y)

och motsvarande för fördelningsfunktionen.

5. 66.45, 62.82, 56.74, 50.00, 43.26, 26.74.

References

Related documents

Du har rätt att efter skriftlig begäran få information om vilka personuppgifter som behandlas om dig eller ditt minderåriga barn (behöver bara vara med ifall det rör

Att dela in mat i pyramider bygger på att äta mest av det som finns i botten, lagom av det som finns i mitten och undvika eller äta lite av det som hamnar i

Därför har Hörselskadades Riksförbund (HRF) nu tagit fram en gratisapp för iPhone/iPad och Android som gör det enkelt att ta reda på om det finns tecken på

ꟷ I nästa steg får du välja ljud, klicka på Join with computer audio.. • Nu är du inne

Anledningen till detta kan vara att personen, på grund av långvarig sjukdom, funkt- ionshinder, hög ålder eller missbruk, inte klarar sin vardag på egen hand.. Vilket stöd kan

Nu har förtidsröstningen startat. Du kan välja på tre

Material: Spänningsaggregat, multimeter, dekadmotstånd, kablar och en lång kabel Rapport: Labben redovisas genom att ni svarar på frågorna i detta labb-PM och.. lämnar in

Välkommen till Arena Energiaskor för en hållbar och resurseffektiv användning av energiaskor.. 15 september 2020