Datorlaboration: Del1

Uppgiftsförklaring

Uppgiften går ut på att samla ihop och sammanställa 4 uppgifter från tidigare labbar så att de går att läsa och kan förstås självständigt utan andra hjälpmedel.

Uppgift 1

I uppgift 1 undersöker vi sammansatta sannolikheter för populationer på nationell nivå. Vi jämför även teoretiska värden mot exakta uträkningar och slutligen räknar vi ut sannolikheten att två personer i en grupp, att det finns två personer med samma födelsedag inom gruppen.

Uppgift 1a)

Uppgift 1a går ut på att beräkna exakt sannolikhet att två slumpvist valda USA-medborgare har samma blodtyp. Sannolikheten för en slumpmässigt utvald medborgare lyder: A: 0.40 AB:0.04 B: 0.11 O: 0.45

Sannolikheterna multipliceras med varandra för att få ut sannolikheten att två personer har samma blodgrupp.

a <- c("A", "AB", "B", "O")
b <- c(0.40, 0.04, 0.11, 0.45)

totallength <- 0 # initialize
totallength[1:(length(b))] <- 0 # initialize
totallength[1] <- b[1]*b[1]
totallength[2] <- b[2]*b[2]
totallength[3] <- b[3]*b[3] 
totallength[4] <- b[4]*b[4] 

table(a,b)

##     b
## a    0.04 0.11 0.4 0.45
##   A     0    0   1    0
##   AB    1    0   0    0
##   B     0    1   0    0
##   O     0    0   0    1

print(totallength)

## [1] 0.1600 0.0016 0.0121 0.2025

sam <- sum(totallength)
print(sam)

## [1] 0.3762

I den här delen beräknar vi sannolikheten från ett replikat på 10000 simulerade “deltagare”, där alla simuleringar har samma initiella parametrar (sannolikheter). Och sedan studerar vi utfallet av simuleringen. Desto fler simuleringar vi tar med, destor närmare kommer vi det uträknade värdet ovan.

Det går säkert att göra slututräkningen på ett elegantera sätt, men då vi har 10000 replikat kommer slutsumman att bli 10000*10000 måste vi dela slutprodukten för att få den på samma form som ovan.

a <- c("A", "AB", "B", "O")
b <- c(0.40, 0.04, 0.11, 0.45)
sample(x = a, size = 2, prob = b, replace = TRUE)

## [1] "O" "O"

sim_data <- sample(
  x = a, size = 10000,
  prob = b, replace = TRUE
)
table(sim_data)

## sim_data
##    A   AB    B    O 
## 3943  392 1136 4529

table(sim_data) / 10000

## sim_data
##      A     AB      B      O 
## 0.3943 0.0392 0.1136 0.4529

totallength <- 0 # initialize
totallength[1:(length(b))] <- 0 # initialize
totallength[1] <- table(sim_data)[1]*table(sim_data)[1]
totallength[2] <- table(sim_data)[2]*table(sim_data)[2]
totallength[3] <- table(sim_data)[3]*table(sim_data)[3] 
totallength[4] <- table(sim_data)[4]*table(sim_data)[4] 


print(totallength/100000000)

## [1] 0.15547249 0.00153664 0.01290496 0.20511841

sim <- sum(totallength/100000000)
print(sim)

## [1] 0.3750325

Och då kan man jämföra det simulerade värdet i ‘sim’ och mot det uträknade värdet tidigare i ‘sam’.

Uppgift 1b)

Uppgift 1b) går ut på att beräkna sannolikheten den exakta sannolikheten att det finns två personer med samma födelsedag. Det lättaste sättet att göra detta på var att loopa resultatet, då sannolikheterna bygger på varandra. Vi bygger loopen på antagandet att det inte finns någon matchning, och går stegvis igenom sannolikheten för varje klasskamrat.

k = 25  # klasskamrater
p <- numeric(k)  # sannolikheter
for (i in 1:k)      {
            q <- 1 - (0:(i - 1))/365  # 1 - (sannolikhet(inga matchningar))
            p[i] <- 1 - prod(q)  }
prob <- p[k]

Här ser vi beräknad sannolikhet och jämför det med det inbyggda kommandot i R för jämförelse av födelsedagar.

prob

## [1] 0.5686997

pbirthday(25)

## [1] 0.5686997

Och avslutande var uppgiften att ta i beaktning skottår. Enda skillnaden blir att vi tar med årets “exakta” dagsantal, istället för det avrundade värdet vi använde tidigare.

k = 25  # klasskamrater
p <- numeric(k)  # sannolikheter
for (i in 1:k)      {
            q <- 1 - (0:(i - 1))/365.2425  # 1 - sannolikhet(inga matchningar)
            p[i] <- 1 - prod(q)  }
prob <- p[k]
prob

## [1] 0.5684532

Uppgift2

I uppgift 2 undersöker vi en sannolikhets-massfunktion, tar fram dess kvantilfunktion och utifrån sagd funktion dess median. Sedan utför vi en simulering på lynx populationen med uträknade värden och jämför simulering mot uträkningar. Slutligen plottar vi en enkel hypergeometrisk fördelning.

Uppgift 2a)

I av Speegles bok så ges en fördelning för en variabel \(X\) som anger antalet ungar till lodjur. Representera denna fördelning med en tabell som ovan och plotta en “stemplot”. Vi tar värdena ur boken nedan och visualiserar deras sannolikheter i en plot nedan.

df <- tibble(x=0:5,p=c(0,0.18, 0.51, 0.27, 0.04,0))
head(df,4)

## # A tibble: 4 × 2
##       x     p
##   <int> <dbl>
## 1     0  0   
## 2     1  0.18
## 3     2  0.51
## 4     3  0.27

ggplot(aes(x=x,y=p), data=df) + 
    geom_segment(aes(xend=x,yend=0)) + 
    geom_point(size=3,shape="circle",color="red") +
    ggtitle("Stemplot av PMF")

Uppgift 2b)

Plotta kvantilfunktionen \(F^{-1}(p)\)

Gör en plot av både CDF \(F(x)\) och kvantilfunktionen \(F^{-1}(p)\) för Lynx-fördelningen ovan. Vi använder funktionen cumsum, som egentligen betyder kummulativ summa, på sannolikheterna vi använde tidigare i uppgiften.

Avläs från plotten av kvantilfunktionen \(F^{-1}(x)\) ett lämpligt värde för fördelningens median. Medianen är enkelt utläst mitt på denna plot.

dF <- mutate(df,cdf=cumsum(p))
dF

## # A tibble: 6 × 3
##       x     p   cdf
##   <int> <dbl> <dbl>
## 1     0  0     0   
## 2     1  0.18  0.18
## 3     2  0.51  0.69
## 4     3  0.27  0.96
## 5     4  0.04  1   
## 6     5  0     1

ggplot(dF,aes(x,y=cdf)) + 
  geom_step() + 
    ggtitle("Plot av CDF") +
    labs(y="F(x)",x="Variabel X")

ggplot(dF,aes(x=cdf,y=x)) + 
  geom_line() + 
    ggtitle("Plot av Kvantilerna") +
    labs(y="X vardet",x="Sannolikhet p")

Medianen blir cirka 1.62

Uppgift 3c)

Lynx simulering.

Låt \(X\) ange kullstorleken hos lodjur som ovan. Simulering med \(N=1000\) replikat värdena av \(X^2\). Plotta data i ett histogram med geom_bar(). Beräkna medelvärdet av observationerna av \(X^2\) och jämför med det teoretiska värdet av \(E(X^2)\) som beräknas som ovan.

Vi gör beräkningar för att räkna ut meddelvärde, varians och standardavvikelse och jämför det teoretiska med det simulerade värdet på \(E(X^2)\) .

sumdf <- df %>% 
  summarise(m1 = sum(x*p), m2=sum(p*x^2), varians=m2-m1^2, stdavvikelse=sqrt(varians))
sumdf

## # A tibble: 1 × 4
##      m1    m2 varians stdavvikelse
##   <dbl> <dbl>   <dbl>        <dbl>
## 1  2.17  5.29   0.581        0.762

litterpmf <- c(0,0.18, 0.51, 0.27, 0.04,0)
sum((0:5)^2 * litterpmf)

## [1] 5.29

X <- sample((0:5)^2, 1000, replace = TRUE, prob = litterpmf)
mean(X)

## [1] 5.141

##Uppgift 3d) Plotta den hypergeometriska fördelningen

Och avslutningsvis på denna uppgiften skulle vi plotta den hypergeometriska fördelningen med nedanstående parametrar. Plotta CDF för den hypergeometriska fördelningen som ger antalet röda kuler ett urval på 13 kulor utan återläggning från en urna med 30 kulor varav 5 är röda.

x_phyper <- seq(0, 30, by = 1) 
y_phyper <- phyper(x_phyper, m = 13, n = 17, k = 5, log=FALSE) 
plot(y_phyper,) 

Uppgift 3

I uppgiften plottar vi två normalfördelade variabler och kämför deras pdf. Vi undersöker även dess medelvärde, median och undersöker dess kvantiler. Slutligen tittar vi på skillnaden mellan värdet på beräknad varians och simulerad varians med två olika fördelningar.

Uppgift 3a)

Plotta täthetsfunktionerna för två normalfördelade variabler.

I den här uppgiften går det bra att använda vilka normalfördelade variabler man än vill använda, man får bara se till att justera grafen så att båda kurvor syns tydligt. Vi använder dnorm() som är en standardfunktion för att ta fram pdf:en av en normalfördelad fördelning.

x1 <- seq(-15, 18, length=1000)
y1 <- dnorm(x1, mean=1, sd=1)

x2 <- seq(-15, 18, length=1000)
y2 <- dnorm(x2, mean=2, sd=5)

plot(x2,y2,type="l",col="red",ylim = c(0, .4), xlab='x', ylab='Sannolikhet')
lines(x1,y1,col="green")

Upg 3b)

Uppskatta \(\mu\) av X, medianen m och IQR

\(\mu\) =1, m=1 p.g.a symmetri, det är där största densiteten ligger samt vart “massan” är fokuserad runt.

IQR i figuren ser vi att i intervallet (-0.1, 2.1) ligger ungefär 50 procent av massan, via en grov, visuell, skattning.

Kvantil1: Q1=-0.1

Kvantil3: Q3=2.1

Upg 3c)

För a=1 och b=3 plotta pdf och cdf för U.

Här beräknar vi en likformigt kummulativ plot och en likformig plot för att se hur de skiljer sig åt, med redan givna parametrar. Dunif() ger den kontinuerliga likformiga tätheten. Punif() ger den kontinuerligt likformiga fördelningen.

x <- seq(0, 4, length=1000)
y <- dunif(x, min = 1, max = 3)
plot(x, y, type = 'l', lwd = 3, ylim = c(0, .5), col='blue',
     xlab='x', ylab='Sannolikhet', main='Likformig density plot')

x_punif <- seq(0, 5, length=1000)
y_punif <- punif(x_punif, min = 1, max = 3) 
plot(x_punif, y_punif, type = 'l', lwd = 3, ylim = c(0, 1), col='blue',
     xlab='x', ylab='Sannolikhet', main='Likformigt kummulativ densitets plot')

Upg 3d)

Verifiera genom simulering att U har varians Var(U) som ges av:

sigma2 = (b-a)^2 / 12

Vi har redan given uträkning för varians ovan, nedan gör vi beräkning för en simulering med runif() som ger slumpmässiga tal på den likformiga fördelningen. 10000 replikat. Vi gör detta för med två olika upptsättningar av parametrar. mean(c^2)-mean(c)^2 ger värdet på variansen på formen: \(E(X^2)\)-\(E(X)^2\)

Uträkning samt simulering för b=3, a=1, b-a=2

c <- runif(10000,1,3)
mean(c^2)-mean(c)^2

## [1] 0.3400615

b=3
a=1
sigma2 = ((b-a)^2) / 12
sigma2

## [1] 0.3333333

Uträkning samt simulering för b=7, a=5, b-a=2

c <- runif(10000,5,7)
mean(c^2)-mean(c)^2

## [1] 0.3305095

b=7
a=5
sigma = ((b-a)^2) / 12
sigma

## [1] 0.3333333

Vi ser att vi får tal som är väldigt nära varandra, och vi kan alltid öka antalet simuleringar för att komma närmare det exakta värdet.

Uppgift 4

Uppgiften går ut på att betrakta en betingad fördelning för att ta fram en punktgraf som vi sedan använder, omformar och förlänger för att se de olika uttryckens fördelningar grafiskt. I slutet på uppgiften använder vi lagen om total varians och jämför skattningen mot ett uträknat värde. Slutligen jämför vi korrelation mellan tre variabler och den betingade fördelningens residual.

Uppgift 4a) - Betingad fördelning

Här studerar vi en enkel modell med två variabler \((X,Y)\) där \(P(X)=\mathrm U[0,6]\) och

\[P(Y|X)=\mathrm N(\mu=6+1.3\cdot \sin(\pi x/2),\sigma=0.4+(3-x)^2).\]

Vi har att \(P(X)\) är likformigt fördelad mellan a=0, b=6.

\(P(Y|X)\) är normalfördelad med väntevärde: 6+1.3sin(pi*x/2) och std=0.4+(3-x)^2

Först samplar vi \(n=500\) observationer, vi lägger in den data vi fick ovan och sedan beräknar vi det betingade väntevärdet och residualen Y-E(Y|X). Detta plottas med en punktgraf för varje N och för varje värde.

Jag slår även samman uppgiften som ursprungligen var upg c) här i slutet, för enkelhetens skull, där jag genomför en “density-plot” som ger de “empiriska densiteterna” för variablerna \(Y\), \(E(Y|X)\) och \(Y-E(Y|X)\). Använder ggplot med geom_density(). Melt används på ‘obsext’ som vi tog fram i i början på den här uppgiften för att omforma och förlänga datan. Vi sammanväver helt enkelt datan för de tre variabler från punktgrafen i 4a.

n <- 500

# Betingat väntevärde och betingad varians (standardavvikelse)
mf <- function(x) 6+1.3*sin(pi*x/2)
sdf <- function(x) 0.4+(3-x)^2
obs <- tibble(
  X=runif(n,min=0,max=6),
  Y=rnorm(n,mean=mf(X),sd=sdf(X))
)

# Lägg till betingade väntevärdet och residualen
obsext <- obs %>% mutate(
  `E(Y|X)`=mf(X),
  `Y-E(Y|X)`=Y-mf(X)
  ) 


obsext %>%
  pivot_longer(cols=c("Y","E(Y|X)","Y-E(Y|X)")) %>%
  ggplot(aes(x=X,y=value,color=name)) +
  geom_point(size=1) +
  scale_colour_manual(values=c("Y" = "red","E(Y|X)" = "blue","Y-E(Y|X)" ="black")) + 
  ggtitle("En fördelning av (X,Y) med sinusfunktion som betingat väntevärde av Y") +
  labs(y="Y",color="")

obsext.m <- melt(obsext)

## No id variables; using all as measure variables

ggplot(obsext.m, aes(value, fill=variable)) +
     geom_density(alpha=0.4)

Uppgift 4b)

Beräkna variansen av \(Y\) med hjälp av lagen om total varians. Kontrollera genom att jämföra med den skattade varians på de simulerade värdena.

Lagen om total varians säger att:

Var(Y ) = E (Var(Y |X)) + Var (E(Y |X))

E(Var(Y|X))=Var(Y-E(Y|X)) - Vilket är första termen för att beräknar Var(Y)

Vi ser variansen, var(Y) i tabellen nedan, där vi utför ovanstående beräkningar.

  zx <- obs %>% mutate(
`E(Y|X)`=mf(X),
`Y-E(Y|X)`=Y-mf(X),
`var(Y)` = var(`Y-E(Y|X)`)+var(`E(Y|X)`)
)
  zx

## # A tibble: 500 × 5
##        X     Y `E(Y|X)` `Y-E(Y|X)` `var(Y)`
##    <dbl> <dbl>    <dbl>      <dbl>    <dbl>
##  1  5.46 17.2      6.98    10.2        17.7
##  2  2.74  4.77     4.81    -0.0382     17.7
##  3  4.98  7.35     7.30     0.0525     17.7
##  4  3.09  4.74     4.71     0.0228     17.7
##  5  1.03  1.80     7.30    -5.49       17.7
##  6  3.00  5.32     4.70     0.616      17.7
##  7  5.70 16.8      6.59    10.2        17.7
##  8  5.23 17.2      7.22     9.99       17.7
##  9  2.60  4.13     4.95    -0.822      17.7
## 10  1.62  6.13     6.73    -0.599      17.7
## # … with 490 more rows

Uppgift 4c)

Visa med hjälp av tre exempel att variabeln \(R=Y-E(Y|X)\) är okorrelerad med variabler av typen \(W=h(X)\). Inför tre nya variabler som funktioner av \(X\), exempelvis \(X\) , \(X^2\) och \(X^3\). Använd simulerade data och skattning med hjälp av cov(x,y).

Okorrelerade om Cov(R,W)=blir nära noll. Och vi ser att vår kovarians blir väldigt liten så de har i stort sett ingen korrelation.\

n <- 20

# Betingat väntevärde och betingad varians (standardavvikelse)

  zd <- obs %>% mutate(
`E(Y|X)`=mf(X),
 `R`=Y-mf(X),
`var(Y)` = var(`R`)+var(`E(Y|X)`)


 

)
 
   ds <- zd %>% 
     mutate(
                       w1=X,
                       w2=X^2,
                       w3=X^3 
          
           
            )%>%
    summarise(cov(`R`,w1),cov(`R`,w2),cov(`R`,w3),)
   ds

## # A tibble: 1 × 3
##   `cov(R, w1)` `cov(R, w2)` `cov(R, w3)`
##          <dbl>        <dbl>        <dbl>
## 1       -0.269        -1.80        -10.8

Och slutligen gör vi en jämförelse mellan tre variabler som vi skattar, och jämför det med residualen. Vi ser att för de tal vi valt är de i stort sett okorrelerade, men den korrelation de har är alltid negativ.