I den här datorlabben kommer vi att arbeta med fördelningar i R, bland annat lära oss göra sannolikhetsberäkningar och simulera från olika fördelningar, dvs dra slumptal från fördelningar.
Installera paket
Den här labben förutsätter inga installerade paket.
Skapa mapp för labben
Skapa en mapp Lab5 i din kursmapp SDA1. Ladda ner Quarto-filen för denna lab genom att högerklicka här och välj ‘Spara länk’ eller något likande från menyn som dyker upp. Spara filen i din Lab5 mapp. Öppna Quarto-filen i RStudio och fortsätt med denna laboration direkt i Quarto-dokumentet, där du också fyller i svaren på dina laborationsövningar.
Fördelningar i R: r, p, d och q-funktionerna 😍
R har en massa statistiska fördelningar inbyggda från start, och ytterligare många, många mer i olika R-paket. För varje fördelning finns det fyra typer av funktioner:
r-funktionen som genererar slumptal från fördelningen, t ex rnorm(n = 5, mean = 2, sd = 3) genererar 5 slumptal från normalfördelningen \(X\sim N(2,3)\). Argumentet sd betyder standard deviation, dvs standardavvikelse. Tio slumptal från Poissonfördelningen \(\mathrm{Pois}(\lambda = 2)\) får man med rpois(n = 10, lambda = 2) osv. r-funktionen har fått sitt namn från random.
p-funktionen som beräknar sannolikheten att slumpvariabeln är mindre än ett visst tal, dvs \(P(X\leq x)\) för något \(x\). p-funktionen har fått sitt namn från engelskans probability. pnorm(q = 1, mean = 2, sd = 3) beräknar sannolikheten att \(X\sim N(2,3)\) är mindre än 1.
d-funktionen som för en diskret variabel \(X\) beräknar sannolikheten \(P(X=x)\). För en kontinuerlig variabel beräknar d-funktionen täthetsfunktionen \(f(x)\) i en given punkt \(x\). Det är från det kontinuerliga fallet som d-funktionen fått sitt namn, d som i density.
Kommandot dnorm(x = 0, mean = 2, sd = 3) beräknar täthetsfunktionens värde i punkten \(x=0\), dvs \(f(0)\) för en \(X\sim N(2,3)\) variabel.
Koden dpois(x = 3, lambda = 2) beräknar sannolikheten \(P(X=3)\) att den diskreta Poissonfördelade variabeln \(X\sim \mathrm{Pois}(\lambda=2)\) antar värdet \(x=3\).
q-funktionen beräknar \(p\)-kvantilen för en fördelning, dvs det värde \(x\) där \(P(X \leq x)=p\). Vi kan t ex beräkna 25% eller 0.25-kvantilen för en normalfördelad variabel \(X\sim N (2,3)\) med kommandot qnorm(p = 0.25, mean = 2, sd = 3).
Visa mig koden!
m =2# mean for the normal distributions =3# standard deviation for the normal distributionx =-1# point where the density is evaluatedx_grid =seq(-10, 10, by =0.1) # a vector of x-values for which we evaluate the density f(x)x_sample =rnorm(n =1000, mean = m, sd = s)dens_norm =dnorm(x_grid, mean = m, sd = s) # density (x) for all values in vector xgridquantile10 =qnorm(0.1, mean = m, sd = s) # 10% percentile# Plotting - a lot of code because I want pretty colors and legends. Sorry.hist(x_sample, breaks =50, freq = F, main ="X ~ N(2,3²)", ylim =c(0,0.15), xlim =c(-10,10), xlab ="x", ylab ="f(x)", col ="cornflowerblue") # freq = F gives us proportionslines(x_grid, dens_norm, lwd =3, col ="orange")density_at_x =dnorm(x, mean = m, sd = s)points(x, density_at_x)lines(x =c(x,x), c(0,density_at_x), col ="gray", lty =2, lwd =3)lines(x =c(-10,x), c(density_at_x,density_at_x), col ="gray", lty =2, lwd =2)x_lower = x_grid[x_grid < quantile10]polygon(x =c(x_lower, rev(x_lower)), border =NA,y =c(rep(0, length(x_lower)), rev(dnorm(x_lower, mean = m, sd = s))), col =rgb(1,0.8431373,0, alpha =0.7))legend("topleft", inset=.01, legend=c("histogram data", "density function", "density evaluation", "P(X<=-1)"),col=c("cornflowerblue", "orange", "gray"), lty =c(NA, 1, 3, NA), lwd =c(NA, 2, 2, NA), fill =c("cornflowerblue", NA, NA, rgb(1,0.8431373,0, alpha =0.7)), border =c(1,0,0,1), cex=1)
dnorm to compute density at -1: f(x) = 0.0806569081730478
pnorm to compute prob that X is lower than -1: P(X<=x) = 0.158655253931457
qnorm to compute 10% quantile of X as -1.8446546966338
Glöm inte Rs hjälp!
Skriv t ex ?rnorm i Console för att se hjälpfilen för normalfördelningen.
Note
Jag har skrivit ut namnen på alla funktionsargumenten ovan, t ex n=5, mean = 2 och sd = 3 i rnorm(n = 5, mean = 2, sd = 3). Det hade dock gått lika bra att skriva rnorm(5, 2, 3), då förstår R att jag menar n=5, mean = 2 och sd = 3 . Det fungerar dock bara om man använder exakt den ordning på argumenten som man ser när man skriver ?rnorm. Annars kan R inte veta vilka argument som ska matchas mot vilka siffror. Om man skriver ut argumentens namn så kan man ha vilken ordning som helst, t ex rnorm(mean = 2, sd = 3, n = 5).
Vektorisera mera
Kommandot dnorm(0, mean = 0, sd = 1) räknar ut täthetsfunktionens värde i punkten x=0 för en standard normalfördelning. Om man vill räkna ut tätheten i flera punkter? Ett sätt är att använda en for-loop som upprepar dnorm(x, mean = 0, sd = 1) för olika x-värden (se Lab1). Men det finns ett enklare sätt! Många av Rs funktioner är vektoriserade. Det betyder att du kan räkna ut funktionen för alla värden i en vektor ‘på en gång’. 😍
Så här beräknar man t ex täthetsfunktionen för tre olika x-värden -1, 0, och 1:
dnorm(x =c(-1,0,1), mean =0, sd =1)
[1] 0.2419707 0.3989423 0.2419707
Ok, snyggt, men om man vill ha olika väntevärden för var och ett av x-värdena, då? Yep, funkar:
dnorm(x =c(-1,0,1), mean =c(0, 0.5, 3), sd =1)
[1] 0.24197072 0.35206533 0.05399097
Normalfördelning med varians eller standardavvikelse?
I funktionen rnorm() anger man fördelningens spridning som en standardavvikelse med argumentet sd. Det motsvarar att man skriver \(N(\mu,\sigma)\) för en normalfördelning i matematisk/symbolisk notation. Många böcker skriver dock normalfördelningen med variansen som andra argument i den symboliska beskrivningen: \(N(\mu,\sigma^2)\). Man får helt enkelt se upp och kontrollera om en normalfördelning skrivs med en standardavvikelse eller varians. I formelsamlingen och på tentan har jag kommer jag skriva \(N(\mu,\sigma)\), som boken gör.
Rita upp sannolikhetsfördelningen för en diskret variabel
Vi ska nu försöka göra ett stapeldiagram för sannolikhetsfördelningen för en binomialfördelad variabel med parametrarna \(n=10\) och \(p=0.4\), dvs för slumpvariabeln \(X\sim \mathrm{Binom}(n,p)\). Vi har ju just sett att vi kan beräkna sannolikheterna \(P(X=x)\) för en massa olika \(x\)-värden genom vektorisering, så låt oss börja där:
n =10p =0.4xvalues =seq(0,n) # en vektor med alla möjliga utfall på X, dvs 0,1,2,...,nprobs =dbinom(x = xvalues, size = n, prob = p) # vektoriserat probs
Notera hur jag först skapade en vektor med alla x-värden (xvalues) mellan \(0\) och \(n=10\) (vilket ju är alla möjliga x-värden för en \(\mathrm{Bin}(10,p)\)-fördelad variabel). Den vektorn matade jag sen in i dbinom() funktionen för att få sannolikheten \(P(X=x)\) för varje x-värde. Vi ser att sannolikheten för \(x=0\) är väldigt låg \(P(X=0)=0.0060466176\). Den sannolikheten är lätt att räkna för hand som kontroll: det enda sätt vi kan få 0 lyckade försök i 10 Bernoulliförsök är om vi misslyckas (vars sannolikhet är \(q=1-p=1-0.4=0.6\)) på alla 10 försök: \(0.6^{10}=0.006046618\). Yes, checks out.
Koden nedan gör nu ett stapeldiagram över sannolikhetsfördelningen. Notera att vi måste använda names() funktionen för sätta ett “namn” på varje sannolikhet (dvs det x-värde som sannolikheten hör ihop med) för att vi ska få \(x\)-värdena utskrivna på \(x\)-axeln.
names(probs) <- xvaluesbarplot(probs, col ="cornflowerblue", ylab ="P(X=x)", xlab ="x", main ="X ~ Binom(10,0.4)")
Hur gör man om slumpvariabeln kan anta ett oändligt antal värden, som t ex Poissonfördelningen där \(x=0,1,2,….\) utan övre gräns? Lösningen är att låta vektorn med x-värden innehålla alla x-värden där sannolikheten \(P(X=x)\) är tillräckligt stor för att spela roll i figuren. Här kan prova sig fram, dvs öka på antalet x-värden tills sannolikheten för större x-värden är nära noll. Så här kan vi t ex göra ett stapeldiagram över sannolikhetsfördelningen för en \(\mathrm{Pois}(\lambda = 2)\) variabel:
xvalues =0:5probs =dpois(x = xvalues, lambda =2)names(probs) <- xvaluesbarplot(probs, col ="orange", xlab ="x", ylab ="P(X=x)", main ="Poisson med lite för få x-värden")
Som du ser var jag lite för snål med antal x-värden, det verkar som om x-värden större än 5 borde ha varit med i plotten eftersom sannolikheten för x=5 är inte riktigt nära noll. Låt oss prova med värden från 0 till 10:
xvalues =0:10probs =dpois(x = xvalues, lambda =2)names(probs) <- xvaluesbarplot(probs, col ="yellow", xlab ="x", ylab ="P(X=x)", main ="Poisson med lagom många x-värden")
Rita upp täthetsfunktionen för en kontinuerlig variabel
Hur kan man rita upp sannolikheterna för en kontinuerlig variabel? Kontinuerliga variabler antar ju alla värden, även decimaltal, så här måste vi göra lite annorlunda. En kontinuerlig variabel har ju också sannolikhet noll för alla \(x\): \(P(X=x)=0\). Vi använder oss därför av täthetsfunktionen \(f(x)\), där arean under täthetsfunktionen mellan två x-värden \(a\) och \(b\) är sannolikheten \(P(a \leq X \leq b)\). Den s k d-funktionen (t ex dnorm för normalfördelning) ger tätheten för en kontinuerlig variabel. För att plotta upp en täthetsfunktion:
skapa en vektor med en massa x-värden med ganska litet avstånd mellan värdena
beräkna d-funktionen för vektorn med x-värden
plotta täthetsfunktionen som en kurva med kommandot plot() eller lines() (används om det redan finns något uppritat i grafen som man vill behålla).
Här plottar vi täthetsfunktionen för en standard normalfördelningen:
xvalues =seq(-3, 3, by =0.01) # by=0.01 talar om att vi vill ha avstånd 0.01 mellan värdena.density_values =dnorm(xvalues, mean =0, sd =1)plot(xvalues, density_values, type ="l", xlab ="x", ylab ="density", main ="N(0,1)",col ="steelblue", lwd =2) # lwd är line width, dvs tjocklek på linje.
I det här fallet var det ganska lätt att välja start- och slutpunkt för x-värdena: vi vet ju från 68-95-99.7 regeln att 99.7% av sannolikhetsmassan ligger mellan -3 och 3 för standard normalfördelning. I andra fall kan det vara svårare att bestämma lämpliga x-värden och man får prova sig fram (om man har en q-funktion, dvs kvantilfunktion, för fördelningen kan man använda den för att beräkna lämpliga start- och slutvärden på x). Hur väljer man avståndet mellan olika x-värden (dvs by = 0.01 i min kod)? Det brukar spela mindre roll. Väljer man för stora avstånd blir kurvan hackig. Väljer man verkligt små avstånd blir det många olika x-värden och det kan ta tid för R beräkna alla täthetsvärden. Ett alternativ sätt att skapa denna vektor är xvalues = seq(-3, 3, length = 1000) där man har kontroll på att vektorn ska ha exakt 1000 värden mellan -3 och 3.
1. Räkna med sannolikhetsfördelningar
💪 Uppgift 1.1
Låt \(X\sim \mathrm{Pois}(\lambda)\) med \(\lambda = 2\). Beräkna \(P(X=2)\) och \(P(X\leq 3)\). [Hint: ?dpois och ?ppois]
# Write your code here
💪 Uppgift 1.2
I valet 2022 fick Liberalerna 4.61% av rösterna. I en undersökning bland totalt 1046 personser angav 30 personer att de skulle rösta på Liberalerna om det var val idag. Kan det vara så att Liberalernas röstningsandel är (ungefär) oförändrad på 4.61%, eller tyder den nya undersökningen på något annat? Undersök detta genom att rita upp ett stapeldiagram över sannolikhetsfördelningen \(P(X=x)\) för \(X \sim \mathrm{Binom}\operatorname{}(n= 1046, p = 0.0461)\) för x-värdena xvalues = 0:100. Är resultatet från valundersökningen ett sannolikt utfall om det skulle vara så att andelen i populationen faktiskt är oförändrad på 4.61%? Om inte, vilket slutsats drar du om Liberalernas faktiska väljarandel i populationen baserat på utfallet från undersökningen?
# Write your code here
💪 Uppgift 1.3
Beräkna \(P(X \leq 30)\) för \(X \sim \mathrm{Binom}(n= 1046, p = 0.0461)\) genom att använda binomialfördelningen. Jämför svaret med samma sannolikhet från en normalapproximation (se F14):
\[
X\sim \mathrm{Binom}(n,p) \text{ approximeras med } X\sim \mathrm{N}\Big(\mu=np, \sigma = \sqrt{n p (1-p)}\Big)
\]
Observera att jag skrivit normalfördelningen med standardavvikelse som andra argument här. Precis som R gör.
# Write your code here
2. Simulera slumpvariabler
💪 Uppgift 2.1
Simulera nu 1000 slumptal från Poissonfördelningen \(\mathrm{Pois}(\lambda)\) med \(\lambda = 2\) och spara i variabeln (vektorn) x. Beräkna medelvärdet (mean(x)) och standardavvikelsen (sd(x)) för slumptalen. Stämmer de värdena ungefär med vad vi kan förvänta oss? [hint: Poissonfördelningens teoretiska väntevärde]. Simulera nu 10000 slumptal från samma fördelning och beräkna återigen medelvärdet. Vad tror du skulle hända om du fortsatte så här och simulerade fler och fler slumptal? (bortsett från att datorns minne skulle ta slut). [Hint: stora talens lag].
# Write your code here
💪 Uppgift 2.2
Antag att du inte känner till funktionen dpois och att det matematiska uttrycket för Poissonsannolikheterna som ges i SDM boken är för skrämmande för dig. Kan du använda slumptalen i Uppgift 2.1 för att uppskatta sannolikheten \(P(X=3)\)?
# Write your code here
💪 Uppgift 2.3
Simulera 1000 slumptal från exponentialfördelningen med parametern rate = 2 (vilket är det som boken kallar \(\lambda\)) och spara i en vektor x. Beräkna medelvärdet av slumptalen och verifiera att det är hyfsat nära det teoretiska väntevärdet för denna exponentialfördelning.
# Write your code here
💪 Uppgift 2.4
Rita ett histogram för slumptalen från förra uppgiften med argumentet breaks = 40 (gör att du får ungefär 40 staplar i histogrammet). Plotta också den teoretiska täthetsfunktionen för \(\mathrm{Expon}(\lambda = 2)\) i samma figur.
Plottips
Använd lines() när du plottar täthetsfunktionen, annars kommer histogrammet att skrivas över. Använd argumentet freq = FALSE i hist()-funktionen så du får andelen observationer i varje bin, och inte antalet. Annars kommer histogrammet och täthetsfunktionen inte att ha samma skala på y-axeln.
# Write your code here
3. Simulering för att beräkna nya fördelningar
En mycket trevlig egenskap med simulering är att det är väldigt enkelt att beräkna fördelningen för en funktion av slumpvariabler. Vi vet ju t ex att om \(X\sim \operatorname{N}(\mu,\sigma)\) så är fördelningen för linjärkombinationen \(Y=c+aX\) också normalfördelad: \(Y\sim \operatorname{N}(c+a\mu, \vert a \vert \sigma )\). Men det är långt ifrån alla fall där vi matematiskt kan bevisa sådana goa resultat. Med simulering kan vi få fram fördelningen för alla funktioner och även fördelningen för en summa av slumpvariabler, och mer. Den kan vi även göra om fördelningen för \(X\) är något annat än normalfördelad. Låt oss testa direkt i en övning!
💪 Uppgift 3.1
Din vän och du ska beställa hem mat. Ni har olika smak och du vill beställa från restaurang X, din vän vill beställa från restaurang Y. Ni vill såklart gärna äta maten samtidigt och lovar varandra att vänta med att äta tills båda fått maten. Låt \(X \sim \mathrm{Expon}(\lambda = 2)\) vara väntetiden i timmar tills restaurang X levererar mat hem till dig, och \(Y \sim \mathrm{Expon}(\lambda = 4)\) är väntetiden för restaurang Y. Antag att \(X\) och \(Y\) är oberoende slumpvariabler. Slumpvariabeln \(Z=60(X-Y)\) mäter då hur många minuter du måste vänta på maten efter att din vän redan fått maten, dvs hur länge din artiga vän måste vänta innan hen kan börja äta. Om \(Z\) är negativ har du alltså fått maten före din vän. Här ser man direkt att \(Z=60(X-Y)\) inte kan vara exponentialfördelad, \(Z\) kan ju vara negativ! Använd simulering för att simulera från fördelningen för \(Z\) och rita upp ett histogram för att representera täthetsfunktionen (använd argumentet freq = FALSE för att få andelar inom varje bin). Använd simuleringen för att uppskatta sannolikheten att du får maten efter din vän. [Hint1: Prova att skriva Z>0 i Console. Hint2: R hanterar TRUE som 1 och FALSE som 0, så t ex sum(c(TRUE, FALSE, TRUE)) blir 2 i R].
# Write your code here
💪 Uppgift 3.2
Uppskatta variansen för \(Z\) genom simulering. Verifiera att uppskattningen är hyfsat nära den teoretiska variansen för \(Z\). [Hint: F13 gav oss formler för variansen av en summa och en differens av oberoende variabler och F15 har egenskaper för exponentialfördelade variabler, t ex deras varians.]
# Write your code here
4. Analys av antalet besök till webbsida
En webbsida registrerade när besökare besökte deras webbsida under ett dygn. Här är en grafisk illustration av antalet besök under dygnet, varje streck är ett besök vid en viss tidpunkt:
Antalet besök för var och ett av dygnets timmar var
Antag att antalet besök per timme kan modelleras som oberoende Poissonfördelade variabler med parameter \(\lambda\). Vi vill nu bestämma det bästa värdet på \(\lambda\) genom att använda data, dvs vi vill skatta (estimera) modellparametern \(\lambda\). En vanlig skattning av \(\lambda\) i en Poissonfördelning är medelvärdet \(\bar x\) av stickprovet, dvs medelvärdet av variabeln visits_per_hour (det är t ex den skattning som är optimal enligt den s k maximum likelihoodmetoden).
💪 Uppgift 4.1
Skatta parametern \(\lambda\) från variabeln visits_per_hour .
# Write your code here
💪 Uppgift 4.2
Företaget som driver webbsidan vill veta sannolikheten att de får fler än 5 besökare under åtminstone någon av morgondagens 24 timmar. Använd din skattade Poissonmodell för att beräkna den sannolikheten. [hint: union och komplement].
# Write your code here
💪 Uppgift 4.3
Undersök om Poissonfördelningen verkar anpassa dessa data bra. Du kan t ex undersöka om det teoretiska väntevärdet och variansen i den skattade Poissonmodellen från Uppgift 4.1 ligger hyfsat nära medelvärdet och stickprovsvariansen i data. Du kan också kolla hur sannolikt det faktiskt är att observera 9 st besök under en timme (som vi ju ser i data mellan kl 17-18) i den skattade Poissonmodellen.
