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.:
binomické: pbinom,qbinom,dbinom,rbinom
chí-kvadrát: pchisq,qchisq,dchisq,rchisq
exponenciálne: pexp,qexp,dexp,rexp
Fisherovo: pf,qf,df,rf
geometrické: pgeom,qgeom,dgeom,rgeom
hypergeometrické: phyper,qhyper,dhyper,rhyper
Poissonovo: ppois,qpois,dpois,rpois
Studentovo t: pt,qt,dt,rt
rovnomerné: punif,qunif,dunif,runif
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
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")
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
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.
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")
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)
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
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))
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.
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.
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.
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.