Základné operácie v R

Globálne nastavenie Chunkov

V nižšie uvedenom Chunku je urobené základné globálne nastavenie Chunkov v celom Notebooku.

• echo nastavuje, či chceme v Notebooku vypisovať jednotlivé kódy R
• message je nastavený na FALSE, čo znamená, že sa budú potláčať pracovné výstupy z R, t.j. napríklad výsledok otvárania knižníc, a niektoré iné výstupy, ktoré ale pre celkový Notebook nemajú takmer žiaden význam a sú vhodné skôr pri ladení kódov.
• warning nastavený na FALSE potláča zobrazenie chybových hlášok

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE
)

Úvod k základným operáciám v R

Tento notebook demonštruje základné operácie v jazyku R so:

• skalárnymi číslami (t.j. 1 číslo),
• textovými (znakovými) reťazcami,
• logickými (boolovskými) hodnotami a premennými,
• (numerickými) vektormi,
• maticami,
• vlastnými číslami a vlastnými vektormi (variančno-kovariančných) matíc.

Tam, kde je to užitočné, sú zahrnuté malé cvičenia.

---

Skaláre (jednočíselné hodnoty)

Numerické skaláre

# Priradenie konštanty do premennej
a <- 8
b <- 3.75

# Arithmetic
sum_ab      <- a + b        # sucet
diff_ab     <- a - b        # rozdiel
prod_ab     <- a * b        # násobenie
quot_ab     <- a / b        # delenie
power_ab    <- a ^ b        # umocňovanie
mod_ab      <- a %% 3       # zbytok po delení tromi (tzv modulo)
int_div_ab  <- a %/% b      # celočíselné delenie


# Rounding
round_b   <- round(b)       # zaokruhlovanie smerom k najblizsiemu celemu cislu
ceil_b    <- ceiling(b)     # najblizsie vyssie cele cislo
floor_b   <- floor(b)       # najblizsie nizsie cele cislo

a; b

[1] 8

[1] 3.75

sum_ab; diff_ab; prod_ab; quot_ab; power_ab; mod_ab; int_div_ab

[1] 11.75

[1] 4.25

[1] 30

[1] 2.133333

[1] 2435.496

[1] 2

[1] 2

round_b; ceil_b; floor_b

[1] 4

[1] 4

[1] 3

Poznámky

• ^ operátor umocňovania.
• %% je modulo, teda zbytok po delení,
• round(x, digits = 0) zaokrúhľovanie na určitý počet desatinných miest (digits=). ak digits = 0, potom ide o celočíselné zaokrúhľovanie

Malé cvičenie

Vypočítajte:

\[\frac{(15^2-4)}{7}\]

(15^2 - 4) / 7

[1] 31.57143

moje vypracovanie:

\[\frac{(18^2-6^2)\cdot3+\sqrt{144}}{(5^2-4)\cdot2}\]

((18^2 - 6^2) * 3 + sqrt(144)) / ((5^2 - 4) * 2)

[1] 20.85714

---

Text

Vytváranie textovými premennými a práca s nimi

first <- "Ema "                   # definovanie obsahu textovej premennej first
last  <-  "Šimkovicová"  # definovanie obsahu text. premennej last
full  <- paste(first, last)               # spojenie dvoch text. premennych do jednej (s medzerou)
full_nospace <- paste0(first, last)       # spojenie bez medzery
csv_line <- paste("apple", "banana", "pear", sep = ",")  # spojenie textov s oddelovacom ,
first; last; full; full_nospace; csv_line   # bodkočiarka tu nahradzuje odskok na novy riadok 

[1] "Ema "

[1] "Šimkovicová"

[1] "Ema  Šimkovicová"

[1] "Ema Šimkovicová"

[1] "apple,banana,pear"

Dĺžka textového reťazca, podreťazec

x <- "R is great!"
nchar(x)                 # počet znakov  v retazci "R is great!"

[1] 11

substr(x, 1, 5)          # podreťazec od 1. do 5. znaku

[1] "R is "

Tip: Knižnica stringr mnohé zaujímavé možnosti práce s textami, ale implicitné knižnice R pokrývajú väčšinu bežných potrieb páce s textami.

---

Logické (boolovské) hodnoty a premenné

Základy

p <- TRUE
q <- FALSE
!p                 # NOT

[1] FALSE

p & q              # AND

[1] FALSE

p | q              # OR

[1] TRUE

xor(p, q)          # exclusive OR - platí len jedno z p,alebo q

[1] TRUE

Logický výsledok porovnávania

3 < 5

[1] TRUE

7 >= 7

[1] TRUE

"cat" == "cat"

[1] TRUE

"cat" != "dog"   # vykricnik je tu v zmysle negacie. Napr.:   !=, !>, !<, !TRUE

[1] TRUE

!TRUE

[1] FALSE

Zložitejšie logické operácie

x <- 10
x > 5 & x < 20      # a sucasne - logicky prienik (sucin)

[1] TRUE

x < 0 | x > 100     # alebo - logicke zjednotenie (sucet)

[1] FALSE

                    # pri zlozitejsich vztahoch pouzivajte zatvorky ()

Zlučovanie viacerých log. premenných do vektora

vals <- c(TRUE, FALSE, TRUE, TRUE)   # definicia vektora s logickymi hodnotami

---

Numerické vektory

Generovanie vektorov

v1 <- c(2, 4, 6, 8)
v2 <- 1:5                  # postupnost 1,2,3,4,5
v3 <- seq(from = 0, to = 1, by = 0.25)  # postupnost s krokom 0.25
v4 <- rep(3, times = 5)    # 3,3,3,3,3  # 5 clenna postupnost trojak
v5 <- runif(5)             # generovanie rovnomerne rozdelenych premennych v intervale [0,1]
v6 <- rnorm(5)             # generovanie normalne rozdelenych premennych
v1; v2; v3; v4; v5

[1] 2 4 6 8

[1] 1 2 3 4 5

[1] 0.00 0.25 0.50 0.75 1.00

[1] 3 3 3 3 3

[1] 0.1655344 0.1104430 0.9601715 0.8733302 0.1476671

Aritmetické operácie s vektormi

v <- c(1, 2, 3, 4)
v + 10           # kazdy prvok vektora zvacsime o 10

[1] 11 12 13 14

v * 2            # kazdy prvok vektora prenasobime 2

[1] 2 4 6 8

(v + 1) / 2

[1] 1.0 1.5 2.0 2.5

exp(v)           # exponencialna funkcia z kazdeho prvku vektora

[1]  2.718282  7.389056 20.085537 54.598150

sum(c(1,2,3),c(1,1,1))          # skalarny sucin - vysledok je skalar

[1] 9

crossprod(c(1,2,3),c(1,1,1))    # skalarny sucin - vysledok je matica 1x1

     [,1]
[1,]    6

c(1,2,3)*c(1,1,1)               # Hadamardov sucin (sucin zodpovedajucich prvkov vektora)

[1] 1 2 3

Matematické operácie s 2 vektormi rovnakého rozmeru

length(c(1,2,3,4,5))

[1] 5

length(v5)            #vektor v5 je definovany vyssie

[1] 5

c(1,2,3,4,5) + v5     # pozor, oba vektory musia mat rovnaky rozmer

[1] 1.165534 2.110443 3.960172 4.873330 5.147667

Indexovanie a výber niektorych prvkov vektora

x <- c(5, 12, 3, 18, 7, 0, 21)
x[1]           # indexovanie - novy jedno-prvkovy vektor - prvy prvok vektora x

[1] 5

x[2:4]         # novy vektor s druhym az stvrtym prvkom vektora x

[1] 12  3 18

x[-1]          # novy vektor - vsetky prvky vektora x okrem prvého

[1] 12  3 18  7  0 21

x[x > 10]      # novy vektor definovany prvkami x vacsimi ako 10

[1] 12 18 21

which(x > 10)  # ktore prvky zodpovedaju podmienke vacsieho ako 10?

[1] 2 4 7

Práca s chýbajúcimi hodnotami

y <- c(1, NA, 3, NA, 5)
is.na(y)

[1] FALSE  TRUE FALSE  TRUE FALSE

mean(y)                 # NA

[1] NA

mean(y, na.rm = TRUE)   # remove NAs

[1] 3

Základné štatistiky a usporiadanie prvkov vektora podľa veľkosti

z <- c(10, 3, 5, 8, 2)
mean(z)                 # priemerna hodnota

[1] 5.6

sd(z)                  # standardna odchylka

[1] 3.361547

max(z)                  # maximalna hodnota

[1] 10

summary(z)              # rychly prehlad zakladnych statistik o vektore

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    2.0     3.0     5.0     5.6     8.0    10.0 

sort(z)                 # rastuce usporiadanie 

[1]  2  3  5  8 10

sort(z, decreasing = TRUE)  # klesajuce

[1] 10  8  5  3  2

Malé cvičenie

Vytvorte vektor w s číslami 1..20 a vypočítajte sumu všetkých párnych čísel.

w <- 1:20
sum(w[w %% 2 == 0])

[1] 110

Moje vypracovanie: Vytvorte vektor w s číslami 1..150 vypočítajte sumu všetkých nepárnych čísel, ktoré sú zároveň deliteľné číslami 3 a 5.

w <- 1:150
sum(w[w %% 2 != 0 & w %% 3 == 0 & w %% 5 == 0])

[1] 375

---

Matice

Vytvorenie matíc

m <- matrix(1:12, nrow = 3, ncol = 4)            # hodnoty sú zadavane po stlpcoch
m_byrow <- matrix(1:12, nrow = 3, byrow = TRUE)  # hodnoty su zadavane po riadkoch
m; m_byrow

     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12

Rozmery matice

dim(m)                   # (rows, cols)

[1] 3 4

m

     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

Adresovanie prvkov matice

m[1, 2]      # riadok 1, stlpec 2

[1] 4

m[ , 3]      # vsetky prvky v tretom stlpci - vysledok matica 3x1

[1] 7 8 9

m[2, ]       # vsetky prvky v druhom riadku - vysledok matica 1*3

[1]  2  5  8 11

m[1:2, 2:3]  # podmatica tvorena riadkami 1, 2 a stlpcami 2, 3

     [,1] [,2]
[1,]    4    7
[2,]    5    8

Maticové operácie

A <- matrix(c(1,2,3,4), nrow = 2)
B <- matrix(c(5,6,7,8), nrow = 2)

A + B        # scitanie matic

     [,1] [,2]
[1,]    6   10
[2,]    8   12

A * B        # Hadamard product - nasobenie po zodpovedajucich prvkoch

     [,1] [,2]
[1,]    5   21
[2,]   12   32

A %*% B      # nasobenie matic

     [,1] [,2]
[1,]   23   31
[2,]   34   46

t(A)         # transpozicia matice A - vymena riadkov a stlpcov

     [,1] [,2]
[1,]    1    2
[2,]    3    4

det(A)       # determinant matice

[1] -2

solve(A)     # inverzia matice (ak je matica regularna - teda inverzia sa da spocitat)

     [,1] [,2]
[1,]   -2  1.5
[2,]    1 -0.5

Zlučovanie vektorov do matíc

C <- cbind(1:3, 4:6)   #  - po stlpcoch 
D <- rbind(1:3, 4:6)   #  - po riadkoch 
C; D

     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6

Vypočítanie zvolenej štatistiky po riadkoch (stĺpcoch) matice

M <- matrix(1:9, nrow = 3)
M

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

apply(M, 1, sum)   # suma po riadkoch

[1] 12 15 18

apply(M, 2, mean)  # priemery po stĺpcoch

[1] 2 5 8

Malé cvičenie

Vytvorte maticu 5x5 s hodnotami po riadkoch 1..25, vypočítajte stĺpcové sumy a súčin matíc (M^t M).

M2 <- matrix(1:25, nrow = 5, byrow = TRUE)
colSums(M2)

[1] 55 60 65 70 75

t(M2) %*% M2

     [,1] [,2] [,3] [,4] [,5]
[1,]  855  910  965 1020 1075
[2,]  910  970 1030 1090 1150
[3,]  965 1030 1095 1160 1225
[4,] 1020 1090 1160 1230 1300
[5,] 1075 1150 1225 1300 1375

Moje vypracovanie: Vytvorte maticu 8x8 s hodnotami po riadkoch 1..64, vypočítajte riadkové a stĺpcové sumy, transponovanú maticu a súčin matíc (M^t M).

M2 <- matrix(1:64, nrow = 8, byrow = TRUE)
rowSums(M2)

[1]  36 100 164 228 292 356 420 484

colSums(M2)

[1] 232 240 248 256 264 272 280 288

t(M2)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    9   17   25   33   41   49   57
[2,]    2   10   18   26   34   42   50   58
[3,]    3   11   19   27   35   43   51   59
[4,]    4   12   20   28   36   44   52   60
[5,]    5   13   21   29   37   45   53   61
[6,]    6   14   22   30   38   46   54   62
[7,]    7   15   23   31   39   47   55   63
[8,]    8   16   24   32   40   48   56   64

t(M2) %*% M2

      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,]  9416  9648  9880 10112 10344 10576 10808 11040
[2,]  9648  9888 10128 10368 10608 10848 11088 11328
[3,]  9880 10128 10376 10624 10872 11120 11368 11616
[4,] 10112 10368 10624 10880 11136 11392 11648 11904
[5,] 10344 10608 10872 11136 11400 11664 11928 12192
[6,] 10576 10848 11120 11392 11664 11936 12208 12480
[7,] 10808 11088 11368 11648 11928 12208 12488 12768
[8,] 11040 11328 11616 11904 12192 12480 12768 13056

---

Vlastné čísla a vlastné vektory matíc

Čo znamenajú vlastné čísla a vlastné vektory?

Nech \(\mathbf A\) je štvorcová matica. Vlastné číslo (eigenvalue) \(\lambda\) a príslušný vlastný vektor (eigenvector) \(\mathbf v \ne \mathbf 0\) spĺňajú vzťah

\[ \mathbf A \mathbf v = \lambda \mathbf v. \]

Násobenie maticou \(\mathbf A\) teda pôsobí na vlastný vektor rovnako ako násobenie číslom \(\lambda\). Výsledok reprezentuje vektor rovnakého smeru, ale rozličných veľkostí:

• ak \(\lambda > 0\), orientácia vektora sa zachová;
• ak \(\lambda < 0\), orientácia sa obráti;
• ak \(\lambda = 0\), výsledkom je nulový vektor.

Dĺžka sa násobí hodnotou \(|\lambda|\). Napríklad pri \(\lambda = 3\) sa vektor predĺži trojnásobne, pri \(\lambda = 0.5\) sa skráti na polovicu. Tento geometrický opis sa týka reálnych vlastných čísel a vektorov.

Jednoduchý príklad: matica 2 × 2

Použijeme symetrickú maticu, teda maticu, pre ktorú platí $A^T

\[ \mathbf A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}. \]

Vlastné čísla nájdeme podľa nasledovného postupu:

\[\begin{eqnarray} \mathbf A \mathbf h & = & \lambda \mathbf h \\ \mathbf A \mathbf h - \lambda \mathbf h & = & \mathbf 0 \\ (\mathbf A - \lambda \mathbf I)\mathbf h & = & \mathbf 0, \end{eqnarray} \]

kde \(\mathbf I\) je jednotková matica.

Netriviálne (nenulové) riešenie tejto sústavy existuje vtedy, keď

\[ \det(\mathbf A-\lambda \mathbf I)=0. \]

V našom prípade

\[ \mathbf A-\lambda \mathbf I = \begin{pmatrix} 2-\lambda & 1 \\ 1 & 2-\lambda \end{pmatrix}, \]

a preto

\[ \det(\mathbf A-\lambda \mathbf I) = \begin{vmatrix} 2-\lambda & 1 \\ 1 & 2-\lambda \end{vmatrix} = (2-\lambda)^2-1. \]

Dostávame charakteristickú rovnicu

\[ (2-\lambda)^2-1 = \lambda^2-4\lambda+3 = (\lambda-3)(\lambda-1) = 0. \]

Vlastné čísla sú koreňmi (riešeniami) kvadratickej rovnice a teda

\[ \lambda_1=3, \qquad \lambda_2=1. \]

Výpočet vlastného vektora pre \(\lambda_1=3\)

Vlastný vektor vypočítame zo vzťahu

\[ (\mathbf A-\lambda_1\mathbf I)\mathbf h_1=\mathbf 0. \]

Po dosadení \(\lambda_1=3\) dostávame

\[ \left[ \begin{pmatrix} 2 & 1\\ 1 & 2 \end{pmatrix} - 3 \begin{pmatrix} 1&0\\ 0&1 \end{pmatrix} \right] \begin{pmatrix} h_{11}\\ h_{21} \end{pmatrix} = \begin{pmatrix} 0\\ 0 \end{pmatrix}. \]

Teda

\[ \begin{pmatrix} -1&1\\ 1&-1 \end{pmatrix} \begin{pmatrix} h_{11}\\ h_{21} \end{pmatrix} = \begin{pmatrix} 0\\ 0 \end{pmatrix}. \]

Získame sústavu rovníc

\[ \begin{eqnarray} -h_{11}+h_{21} &=&0,\\ h_{11}-h_{21} &=&0, \end{eqnarray} \] pričom zavedieme ešte dodatočnú podmienku normalizácie vlastných vektorov, teda

\[\sqrt{h_{11}^2 + h_{21}^2} = 1\]

Prvé dve rovnice vyjadrujú tú istú podmienku:

\[ h_{11}=h_{21}. \]

Teda z normalizačnej podmienky dostávame

\[\sqrt{h_{11}^2 + h_{21}^2} = 1 \Longleftrightarrow h_{11} = h_{21} = \pm \sqrt{\frac{1}{2}}\]

Môžeme preto napríklad zvoliť

\[ h_{11}=+\sqrt{\frac{1}{2}}, \]

a potom

\[ h_{21}=+\sqrt{\frac{1}{2}}. \]

Vlastný vektor zodpovedajúci vlastnému číslu \(\lambda_1=3\) môže byť teda

\[ \mathbf h_1= \begin{pmatrix} +\sqrt{\frac{1}{2}}\\ +\sqrt{\frac{1}{2}} \end{pmatrix}. \]

Každý nenulový násobok tohto vektora je však tiež vlastným vektorom. Napríklad

\[ \begin{pmatrix} +1\\ +1 \end{pmatrix} \]

predstavuje ten istý smer.

Výpočet vlastného vektora pre \(\lambda_2 = 1\)

Keďže platí

\[ \begin{aligned} \left[ \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} - 1 \cdot \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \right] \begin{pmatrix} h_{12} \\ h_{22} \end{pmatrix} &= \begin{pmatrix} 0 \\ 0 \end{pmatrix} \\[1em] \Longleftrightarrow \qquad \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} h_{12} \\ h_{22} \end{pmatrix} &= \begin{pmatrix} 0 \\ 0 \end{pmatrix}. \end{aligned} \]

dostávame sústavu rovníc

\[ \begin{eqnarray} h_{12}+h_{22} &=&0,\\ h_{12}+h_{22} &=&0. \end{eqnarray} \]

Z toho vyplýva

\[ h_{22}=-h_{12}. \]

Ak zvolíme dodatočnú normalizačnú podmienku

\[\sqrt{h_{12}^2 + h_{22}^2} = 1,\] potom ak

\[ h_{12}=\sqrt{\frac{1}{2}}, \]

pre druhý prvok vektora bude platiť

\[ h_{22}=-\sqrt{\frac{1}{2}}. \]

Vlastný vektor zodpovedajúci vlastnému číslu \(\lambda_2=1\) teda môžeme zapísať ako

\[ \mathbf h_2= \begin{pmatrix} \sqrt{\frac{1}{2}}\\ -\sqrt{\frac{1}{2}} \end{pmatrix}. \]

Opäť platí, že ľubovoľný nenulový násobok tohto vektora je tiež vlastným vektorom.

---

Overenie pre prvý vlastný vektor

Pre prvý vlastný vektor platí

\[ \mathbf A\mathbf h_1 = \begin{pmatrix} 2&1\\ 1&2 \end{pmatrix} \begin{pmatrix} \sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{2}} \end{pmatrix} = 3 \begin{pmatrix} \sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{2}} \end{pmatrix} = \lambda_1 \mathbf h_1 \]

Teda

\[ \mathbf A\mathbf h_1=\lambda_1\mathbf h_1. \]

Overenie pre druhý vlastný vektor

Podobne pre druhé vlastné číslo \(\lambda_2 = 1\) a druhý vlastný vektor \(\mathbf h_2\) platí

\[ \mathbf A\mathbf h_2 = \begin{pmatrix} 2&1\\ 1&2 \end{pmatrix} \begin{pmatrix} \sqrt{\frac{1}{2}}\\ -\sqrt{\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} 2\sqrt{\frac{1}{2}}-\sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{2}}-2\sqrt{\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} \sqrt{\frac{1}{2}}\\ -\sqrt{\frac{1}{2}} \end{pmatrix} = 1\cdot\mathbf h_2 = \lambda_2 \mathbf h_2. \]

Násobenie maticou \(\mathbf A\) teda vektor \(\mathbf h_2\) vôbec nezmení – jeho smer aj dĺžka zostanú rovnaké, pretože \(\lambda_2 = 1\).

---

Matice vlastných čísel a vlastných vektorov

Dostávame teda dvojice vlastných čísel a vlastných vektorov:

\[ \lambda_1=3, \qquad \mathbf h_1= \begin{pmatrix} \sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{2}} \end{pmatrix}, \]

a

\[ \lambda_2=1, \qquad \mathbf h_2= \begin{pmatrix} \sqrt{\frac{1}{2}}\\ -\sqrt{\frac{1}{2}} \end{pmatrix}. \]

Keďže matica \(\mathbf A\) je symetrická, vlastné vektory zodpovedajúce rôznym vlastným číslam sú navzájom ortogonálne (navzájom kolmé a po dvojiciach ich skalárne súčiny sa rovnajú nule). Skutočne,

\[ \mathbf h_1^T\mathbf h_2 = \left( \sqrt{\frac{1}{2}},\sqrt{\frac{1}{2}} \right) \begin{pmatrix} \sqrt{\frac{1}{2}}\\ -\sqrt{\frac{1}{2}} \end{pmatrix} = \frac{1}{2} - \frac{1}{2} = 0. \]

Vlastné vektory teraz zapíšeme do stĺpcov matice \(\mathbf H\) a vlastné čísla na diagonálu matice \(\boldsymbol\Lambda\):

\[ \mathbf H = \left(\mathbf h_1, \mathbf h_2\right) = \begin{pmatrix} \sqrt{\frac{1}{2}} & \sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{2}} & -\sqrt{\frac{1}{2}} \end{pmatrix}, \qquad \boldsymbol\Lambda = \begin{pmatrix} 3 & 0\\ 0 & 1 \end{pmatrix}. \]

Obe rovnice \(\mathbf A\mathbf h_1=\lambda_1\mathbf h_1\) a \(\mathbf A\mathbf h_2=\lambda_2\mathbf h_2\) potom môžeme zapísať naraz jedinou maticovou rovnicou

\[ \mathbf A\mathbf H = \mathbf H\boldsymbol\Lambda . \]

Keďže vlastné vektory sú normalizované a navzájom ortogonálne, platí \(\mathbf H^T\mathbf H = \mathbf I\), teda \(\mathbf H^{-1} = \mathbf H^T\). Ak rovnicu \(\mathbf A\mathbf H = \mathbf H\boldsymbol\Lambda\) vynásobíme sprava maticou \(\mathbf H^T\), dostaneme spektrálny rozklad matice \(\mathbf A\):

\[ \mathbf A = \mathbf H\boldsymbol\Lambda\mathbf H^T . \]

Overenie pre náš príklad:

\[ \mathbf H\boldsymbol\Lambda\mathbf H^T = \begin{pmatrix} 3\sqrt{\frac{1}{2}} & \sqrt{\frac{1}{2}}\\ 3\sqrt{\frac{1}{2}} & -\sqrt{\frac{1}{2}} \end{pmatrix} \begin{pmatrix} \sqrt{\frac{1}{2}} & \sqrt{\frac{1}{2}}\\ \sqrt{\frac{1}{2}} & -\sqrt{\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} \frac{3}{2}+\frac{1}{2} & \frac{3}{2}-\frac{1}{2}\\ \frac{3}{2}-\frac{1}{2} & \frac{3}{2}+\frac{1}{2} \end{pmatrix} = \begin{pmatrix} 2 & 1\\ 1 & 2 \end{pmatrix} = \mathbf A . \]

Spektrálny rozklad budeme neskôr používať pri variančno-kovariančných maticiach: matica sa „rozloží“ na smery (vlastné vektory) a na to, aká veľká časť variancie pripadá na každý smer (vlastné čísla).

---

Výpočet pomocou R funkcie eigen()

Funkcia eigen() je súčasťou základného R; nepotrebujeme ďalšiu knižnicu. Argument symmetric = TRUE použijeme, keď vieme, že matica je symetrická. Pri všeobecnej matici môžeme napísať jednoducho eigen() a R symetriu skontroluje samo.

# Vytvorenie matice z analytického príkladu
A_eig <- matrix(c(2, 1,
                  1, 2),
                nrow = 2, byrow = TRUE)

A_eig

     [,1] [,2]
[1,]    2    1
[2,]    1    2

# Výpočet vlastných čísel a vlastných vektorov
eig_A <- eigen(A_eig, symmetric = TRUE) #vyuzitie funkcie eigen()
eig_A                                   # vysledok sa ulozil do listu, ktory sme nazvali eig_A

eigen() decomposition
$values
[1] 3 1

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

 #Dalej pristupime k rozbaleniu obsahu listu eig_A

lambda_A <- eig_A$values    # vektor vlastných čísel - ulozenie do lambda_A
V_A <- eig_A$vectors        # vlastné vektory v stĺpcoch ulozenie do V_A

lambda_A

[1] 3 1

V_A

          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

Výsledok eig_A je list. Znak $ vyberá jeho pomenovanú zložku:

Zápis v R	Význam
eig_A$values	Všetky vlastné čísla
eig_A$vectors	Vlastné vektory uložené v stĺpcoch
eig_A$values[1]	Prvé vlastné číslo
eig_A$vectors[, 1]	Vlastný vektor prislúchajúci prvému vlastnému číslu

Pri reálnej symetrickej matici R vráti vlastné čísla zostupne; tu teda 3 a 1. Vlastné vektory majú jednotkovú dĺžku. Môžu sa preto líšiť od našich ručne zvolených vektorov mierkou alebo znamienkom. Vektory \(v\) aj \(-v\) sú rovnako správne: vlastný vektor nie je jednoznačný.

Ak potrebujeme iba vlastné čísla:

eigen(A_eig, symmetric = TRUE, only.values = TRUE)$values

[1] 3 1

Overenie vzťahu Av = λv

V R použijeme na násobenie matice a vektora operátor %*%. Operátor * používame pri násobení vektora skalárom.

lambda_1 <- lambda_A[1]
v_1 <- V_A[, 1]

lava_strana <- as.vector(A_eig %*% v_1)
prava_strana <- lambda_1 * v_1

cat("Ľavá strana je ", paste(lava_strana, collapse = ", "), 
    " a pravá strana je ", paste(prava_strana, collapse = ", "), "!", sep = "")

Ľavá strana je 2.12132034355964, 2.12132034355964 a pravá strana je 2.12132034355964, 2.12132034355964!

Všetky vlastné páry naraz overíme v maticovom tvare \(\mathbf A\mathbf H = \mathbf H\boldsymbol\Lambda\). Overíme aj ortonormalitu \(\mathbf H^T\mathbf H = \mathbf I\) a spektrálny rozklad \(\mathbf A = \mathbf H\boldsymbol\Lambda\mathbf H^T\). Pri výpočtoch na počítači vznikajú drobné zaokrúhľovacie chyby (rádovo \(10^{-16}\)), preto výsledky zaokrúhlime.

H_A <- V_A                       # vlastne vektory v stlpcoch
Lambda_A <- diag(lambda_A)       # diagonalna matica vlastnych cisel

round(A_eig %*% H_A - H_A %*% Lambda_A, 12)   # nulova matica:  A H = H Lambda

     [,1] [,2]
[1,]    0    0
[2,]    0    0

round(t(H_A) %*% H_A, 12)                     # jednotkova matica:  H^T H = I

     [,1] [,2]
[1,]    1    0
[2,]    0    1

H_A %*% Lambda_A %*% t(H_A)                   # spat dostaneme maticu A

     [,1] [,2]
[1,]    2    1
[2,]    1    2

Geometrické znázornenie

Porovnáme dva vlastné vektory. Pri vlastnom vektore zostáva jeho obraz (transformácia) na tej istej priamke.

local({
  povodne <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
  on.exit(par(povodne))

  priklady <- list(c(1, 1), c(1, -1))
  nazvy <- c("Vlastný vektor: lambda = 3",
             "Vlastný vektor: lambda = 1")

  for (j in seq_along(priklady)) {
    x_plot <- priklady[[j]]
    Ax_plot <- as.vector(A_eig %*% x_plot)
    plot(NA, xlim = c(-1.5, 3.5), ylim = c(-1.5, 3.5),
         asp = 1, xlab = "Prvá súradnica", ylab = "Druhá súradnica",
         main = nazvy[j], cex.main = 0.85)
    abline(h = 0, v = 0, col = "grey80")
    arrows(0, 0, Ax_plot[1], Ax_plot[2],
           col = "darkorange2", lwd = 4, length = 0.12)
    arrows(0, 0, x_plot[1], x_plot[2],
           col = "steelblue4", lwd = 2, lty = 2, length = 0.10)
    legend("topleft", legend = c("x", "Ax"),
           col = c("steelblue4", "darkorange2"),
           lty = c(2, 1), lwd = c(2, 4), bty = "n", cex = 0.8)
  }
})

Vektor x (modrý) a jeho obraz Ax (oranžový).

V pravom obrázku sa šípky prekrývajú, pretože \(Ax = x\). Cyklus for opakuje rovnaký postup pre oba vektory. Funkcia local() udržiava pomocné premenné grafu oddelene od ostatných výpočtov.

Matica 3 × 3 - príklad v R

B_eig <- matrix(c(2, 1, 0,
                  1, 2, 0,
                  0, 0, 4), nrow = 3, byrow = TRUE)
eig_B <- eigen(B_eig, symmetric = TRUE)
eig_B$values                # 4, 3, 1

[1] 4 3 1

eig_B$vectors

     [,1]      [,2]       [,3]
[1,]    0 0.7071068  0.7071068
[2,]    0 0.7071068 -0.7071068
[3,]    1 0.0000000  0.0000000

# Súčty štvorcov prvkov jednotlivých stĺpcov sú 1.
colSums(eig_B$vectors^2)

[1] 1 1 1

# Overenie všetkých troch vlastných párov.
max(abs(B_eig %*% eig_B$vectors -
        eig_B$vectors %*% diag(eig_B$values)))

[1] 5.551115e-16

Pre štvorcovú maticu platí: súčet vlastných čísel sa rovná stope matice (súčtu diagonálnych prvkov) a ich súčin sa rovná determinantu. Vlastné čísla pritom počítame aj s ich násobnosťami.

c(sucet_vlastnych_cisel = sum(eig_B$values),
  stopa_matice = sum(diag(B_eig)))             # obe hodnoty sú 8

sucet_vlastnych_cisel          stopa_matice 
                    8                     8 

c(sucin_vlastnych_cisel = prod(eig_B$values),
  determinant = det(B_eig))                   # obe hodnoty sú 12

sucin_vlastnych_cisel           determinant 
                   12                    12 

Malé cvičenie

Pre variančno kovariančnú maticu počítanú z daily close cien akcií za posledný rok, vypočítajte vlastné čísla a vlastné vektory.

Asset	T	JNK	TQQQ
T	0.000134	0.000008	0.000118
JNK	0.000008	0.000016	0.000062
TQQQ	0.000118	0.000062	0.001185

# vytvorenie matice
C <- matrix(c(
  0.000134, 0.000008, 0.000118,  # riadok 1 (T)
  0.000008, 0.000016, 0.000062,  # riadok 2 (JNK)
  0.000118, 0.000062, 0.001185   # riadok 3 (TQQQ)
), nrow = 3, byrow = TRUE)

Vypočítajte vlastné čísla a vlastné vektory a vypíšte ich

# vypocet vlastnych cisel a vlastnych vektorov
eig_C <- eigen(C, symmetric = TRUE)

# vypis vlastnych cisel
lambda_C <- eig_C$values
lambda_C

[1] 1.201381e-03 1.209261e-04 1.269271e-05

# vypis vlastnych vektorov
V_C <- eig_C$vectors
V_C

            [,1]        [,2]        [,3]
[1,] -0.11011897  0.99378875  0.01605395
[2,] -0.05265593  0.01029636 -0.99855963
[3,] -0.99252263 -0.11080569  0.05119505

Moje vypracovanie

Pre variančno kovariančnú maticu počítanú z daily close cien akcií za posledný rok, vypočítajte vlastné čísla a vlastné vektory.

Asset	AAPL	MSFT	NVDA	TSLA
AAPL	0.000245	0.000085	0.000145	0.000095
MSFT	0.000085	0.000195	0.000115	0.000075
NVDA	0.000145	0.000115	0.000685	0.000225
TSLA	0.000095	0.000075	0.000225	0.000895

Vypočítajte vlastné čísla a vlastné vektory a vypíšte ich

# Vytvorenie variancno-kovariancnej matice
C <- matrix(c(
  0.000245, 0.000085, 0.000145, 0.000095,
  0.000085, 0.000195, 0.000115, 0.000075,
  0.000145, 0.000115, 0.000685, 0.000225,
  0.000095, 0.000075, 0.000225, 0.000895
), nrow = 4, ncol = 4, byrow = TRUE)

# Pomenovanie riadkov a stlpcov
akcie <- c("AAPL", "MSFT", "NVDA", "TSLA")
dimnames(C) <- list(akcie, akcie)

# Vypis matice
print(C)

         AAPL     MSFT     NVDA     TSLA
AAPL 0.000245 0.000085 0.000145 0.000095
MSFT 0.000085 0.000195 0.000115 0.000075
NVDA 0.000145 0.000115 0.000685 0.000225
TSLA 0.000095 0.000075 0.000225 0.000895

# Kontrola symetrie matice
print(isSymmetric(C))

[1] TRUE

# Vypocet vlastnych cisel a vlastnych vektorov
eig_C <- eigen(C, symmetric = TRUE)

# Vypis vlastnych cisel
lambda_C <- eig_C$values
names(lambda_C) <- paste0("lambda_", 1:4)

cat("Vlastne cisla:\n")

Vlastne cisla:

print(lambda_C)

    lambda_1     lambda_2     lambda_3     lambda_4 
0.0010911918 0.0005658695 0.0002316186 0.0001313201 

# Vypis vlastnych vektorov
V_C <- eig_C$vectors
rownames(V_C) <- akcie
colnames(V_C) <- paste0("v", 1:4)

cat("Vlastne vektory:\n")

Vlastne vektory:

print(round(V_C, 6))

            v1        v2        v3        v4
AAPL -0.199718 -0.201305  0.759024  0.586065
MSFT -0.156430 -0.156404  0.542895 -0.810144
NVDA -0.554578 -0.752024 -0.355969  0.013724
TSLA -0.792518  0.607843 -0.049340  0.002615

# Graf vlastnych vektorov
matplot(
  x = 1:4,
  y = V_C,
  type = "b",
  pch = 19,
  lty = 1,
  col = c("blue", "red", "green", "purple"),
  xlab = "Akcie",
  ylab = "Hodnota vlastneho vektora",
  main = "Graf vlastnych vektorov",
  xaxt = "n",
  ylim = c(-1, 1.6)      # miesto nad krivkami, aby legenda nezakryvala body
)
abline(h = 0, col = "grey80")

# Oznacenie akcii na osi X
axis(
  side = 1,
  at = 1:4,
  labels = akcie
)

# Legenda
legend(
  "topright",
  legend = c(
    "Prvy vlastny vektor",
    "Druhy vlastny vektor",
    "Treti vlastny vektor",
    "Stvrty vlastny vektor"
  ),
  col = c("blue", "red", "green", "purple"),
  lty = 1,
  pch = 19,
  cex = 0.8
)

Výsledky ešte overíme a zistíme, akú časť celkovej variancie vysvetľuje každé vlastné číslo.

# Overenie C V = V Lambda pre vsetky styri vlastne pary naraz
max(abs(C %*% V_C - V_C %*% diag(lambda_C)))

[1] 4.336809e-19

# Sucet vlastnych cisel = stopa matice = sucet rozptylov styroch akcii
c(sucet_vlastnych_cisel = sum(lambda_C), stopa_matice = sum(diag(C)))

sucet_vlastnych_cisel          stopa_matice 
              0.00202               0.00202 

# Podiel celkovej variancie pripadajuci na jednotlive vlastne cisla (v %)
podiel_C <- 100 * lambda_C / sum(lambda_C)
round(podiel_C, 1)

lambda_1 lambda_2 lambda_3 lambda_4 
    54.0     28.0     11.5      6.5 

Interpretácia. Prvé vlastné číslo vysvetľuje 54 % celkovej variancie, prvé dve spolu 82 %. Všetky zložky prvého vlastného vektora majú rovnaké znamienko – ide o spoločný pohyb všetkých štyroch akcií. Najväčšie váhy v ňom majú TSLA a NVDA, teda akcie s najväčším rozptylom. V druhom vlastnom vektore majú NVDA a TSLA opačné znamienka: tento smer opisuje situácie, keď jedna z týchto akcií rastie a druhá klesá.

---

Malé cvičenie

Moje vypracovanie: Z denných close cien štyroch európskych akciových indexov vypočítajte denné výnosy a ich variančno-kovariančnú maticu. Nájdite jej vlastné čísla a vlastné vektory, overte ich vlastnosti a vysvetlite, čo hovoria o riziku (variancii) kombinácií týchto indexov.

1. Údaje: denné close ceny a výnosy

Použijeme dataset EuStockMarkets, ktorý je súčasťou základného R (netreba nič sťahovať ani inštalovať). Obsahuje denné close ceny indexov DAX (Nemecko), SMI (Švajčiarsko), CAC (Francúzsko) a FTSE (Veľká Británia) z rokov 1991 – 1998.

Variančno-kovariančnú maticu nepočítame z cien, ale z výnosov. Ceny majú dlhodobý trend a ich kovariancie by zachytávali hlavne spoločný rast v čase, nie to, ako sa indexy pohybujú spolu zo dňa na deň. Denný logaritmický výnos je

\[ r_t = \ln P_t - \ln P_{t-1}. \]

ceny <- EuStockMarkets            # denne close ceny (casovy rad so 4 stlpcami)
dim(ceny)                         # pocet dni x pocet indexov

[1] 1860    4

head(ceny, 3)

         DAX    SMI    CAC   FTSE
[1,] 1628.75 1678.1 1772.8 2443.6
[2,] 1613.63 1688.5 1750.5 2460.2
[3,] 1606.51 1678.6 1718.0 2448.2

vynosy <- as.matrix(diff(log(ceny)))   # denne logaritmicke vynosy ako obycajna matica
n <- nrow(vynosy)                      # pocet pozorovani (dni)
p <- ncol(vynosy)                      # pocet premennych (indexov)
c(n = n, p = p)

   n    p 
1859    4 

round(head(vynosy, 3), 5)

          DAX      SMI      CAC     FTSE
[1,] -0.00933  0.00618 -0.01266  0.00677
[2,] -0.00442 -0.00588 -0.01874 -0.00489
[3,]  0.00900  0.00327 -0.00578  0.00903

2. Variančno-kovariančná matica

Maticu vypočítame najprv „ručne“ pomocou maticových operácií a potom funkciou cov(). Ak od každého stĺpca matice výnosov odpočítame jeho priemer, dostaneme centrovanú maticu \(\mathbf X_c\) a platí

\[ \mathbf S = \frac{1}{n-1}\,\mathbf X_c^T\mathbf X_c . \]

priemery <- colMeans(vynosy)               # priemerny denny vynos kazdeho indexu
Xc <- sweep(vynosy, 2, priemery)           # od kazdeho stlpca odpocitame jeho priemer
S_rucne <- t(Xc) %*% Xc / (n - 1)          # kovariancna matica pomocou matic

S <- cov(vynosy)                           # to iste pomocou funkcie cov()
S

              DAX          SMI          CAC         FTSE
DAX  1.061072e-04 6.699564e-05 8.345130e-05 5.241794e-05
SMI  6.699564e-05 8.556317e-05 6.285881e-05 4.304517e-05
CAC  8.345130e-05 6.285881e-05 1.216802e-04 5.693174e-05
FTSE 5.241794e-05 4.304517e-05 5.693174e-05 6.332543e-05

max(abs(S - S_rucne))                      # rozdiel je prakticky nulovy

[1] 3.388132e-19

isSymmetric(S)                             # kovariancna matica je symetricka

[1] TRUE

Na diagonále sú rozptyly (variancie) denných výnosov jednotlivých indexov, mimo diagonály kovariancie dvojíc indexov. Všetky kovariancie sú kladné, indexy sa teda v priemere pohybujú rovnakým smerom.

3. Vlastné čísla a vlastné vektory

Vlastný vektor je určený až na znamienko (\(\mathbf v\) aj \(-\mathbf v\) sú správne). Pre ľahšie čítanie preto prvý vektor otočíme tak, aby súčet jeho zložiek bol kladný.

eig_S <- eigen(S, symmetric = TRUE)

lambda_S <- eig_S$values
names(lambda_S) <- paste0("lambda_", 1:p)
lambda_S

    lambda_1     lambda_2     lambda_3     lambda_4 
2.845255e-04 3.881170e-05 2.796618e-05 2.537260e-05 

V_S <- eig_S$vectors
if (sum(V_S[, 1]) < 0) V_S[, 1] <- -V_S[, 1]   # otocenie znamienka prveho vektora
dimnames(V_S) <- list(colnames(vynosy), paste0("v", 1:p))
round(V_S, 4)

         v1      v2      v3      v4
DAX  0.5553  0.1770  0.3893  0.7132
SMI  0.4537  0.7139  0.0053 -0.5333
CAC  0.5896 -0.6742  0.1965 -0.3990
FTSE 0.3716 -0.0664 -0.8999  0.2183

4. Overenie vlastností

Pre variančno-kovariančnú maticu \(\mathbf S\) s maticou vlastných vektorov \(\mathbf V\) a diagonálnou maticou vlastných čísel \(\boldsymbol\Lambda\) overíme:

1. \(\mathbf S\mathbf v_i = \lambda_i\mathbf v_i\) pre všetky páry naraz, t. j. \(\mathbf S\mathbf V = \mathbf V\boldsymbol\Lambda\),
2. vlastné vektory sú ortonormálne: \(\mathbf V^T\mathbf V = \mathbf I\),
3. spektrálny rozklad: \(\mathbf S = \mathbf V\boldsymbol\Lambda\mathbf V^T\),
4. súčet vlastných čísel = stopa matice, súčin vlastných čísel = determinant,
5. všetky vlastné čísla sú kladné, teda matica je kladne definitná.

Lambda_S <- diag(lambda_S)                       # diagonalna matica vlastnych cisel

# 1. S V = V Lambda
max(abs(S %*% V_S - V_S %*% Lambda_S))

[1] 1.084202e-19

# 2. V^T V = I
round(t(V_S) %*% V_S, 10)

   v1 v2 v3 v4
v1  1  0  0  0
v2  0  1  0  0
v3  0  0  1  0
v4  0  0  0  1

# 3. S = V Lambda V^T
max(abs(S - V_S %*% Lambda_S %*% t(V_S)))

[1] 9.486769e-20

# 4. stopa a determinant
c(sucet_vlastnych_cisel = sum(lambda_S), stopa_matice = sum(diag(S)))

sucet_vlastnych_cisel          stopa_matice 
          0.000376676           0.000376676 

c(sucin_vlastnych_cisel = prod(lambda_S), determinant = det(S))

sucin_vlastnych_cisel           determinant 
         7.835779e-18          7.835779e-18 

# 5. kladna definitnost
all(lambda_S > 0)

[1] TRUE

Stopa kovariančnej matice je súčet rozptylov všetkých indexov, teda celková variancia. Vlastné čísla túto celkovú varianciu rozdeľujú na časti, z ktorých každá patrí jednému smeru (vlastnému vektoru). Variančno-kovariančná matica nemôže mať záporné vlastné čísla, pretože rozptyl nemôže byť záporný. Nulové vlastné číslo by znamenalo, že niektorá kombinácia indexov má nulový rozptyl (jeden index by sa dal presne vyjadriť pomocou ostatných). Tu sú všetky vlastné čísla kladné.

5. Podiel vysvetlenej variancie

podiel_S <- lambda_S / sum(lambda_S)
round(100 * rbind(podiel = podiel_S, kumulativne = cumsum(podiel_S)), 1)

            lambda_1 lambda_2 lambda_3 lambda_4
podiel          75.5     10.3      7.4      6.7
kumulativne     75.5     85.8     93.3    100.0

stlpce <- barplot(100 * podiel_S,
                  names.arg = parse(text = paste0("lambda[", 1:p, "]")),
                  col = "steelblue4", border = NA, ylim = c(0, 100),
                  ylab = "Podiel celkovej variancie (%)",
                  main = "Vlastné čísla kovariančnej matice výnosov")
text(stlpce, 100 * podiel_S, paste0(round(100 * podiel_S, 1), " %"), pos = 3)

Podiel celkovej variancie pripadajúci na jednotlivé vlastné čísla.

Prvé vlastné číslo vysvetľuje až 75.5 % celkovej variancie denných výnosov. Zvyšné tri vlastné čísla sú výrazne menšie a spolu vysvetľujú len 24.5 %.

6. Interpretácia vlastných vektorov

barplot(t(V_S[, 1:2]), beside = TRUE,
        col = c("steelblue4", "darkorange2"), border = NA,
        ylim = c(-1, 1), ylab = "Hodnota zložky vlastného vektora",
        main = "Prvý a druhý vlastný vektor")
abline(h = 0, col = "grey50")
legend("bottomleft", legend = c("v1", "v2"),
       fill = c("steelblue4", "darkorange2"), border = NA, bty = "n")

Zložky prvých dvoch vlastných vektorov kovariančnej matice.

• Prvý vlastný vektor má všetky zložky s rovnakým znamienkom a podobnej veľkosti. Opisuje spoločný pohyb celého európskeho trhu – keď rastie jeden index, zvyčajne rastú aj ostatné. Tento „trhový“ smer nesie najväčšiu časť rizika.
• Druhý vlastný vektor má veľké zložky pri SMI a CAC s opačnými znamienkami. Opisuje situácie, keď sa švajčiarsky a francúzsky trh pohybujú odlišne.
• Ďalšie vlastné vektory opisujú menšie, „špecifické“ pohyby (napr. tretí vektor je tvorený hlavne indexom FTSE).

7. Vlastné čísla ako rozptyl kombinácie indexov

Ak vytvoríme kombináciu výnosov \(\mathbf w^T\mathbf r\), kde \(\mathbf w\) je vektor s dĺžkou 1, jej rozptyl je

\[ \operatorname{var}(\mathbf w^T\mathbf r) = \mathbf w^T\mathbf S\,\mathbf w . \]

Pre \(\mathbf w = \mathbf v_i\) dostaneme \(\mathbf v_i^T\mathbf S\mathbf v_i = \lambda_i\mathbf v_i^T\mathbf v_i = \lambda_i\). Navyše žiadny vektor s dĺžkou 1 nemôže mať rozptyl väčší ako \(\lambda_1\) ani menší ako \(\lambda_p\). Overíme to na 10 000 náhodných smeroch.

rozptyl_smeru <- function(w) as.numeric(t(w) %*% S %*% w)   # w^T S w

# rozptyl pozdlz vlastnych vektorov sa rovna vlastnym cislam
sapply(1:p, function(i) rozptyl_smeru(V_S[, i]))

[1] 2.845255e-04 3.881170e-05 2.796618e-05 2.537260e-05

lambda_S

    lambda_1     lambda_2     lambda_3     lambda_4 
2.845255e-04 3.881170e-05 2.796618e-05 2.537260e-05 

# kontrola priamo z dat: rozptyl casoveho radu  r_t^T v_1
var(as.vector(vynosy %*% V_S[, 1]))

[1] 0.0002845255

# nahodne smery s dlzkou 1
set.seed(2026)
W <- matrix(rnorm(10000 * p), ncol = p)
W <- W / sqrt(rowSums(W^2))                   # kazdy riadok normalizujeme na dlzku 1
rozptyly_W <- apply(W, 1, rozptyl_smeru)
c(min_nahodne = min(rozptyly_W), lambda_min = min(lambda_S))

 min_nahodne   lambda_min 
2.543713e-05 2.537260e-05 

c(max_nahodne = max(rozptyly_W), lambda_max = max(lambda_S))

 max_nahodne   lambda_max 
0.0002841523 0.0002845255 

hist(rozptyly_W, breaks = 50, col = "grey80", border = "white",
     xlim = range(c(0, lambda_S)),
     xlab = "Rozptyl kombinácie výnosov", ylab = "Počet smerov",
     main = "Rozptyl náhodných smerov")
abline(v = lambda_S[c(1, p)], col = c("darkorange2", "steelblue4"), lwd = 3)
legend("top", legend = c("lambda_1 (maximum)", "lambda_4 (minimum)"),
       col = c("darkorange2", "steelblue4"), lwd = 3, bty = "n")

Rozptyl kombinácií výnosov pre 10 000 náhodných smerov s dĺžkou 1.

Prvý vlastný vektor je teda smer s najväčším rizikom (rozptylom), posledný vlastný vektor smer s najmenším rizikom. Posledný vektor kombinuje indexy s kladnými aj zápornými znamienkami – riziko sa znižuje tým, že sa pohyby jednotlivých indexov navzájom vyrušia (diverzifikácia).

8. Nové nekorelované premenné

Ak výnosy premietneme do smerov vlastných vektorov, \(\mathbf Z = \mathbf X_c\mathbf V\), dostaneme nové premenné, ktoré sú navzájom nekorelované a ich rozptyly sú práve vlastné čísla: \(\operatorname{cov}(\mathbf Z) = \boldsymbol\Lambda\). Na tomto princípe stojí analýza hlavných komponentov (PCA), ktorú v R počíta funkcia prcomp().

Z <- Xc %*% V_S                          # vynosy v suradniciach vlastnych vektorov

# kovariancna matica novych premennych v jednotkach 10^-6:
# mimo diagonaly su nuly, na diagonale vlastne cisla
round(cov(Z) * 1e6, 4)

         v1      v2      v3      v4
v1 284.5255  0.0000  0.0000  0.0000
v2   0.0000 38.8117  0.0000  0.0000
v3   0.0000  0.0000 27.9662  0.0000
v4   0.0000  0.0000  0.0000 25.3726

round(lambda_S * 1e6, 4)

lambda_1 lambda_2 lambda_3 lambda_4 
284.5255  38.8117  27.9662  25.3726 

# porovnanie s hotovou funkciou prcomp()
pca <- prcomp(vynosy)
max(abs(pca$sdev^2 - lambda_S))                    # rovnake vlastne cisla

[1] 7.589415e-19

max(abs(abs(pca$rotation) - abs(V_S)))             # rovnake vektory (az na znamienko)

[1] 5.218048e-15

9. Geometrické znázornenie pre dva indexy

Podobne ako pri matici \(\mathbf A\) na začiatku sa pozrieme na vlastné vektory graficky. Zoberieme len dva indexy, DAX a CAC, a ich kovariančnú maticu \(2 \times 2\). Šípky ukazujú smery vlastných vektorov; ich dĺžka je \(2\sqrt{\lambda_i}\), teda dve smerodajné odchýlky v danom smere.

S2 <- S[c("DAX", "CAC"), c("DAX", "CAC")]       # podmatica 2x2
eig_S2 <- eigen(S2, symmetric = TRUE)
eig_S2$values

[1] 1.977075e-04 3.007995e-05

eig_S2$vectors

          [,1]       [,2]
[1,] 0.6734604 -0.7392233
[2,] 0.7392233  0.6734604

plot(vynosy[, "DAX"], vynosy[, "CAC"], asp = 1, pch = 16, cex = 0.6,
     col = adjustcolor("grey40", alpha.f = 0.35),
     xlab = "Denný výnos DAX", ylab = "Denný výnos CAC",
     main = "Vlastné vektory určujú hlavné smery rozptylu")
abline(h = 0, v = 0, col = "grey80")
farby <- c("darkorange2", "steelblue4")
for (i in 1:2) {
  dlzka <- 2 * sqrt(eig_S2$values[i])
  arrows(priemery["DAX"], priemery["CAC"],
         priemery["DAX"] + dlzka * eig_S2$vectors[1, i],
         priemery["CAC"] + dlzka * eig_S2$vectors[2, i],
         col = farby[i], lwd = 4, length = 0.12)
}
legend("topleft", legend = c("1. vlastný vektor", "2. vlastný vektor"),
       col = farby, lwd = 4, bty = "n")

Denné výnosy DAX a CAC a vlastné vektory ich kovariančnej matice.

Mrak bodov je pretiahnutý pozdĺž prvého vlastného vektora – v tomto smere majú výnosy najväčší rozptyl. Druhý vlastný vektor je naň kolmý a rozptyl v jeho smere je oveľa menší (pomer vlastných čísel je 6.6 : 1).

---