Rozdelenia pravdepodobnosti

Pre väčšinu rozdelení vieme v R počítať hodnoty nasledovných funkcií (* označuje rozdelenie): - d*: hustota - p*: distribučná funkcia - q*: kvantilová funkcia - r*: realizácia náhodného výberu

Teda napríklad pre normálne rozdelenie so strednou hodnotou mean a smerodajnou odchylkou sd máme:

dnorm(x, mean, sd) #hustota v bode x

pnorm(q, mean, sd) #distribucna funkcia v bode q

qnorm(p, mean, sd) #kvantil v bode p

rnorm(n, mean, sd) #vygenerovanie n nahodnych realizacii z rozdelenia

x <- seq(-10, 10, by=.05)
y <- dnorm(x)
plot(x, y, type="l", col="red", lwd=2)

y <- dnorm(x, mean=2.5, sd=3)
plot(x, y, type="l", col="blue", lwd=2)

y <- pnorm(x, mean=3, sd=4)
plot(x, y, type="l", col="green", lwd=2)

Ďalšie rozdelenia fungujú rovnako, napr.:

Príklad: Biased Coin Design

Testujeme vedľajšie účinky lieku nasledovným spôsobom: prvému pacientovi dáme najnižšiu dávku, akú máme k dispozícii. Ak mal pacient \(i\) nejaké vedľajšie účinky, pacientovi \(i+1\) dávku znížime (ak môžeme). Ak pacient \(i\) nemal žiadne vedľajšie účinky, rozhodujeme sa podľa výsledku hodu mincou: ak padne hlava, pacient \(i+1\) dostane rovnakú dávku ako pacient \(i\), ak padne znak, pacientovi \(i+1\) dávku zvýšime (ak môžeme).

Napíšte program, ktorého vstupom bude vektor pravdepodobností vedľajších účinkov \(p\) pre jednotlivé dávky, pravdepodobnosť padnutia hlavy na minci \(b\) a počet pacientov \(n\) a ktorého výstupom bude vektor dĺžky \(n\) s dávkami lieku pridelenými jednotlivým pacientom. Predpokladáme, že máme k dispozícii dávky \((1,2,\ldots,k)\).

doses <- function(p,b,n){
k <- length(p) #dlzka vektora p
dose <- rep(0,n) #inicializacia vektora davok
dose[1] <- 1 #zaciname najnizsou davkou
for (i in 2:n){
  u <- dose[i-1] #davka pacienta i-1
  tox <- rbinom(1, size=1, prob=p[u]) #pravd. ze pacient i-1 mal vedlajsie ucinky
  head <- rbinom(1, size=1, prob=b) #pravd. padnutia hlavy
  if (tox==1){ #pacient i-1 mal vedlajsie ucinky
    if(dose[i-1]!=1) dose[i] <- dose[i-1] - 1 #pacient i-1 nedostal najnizsiu moznu davku
    else dose[i] <- dose[i-1]
  }
  if (tox==0){ #pacient i-1 nemal vedlajsie ucinky
    if (head==0){
      if (dose[i-1]!=k) dose[i] <- dose[i-1] + 1#pacient i-1 nedostal najvyssiu moznu davku
      else dose[i] <- dose[i-1]
    } else dose[i] <- dose[i-1]
  }
}
return(dose)
}
doses(c(0.1, 0.2, 0.3, 0.4), 0.5, 6)
## [1] 1 1 1 1 1 2

Viacrozmerné normálne rozdelenie

Pre dvoj- a viacrozmerné rozdelenia už väčšinou nemáme podobné jednoduché funkcie ako pre jednorozmerné rozdelenia. Výnimkou je viacrozmerné normálne rozdelenie, z ktorého môžeme generovať pomocou funkcie rmvnorm v balíku mvtnorm:

library(mvtnorm)
## Warning: package 'mvtnorm' was built under R version 4.1.1
x <- rmvnorm(n=50000, mean=c(1,2), sigma=diag(2))
colMeans(x)
## [1] 0.9941614 1.9981545
var(x)
##              [,1]         [,2]
## [1,]  1.013472656 -0.004656853
## [2,] -0.004656853  0.998747864
plot(x, col=rainbow(10), pch=16)

Vidíme, že sa vygeneroval mrak bodov, ktorý má najväčšiu hustotu v okolí strednej hodnoty. Výberový priemer a výberová kovariančná matica sa “podobá” tým, ktoré sme zadali ako mean a sigma, ale samozrejme nie sú úplne rovnaké, keďže ide o náhodné generovanie. Samotnú hustotu môžeme vizualizovať neasledovne:

library(MASS)
fit2 <- kde2d(x[ ,1], x[ ,2])
persp(fit2$x, fit2$y, fit2$z,col="magenta")

persp3d(fit2$x, fit2$y, fit2$z, theta=45, phi=25, col="blue")

Základné štatistické funkcie

Uvažujme dáta, v ktorých sa merala veľkosť chodidla a výška chlapcov a dievčat rôzneho veku. Pred analýzou z nich najskôr odstránime chýbajúce pozorovania, nastavíme pohlavie (gender) ako faktorovú premennú a premenujeme jej úrovne.

foot <- read.csv("http://www.iam.fmph.uniba.sk/ospm/Filova/rm/anthropometry.csv")
head(foot)
##        age gender foot_length height
## 1 4.166667 female         163  103.3
## 2 4.250000   male         163  103.9
## 3 4.416667 female         171  111.2
## 4 3.833333 female         163   99.7
## 5 3.416667 female         167   99.7
## 6 3.666667   male         164   98.7
foot <- foot[complete.cases(foot),]
foot$gender <- as.factor(foot$gender)
levels(foot$gender)
## [1] "female" "male"
levels(foot$gender) <- c(1, 2) 
summary(foot)
##       age         gender    foot_length        height     
##  Min.   : 2.000   1:1880   Min.   :116.0   Min.   : 81.3  
##  1st Qu.: 6.583   2:1951   1st Qu.:185.0   1st Qu.:118.6  
##  Median :10.750            Median :220.0   Median :141.3  
##  Mean   :10.443            Mean   :214.3   Mean   :138.8  
##  3rd Qu.:14.000            3rd Qu.:240.0   3rd Qu.:158.8  
##  Max.   :20.000            Max.   :311.0   Max.   :194.4

Ďalej si dáta vykreslíme, aby sme získali lepšiu predstavu o ich priebehu:

boxplot(foot_length ~ gender, data=foot, col=foot$gender)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.1

t-test

Prirodzene nás môže zaujímať, či je významný rozdiel medzi dievčatami a chlapcami. Túto hypotézu môžeme otestovať (za predpokladu normality) pomocou t-testu:

t.test(foot_length ~ gender, data=foot)
## 
##  Welch Two Sample t-test
## 
## data:  foot_length by gender
## t = -8.1957, df = 3706.3, p-value = 3.397e-16
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -11.760221  -7.219743
## sample estimates:
## mean in group 1 mean in group 2 
##        209.4654        218.9554
t.test(height ~ gender, data=foot)
## 
##  Welch Two Sample t-test
## 
## data:  height by gender
## t = -3.2316, df = 3801.4, p-value = 0.001241
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -4.141379 -1.013804
## sample estimates:
## mean in group 1 mean in group 2 
##        137.4496        140.0272

Poznámka: na test hypotézy, že viacrozmerné pozorovanie (napríklad výška a dĺžka chodidla súčasne) sa medzi dvoma skupinami líši, používame tzv. Hotellingov T-test, s ktorým sa stretnete na magisterskom štúdiu.

Kovariančná a korelačná matica

Všetky tri spojité pozorovania (vek, výšku a dĺžku chodidla) pre dané dieťa môžeme považovať za realizáciu trojrozmerného náhodného vektora, a teda môžeme počítať výberovú kovariančnú a výberovú korelačnú maticu.

Výberovú kovariančnú maticu pre náhodný výber rozsahu \(n\) definujeme ako

\[ S^2=\frac{1}{n-1}\sum_{i=1}^n (x_i-\bar{x})(x_i-\bar{x})^T \]

Na diagonále máme teda výberové disperzie a mimo diagonály výberové kovariancie.

Výberová korelačná matica má prvky

\[ \{R\}_{ij}=\frac{s_{ij}}{\sqrt{s_{ii}s_{jj}}} \]

V <- var(foot[ , c(1,3,4)])
V
##                   age foot_length   height
## age          19.55473    145.2116 104.6565
## foot_length 145.21158   1316.3112 868.4350
## height      104.65645    868.4350 613.3581
R <- cor(foot[ , c(1,3,4)])
R
##                   age foot_length    height
## age         1.0000000   0.9050990 0.9556152
## foot_length 0.9050990   1.0000000 0.9664981
## height      0.9556152   0.9664981 1.0000000
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.1.2
## corrplot 0.92 loaded
corrplot.mixed(R, lower = "number", lower.col = "black", upper = "ellipse", tl.pos="lt")

Lineárna regresia

Na lepšie popísanie vzťahu medzi premennými použijeme lineárnu regresiu. V lineárnej regresii chceme vyjadriť závislosť jednej vysvetľovanej premennej od jednej alebo viacerých vysvetľujúcich premenných.

Jednoduché vizuálne vysvetlenie je napríklad na http://shiny.webpopix.org/sia/polyreg/ alebo v balíku animation:

library(animation)
example(put.points.demo)
example(least.squares)

Skúsime najskôr fitovať model závislosti výšky od veku:

fit <- lm(height ~ age, data=foot)
summary(fit)
## 
## Call:
## lm(formula = height ~ age, data = foot)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -29.7533  -4.4534  -0.0474   4.4605  29.1365 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 82.87175    0.30240   274.1   <2e-16 ***
## age          5.35198    0.02667   200.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.297 on 3829 degrees of freedom
## Multiple R-squared:  0.9132, Adjusted R-squared:  0.9132 
## F-statistic: 4.028e+04 on 1 and 3829 DF,  p-value: < 2.2e-16
confint(fit)
##                 2.5 %    97.5 %
## (Intercept) 82.278873 83.464619
## age          5.299697  5.404257
t <- fit$coef
plot(foot$age, foot$height, pch=20, col=foot$gender)
abline(coef=t, col="blue", lwd=2)

V nasledovnom modeli ukážeme, ako závisí dĺžka chodidla od veku a výšky:

fit2 <- lm(foot_length ~ age + height, data=foot)
summary(fit2)
## 
## Call:
## lm(formula = foot_length ~ age + height, data = foot)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.379  -5.655  -0.304   5.735  43.592 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -5.314      1.699  -3.127  0.00178 ** 
## age           -1.749      0.112 -15.612  < 2e-16 ***
## height         1.714      0.020  85.708  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.032 on 3828 degrees of freedom
## Multiple R-squared:  0.9381, Adjusted R-squared:  0.938 
## F-statistic: 2.899e+04 on 2 and 3828 DF,  p-value: < 2.2e-16
open3d()
## null 
##    3
x <- foot$foot_length
y <- foot$age
z <- foot$height
rgl.spheres(x, y, z, r = 0.8, color = foot$gender) 
aspect3d(1,1,1)
axes3d()
fit <- lm(z ~ x + y)
coefs <- coef(fit)
a <- coefs["x"]; b <- coefs["y"]; c <- -1; d <- coefs["(Intercept)"]
rgl.planes(a, b, c, d, alpha=0.5, color = "blue")

Samozrejme má zmysel uvažovať nad tým, či nezostaviť dva rôzne modely, zvlášť pre chlapcov a zvlášť pre dievčatá, najmä keď sme pomocou t-testu zistili, že je medzi nimi signifikantný rozdiel:

foot1 <- foot[foot$gender==1, ]
foot2 <- foot[foot$gender==2, ]
fit1 <- lm(foot_length ~ age, data=foot1)
fit2 <- lm(foot_length ~ age, data=foot2)
(t1 <- fit1$coef)
## (Intercept)         age 
##  142.495749    6.418219
(t2 <- fit2$coef)
## (Intercept)         age 
##  131.666843    8.351928
plot(foot$age, foot$foot_length, pch=20, col=foot$gender)
abline(coef=t1, col=1, lwd=2)
abline(coef=t2, col=2, lwd=2)

Bonus: Chíkvadrát test dobrej zhody

V športovom klube sa pýtali náhodne vybraných ľudí, ktorého kandidáta budú voliť za predsedu klubu. 35 opýtaných chce voliť kandidáta A, 24 kandidáta B, 17 kandidáta C a 19 kandidáta D. Otestujte H0, ze všetci kandidáti majú rovnakú podporu.

chisq.test(c(35,24,17,19))          
## 
##  Chi-squared test for given probabilities
## 
## data:  c(35, 24, 17, 19)
## X-squared = 8.2, df = 3, p-value = 0.04205

Testujme teraz na základe rovnakých dát H0, že kandidáti majú podporu 1/3, 1/3, 1/6 a 1/6. Všimnite si štruktúru výstupu funkcie chisq.test.

p0 <- c(1/3,1/3,1/6,1/6)
freqs <- c(25,32,18,20)
vysl <- chisq.test(freqs, p=p0)
vysl
## 
##  Chi-squared test for given probabilities
## 
## data:  freqs
## X-squared = 2.8, df = 3, p-value = 0.4235
vysl$expected
## [1] 31.66667 31.66667 15.83333 15.83333
vysl$residuals
## [1] -1.18469776  0.05923489  0.54451008  1.04713477

Bonus: Buffonova ihla

Na cvičeniach z pravdepodobnosti ste počítali úlohu o Buffonovej ihle, ktorej zadanie bolo nasledovné:

V rovine sú zakreslené rovnobežné priamky s rozostupmi \(t\). Určte pravdepodobnosť toho, že ak na túto rovinu náhodne hodíme ihlu dĺžky \(l < t\), tak pretne niektorú z priamok. Hľadaná pravdepodobnosť je \(p=\frac{2l}{\pi t}\). Vyjadrením \(\pi\) dostaneme \(\pi=\frac{2l}{pt}\).

Keby sme teraz na sústavu priamok hodili n ihiel a pre každú zaznamenali, či pretla alebo nepretla niektorú z priamok, mohli by sme odhadnúť \(\pi\) ako pomer počtu priamok, ktoré preťali priamku ku všetkým ihlám. Samozrejme, takýto odhad bude tým presnejší, čím väčšie je n.

Nasledujúca funkcia odhadne touto metódou \(\pi\) pre fixné n, l a t (môžete skúšať meniť parametre a pozorovať, aké odhady dostávate).

estPi<- function(n, l=1, t=2) {
 m <- 0
 for (i in 1:n) {
   x <- runif(1)
   theta <- runif(1, min=0, max=pi/2)
   if (x < l/2 * sin(theta)) {
     m <- m +1
   }
 }
 return(2*l*n/(t*m))
}

Pozrime sa teraz pozrieť na to, ako sa zlepšuje odhad s rastúcim n. Skúsme použiť predchádzajúcu funkciu a nasimulovať priebeh odhadov, keď budeme brať n od 8000 do 120000:

n=8000
r=15
mat=rep(NA,r)
size=rep(NA,r)
for (i in 1:r) {
 size[i]<-n*i
 mat[i]<-estPi(n*i,l=1,t=2)
}
matrix<-expand.grid(size)
matrix[,2]<-mat
names(matrix)<-list("n","pi")
plot(matrix,type="b");abline(h=pi,col="red",lty=2)

Pekná animácia tohto experimentu je v balíku animation:

library(animation)
ani.options(nmax = 100, interval = 0.5)
par(mar = c(3, 2.5, 0.5, 0.2), pch = 20, mgp = c(1.5, 0.5, 0))
buffon.needle(mat = matrix(c(1, 2, 1, 3), 2))

Úlohy na precvičenie

  1. Máme n nezávislých javov, ktoré nastávajú s pravdepodobnosťami \((p_1,...,p_n)\). Napíšte funkciu, ktorej zadáme vektor \((p_1,...,p_n)\) a ktorá vypočíta pravdepodobnosť, že nastane práve jeden z týchto javov.

  2. Nech pravdepodobnosť prenosu choroby z komára na človeka je pri jednom uštipnutí 1/500. Napíšte funkciu, ktorá vypočíta pravdepodobnosť nákazy po uštipnutí n komármi.

  3. Vygenerujte 1000 bodov z trojrozmerného normálneho rozdelenia s parametrami \(\mu=(1,2,3)'\) a \[ \Sigma=\begin{bmatrix} 1 & 0 & 0\\ 0 & 2 & 0\\ 0 & 0 & 3 \end{bmatrix}. \] Vypočítajte výberový priemer a výberovú kovariančnú maticu takto vygenerovaných bodov.

  4. Vyberte si nejaký dátový súbor dostupný v R cez príkaz data() a vykreslite korelačnú maticu jednotlivých premenných. Skúste si vybrať niektoré premenné a zostaviť lineárny model.