NMST539 | Cvičenie 2

zdrojový Rmd súbor (kódovanie UTF8)

---

Outline druhého cvičenia:

• vlastne čísla v vlastné vektory matíc;
• náhodné vektory, podmienené a marginálne rozdelenia;
• mnohorozmerné normálne rozdelenie (R knižnica ‘mvtnorm’);

Cieľom druhého cvičenia z mnohorozmernej analýzy je zopakovať si niektoré dôležité a pre mnohorozmerné štatistické postupy podstatné znalosti týkajúce sa matíc a rozkladu matíc na vlastné vektory a vlastné čísla.

Stručne si pripomenieme niektoré základné vlastnosti náhodných vektorov. Uvedených bude niekoľko teoretických aj praktických príkladov na počítanie združených, marginálných a podmienených rozdeleni.

Štatisticky program R je k dispozícii (GNU public licence) na adrese https://www.r-project.org

RStudio (tzv. “user-friendly” interface - jeden z mnohých, ktoré na internete sú dostupné): RStudio.

Užitočné návody a manuály pre prácu s programom R:

• Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
• Komárek, A.: Základy práce s R. (PDF súbor)
• Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
• De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)

Odborná literatúra:

• Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 201
• Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.

---

Samostatne

---

• Pripomeňte a zopakujte si základné vlastnosti matíc, vlastných čísel a vlastných vektorov.
• Pokúste sa pomocou programu R naprogramovať vlastnú funkciu, ktorá pre jednoduchú štvorcovú maticu typu $2 \times 2$, alebo $3 \times 3$ nájde príslušné vlastné vektory a vlastné čísla.
• V programe R je pre spektrálny rozklad (t.j. nájdenie vlastných čísel a vlastných vektorov štvorcovej matice) definovaná funkcia eigen(). Pomocou helpu sa podívajte, ako je funkcia implementována a jej fungovanie si vyskúšajte na nejakých príkladoch.

---

Samostatne

---

• Pokúste sa premyslieť si, ako by ste pomocou vyššie uvedeného postupu vizuálne dokázali overili závislosť/nezávislosť dvoch náhodných veličín v nejakom náhodnom výbere $(X_i, Y_i)^\top$, pre $i = 1, \dots, n$.

• Uvedomte si, že klasický $xy$ scatterplot býva často hodne zavádzajúci v posudzovani závislosti a nezávislosti dvoch náhodných veličín. Hodne evidentné to je napríklad pre exponenciálne rozdelenie - viď ilustráciu nižsie:

  n <- 100
  x <- rexp(n, rate = 1)
  y <- rexp(n, rate = 10)
  z <- x + rnorm(n, 0, 2)
  
  par(mfrow = c(1,2))
  plot(y ~ x, pch = 21, bg = "lightblue")
  plot(z ~ x, pch = 21, bg = "lightblue")

  

  Zatiaľ čo prvý obrázok zobrazuje dve nezávisle zložky $X$ a $Y$ - napriek tomu, že dochádza k evidentnému zhlukovaniu bodov v počiatku, druhý scatterplot zobrazuje závisle zložky $X$ a $Z$ - ktoré ale vizuálne by sme asi skôr mali tendenciu posúdiť, ako nezávisle. Celkový obraz sa ale výrazne zmeni, ak namiesto scatterplotov pre pôvodné data $X_i$, $Y_i$ a $Z_i$ (pre $i = 1, \dots, n$) vykreslíme scatterplot pre $F_{X}(X_i)$, $F_{Y}(Y_i)$ a $F_{Z}(Z_i)$, kde $F_{X}$, $F_{Y}$ a $F_{Z}$ sú distribučné $X$, $Y$ a $Z$.

---

library("mvtnorm")

Poznámka

---

Pre umožnenie replikácie experimentu je vhodné používať nastavenie seedu pomocou príkazu set.seed(). Z vedeckého hľadiska je toho dokonca striktná požiadavka, aby bolo možné experiment, resp. štatisticku analýzu dodatočne zopakovať/rekonštruovať a získať ekvivalentné výsledky.

Porovnajte opakované použitie príkazu

(sample <- rnorm(10))

##  [1] -1.7685605 -0.7717708  0.4079737 -0.7056028 -0.7183402 -0.3604600
##  [7] -0.3311834  0.1055884 -0.7350938 -2.2298463

s opakovaným použitím dvojice príkazov

set.seed(1234)
(sample <- rnorm(10))

##  [1] -1.2070657  0.2774292  1.0844412 -2.3456977  0.4291247  0.5060559
##  [7] -0.5747400 -0.5466319 -0.5644520 -0.8900378

---

---

Príklad

Náhodný výber z dvojrozmerného normálneho rozdelenia

$\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_n \sim N_{2}\Big(\boldsymbol{\mu} = \Big(\begin{array}{c}2\\3\end{array}\Big), \Sigma = \Big(\begin{array}{cc}10^2 & 6^2\\ 6^2 & 6^2\end{array}\Big) \Big)$.

pre $\boldsymbol{X}_{i} = (X_{i 1}, X_{i 2})^\top$.

Pri volaní príkazu rmvnorm() je nutné špecifikovať požadovaný vektor stredných hodnôt a variančnú-kovariančnú maticu (t.j. v jednorozmernom prípade parameter rozptylu). Avšak v prípade jednorozmerného náhodného generátoru z normálneho rozdelenia (príkaz rnorm()) je nutné špecifikovať parameter smerodatnej chyby (na rozdiel od parameteru rozptylu) - porovnajte príslušný help k dvojici príkazov?rnorm a ?rmnorm. Je preto nutné dostatočne dobre poznať implementáciu používaných príkazov, aby nedošlo k omylom pri generováni.

• Aké je marginálne rozdelenie náhodných zložiek $X_{i 1}$ a $X_{i 2}$? Tie môžeme samozrejme nagenerovať (samostatne) použitím jednorozmerného generátoru - funkcie rnorm().
• Je možné využiť tento jednorozmerný náhodný generátor rnorm() k simulácii jednotlivých zložiek a následne vygenerované zložky vhodne zkombinovať aby sme dosiahli požadované dvojrozmerné normálne rozdelenie?
• Overte, že variančná-kovariančná matica $\Sigma$ je pozitívne definitná.
• Porozmýšlajte nad príkladom dvojrozmerného rozdelenia, kde obe zložky majú marginálne normálne rozdelenie (s vhodnými parametrami), ale združené rozdelenie nie je normálne. Dokážete simulovať z takého rozdelenia?

Nasledujúcim príkazom vygenerujem náhodný výver z daného dvojrozmerného normálneho rozdelenia s rozsahom $n = 1000$:

set.seed(1234)
s1 <- rmvnorm(1000, c(2, 3), matrix(c(10^2, 6^2, 6^2, 6^2),2,2))

Jednoduchý $xy$ scatterplot pomocou príkazu plot():

plot(s1, pch=21, xlim=c(-35, 35), ylim=c(-35,35), xlab="Marginal No. 1", ylab="Marginal No. 2", bg = "lightblue")

Pomocou výberových charakteristík - vektor výberových stredných hodnôt a výberová variančná-kovariančná matica - môžeme zobraziť hĺbkové kontúry daného rozdelenia (čo sa tým myslí?) Nakoľko je náhodný výber generovaný zo znamého teoretického rozdelenia (vyššie), je možné využiť aj teoretický predpis príslušnej dvojrozmernej hustoty a napočítať teoretické kontúry.

V nasledujúcom využijeme výberové charakteristiky spočítane z nagenerovaného náhodného výberu:

center <- apply(s1, 2, mean)   ## sample mean vector
sigma <- cov(s1)               ## sample variance-covariance matrix

Hustota dvojrozmerného normálneho rozdelenia je daná predpisom:

$f(\boldsymbol{x}) = f(x_{1}, x_{2}) = \frac{1}{2 \pi |\Sigma|^{1/2}} exp \left\{ -\frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) \right\}$,

kde $\boldsymbol{x} = (x_{1}, x_{2})^\top \in \mathbb{R}^2$, $\boldsymbol{\mu} = (\mu_{1}, \mu_{2})^\top$ je vektor stredných hodnôt a $\Sigma$ symetrická a pozitívne definitná variančná-kovariančná matica.

Pre spočítanie príslušných kontúr potrebujeme inverznú maticu k výberovej variančnej-kovariančnej matici $\widehat{\Sigma}$. Príslušné kontúry budeme získavať pomocou vlastnej funkcie ellipse():

sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

ellipse <- function(s,t){
                    u<-c(s,t)-center; 
                    e <- (-1) * (u %*% sigma.inv %*% u) / 2; 
                    exp(e)/(2 * pi * sqrt(det(sigma.inv)))}

V závislosti na zvolenom gride bodov spočítame niekoľko bodov na príslušných kontúrach pomocou funkcie ellipse():

n <- 60
x <- (0:(n-1))*2 - 50
y <- (0:(n-1))*2 - 50

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))

Výsledný graf:

plot(s1, pch=21, xlim=c(-35, 35), ylim=c(-35,35), xlab="Margin No. 1", ylab="Margin No. 2", bg = "lightgray", cex = 0.3)
contour(x,y,matrix(z,n,n), levels=(0:15), col = terrain.colors(16), add=TRUE, lwd = 1)

Samostatne

---

• Dokážete intuitivne vysvetliť význam jednotlivých parametrov v predpise dvojrozmernej hustoty normálneho rozdelenia v súvislosti s výsledným/očakávaným scatterplotom bodov generovaných z tohto rozdelnia?.




• Ako by ste intuitivne vysvetlili zmysel a úlohu hĺbkových kontúr? Napadá vás nejaká súvislosť s jednorozmernými kvantilmi?

---

library("mixtools")

rm(list = ls())
library(mixtools)
library(mvtnorm) 

set.seed(1234)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Margin No. 1", ylab="Margin No. 2")
ellipse(mu=colMeans(p), sigma=cov(p), alpha = .05, npoints = 250, col="red") 

Samostatne

---

Manuálna rekonštrukcia kontúry z predchádzajúceho obrázku:

alpha <- 0.05
angles <- seq(0, 2*pi, length = 200)
cpoints <- cbind(cos(angles), sin(angles))

q <- sqrt(qchisq(1 - alpha, df = 2))
mu <- colMeans(p)
sigma <- cov(p)
sigmaEig <- eigen(sigma)
sigmaSqr <- sigmaEig$vectors %*% diag(sqrt(sigmaEig$values)) %*% t(sigmaEig$vectors)
Epoints <-  q * (cpoints %*% sigmaSqr)
Epoints[,1] <- Epoints[,1] + mu[1]
Epoints[,2] <- Epoints[,2] + mu[2]

plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Margin No.1", ylab="Margin No.2")
points(Epoints[,1], Epoints[,2], col = "red", pch = "*")

---

library("MASS")

### two independent samples from two-dimensional normal each
bivn1 <- mvrnorm(1000, mu = c(0, 0), Sigma = matrix(c(1, .5, .5, 1), 2))
bivn2 <- mvrnorm(1000, mu = c(0, 0), Sigma = matrix(c(1.5, 1.5, 1.5, 1.5), 2))

### two dimensional kernel density estimate
bivn.kde1 <- kde2d(bivn1[,1], bivn1[,2], n = 50)
bivn.kde2 <- kde2d(bivn2[,1], bivn2[,2], n = 50)

par(mfrow = c(1,2))
persp(bivn.kde1, phi = 45, theta = 30, shade = .1, border = NA)
#par(new=TRUE)
persp(bivn.kde2, phi = 45, theta = 30, shade = .1, border = NA)

Poznámka

---

Je potrebné si dostatočne uvedomiť a pamätať na to, že zo združeného rozdelenia náhodného vektoru je možné (väčšínou) jednoducho a priamočiaro získať ľubovolné marginálne rozdelenie. Na druhú stranu, aj znalosť všetkých možných marginalov (okrem marginálneho rozdelenia, ktoré je vlastne už združeným rozdelenim) nestači k tomu, aby bolo možné určiť/spočítať celkové združené rozdelenie.

Vzhľadom k tomuto princípu môžeme pomocou jednorozmerného generátora generovať jednotlivé normálne rozdelené zložky, ale nedokážeme kontrolovať celkovú závislostnú štruktúru.

library(ggplot2)

rm(list = ls())
set.seed(1234)

n <- 1000
sample1 <- rnorm(n, mean=2, sd = 2)
sample2 <- 0.2 * sample1 + rnorm(n)
data <- data.frame(x = sample1, y = sample2, group="A")

sample3 <- rnorm(n, mean=0, sd = 1)
sample4 <- 3*sample3 - 1 + rnorm(n)
data <- rbind(data, data.frame(x = sample3, y = sample4, group="B"))

ggplot(data, aes(x=x, y=y,colour=group)) + geom_point(size=1.5, alpha=.6) 

Overenie závislosti/nezávislosti pomocou distribučných funkcii:

xPoints <- pnorm(sample1, mean(sample1), sd(sample1))
yPoints <- pnorm(sample2, mean(sample2), sd(sample2))

plot(xPoints, yPoints, pch = 21, bg = "lightblue")
lines(lowess(yPoints ~ xPoints), col = "red", lwd = 2)

Porovnajte s nasledujúcim obrázkom, pre sample1 a sample3 - teda dva nezávislé náhodné výbery.

xPoints <- pnorm(sample1, mean(sample1), sd(sample1))
yPoints <- pnorm(sample3, mean(sample3), sd(sample3))

plot(xPoints, yPoints, pch = 21, bg = "lightblue")
lines(lowess(yPoints ~ xPoints), col = "red", lwd = 2)

Ot8zky

---

• Korespondujú obrázky vyššie s tým, čo by ste v nich očakávali vzhľadom k tomu, čo už viete a čo viete o datach?
• Ako by ste získali analogicke obrázky, ak by teoretické rozdelenie náhodných výberov nebolo známe?
  (Hint: pre jednoduchosť uvažujte nejake jednorozmerné rozdelenie)




emp <- function(x, X){
   return(sum(X <= x)/length(X))
}

n <- 1000

X <- rnorm(n)
Y <- rnorm(n)
Z <- X + Y

x <- NULL
y <- NULL
z <- NULL

for (i in 1:length(X)){
  x <- c(x, emp(X[i], X))
  y <- c(y, emp(Y[i], Y))
  z <- c(z, emp(Z[i], Z))
  }

par(mfrow = c(1,2))
plot(ecdf(X))
points(x ~ X, col = "red", cex = 1.2)
plot(ecdf(Y))
points(y ~ Y, col = "red", cex = 1.2)



par(mfrow = c(1,2))
plot(x ~ y)
lines(lowess(x ~ y), col = "red", lwd = 2)
plot(z ~ x)
lines(lowess(z ~ x), col = "red", lwd = 2)

---

---