Diskriminačné metódy ako napr. logistická regresia modelujú podmienenú pravdepodobnosť Pr[Y=y|\vec X= \mathbf{x}] priamo.
Generatívne metódy idú na problém nepriamo:
Ak je rozdelenie \vec X normálne, oba prístupy dávajú veľmi podobný výsledok, generatívny je však
Generatívny prístup je jednoducho a prirodzene rozšíriteľný na viac ako 2 triedy.
Nech
Potom podľa Bayesovej vety platí vzťah pre pravdepodobnosť, že pozorovanie \vec X=\mathbf{x} patrí do k-tej triedy, Pr[Y=c_k|\vec X=\mathbf{x}] = \frac{\pi_k f_k(\mathbf{x})}{\sum_{\ell=1}^K \pi_\ell f_\ell(\mathbf{x})}. Nazýva sa aj aposteriórnou pravdepodobnosťou.
To znamená, že namiesto priameho odhadu Pr[Y=c_k|\vec X=\mathbf{x}] (ako pri diskriminatívnom prístupe) nám stačí odhadnúť \pi_k, f_k a dosadiť do vzťahu.
V nasledujúcom si preberieme tri (generatívne) klasifikátory, ktoré používajú rôzne odhady hustoty f_k.
V lineárnej diskriminačnej analýze (LDA) predpokladáme, že f_k je hustota normálneho rozdelenia, f_k(x) = \frac{1}{\sqrt{2\pi}\sigma_k}\exp\left(-\frac{(x-\mu_k)^2}{2\sigma_k^2}\right), kde \mu_k a \sigma_k^2 sú stredná hodnota a rozptyl prediktoru X v k-tej triede.
Ďalej predpokladáme, že \sigma_1=\ldots=\sigma_K=\sigma.
Po dosadení takto špecifikovanej hustoty do Bayesovej vety dostávame Pr[Y=c_k|X=x] = \frac{\pi_k \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu_k)^2}{2\sigma^2}\right)}{\sum_{\ell=1}^K \pi_\ell \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu_\ell)^2}{2\sigma^2}\right)}.
Keďže (bayesovský) klasifikátor zaradí pozorovanie X=x do tej triedy, ktorej zodpovedá najvyššia aposteriórna pravdepodobnosť, stačí pravdepodobnosť nahradiť jednoduchšou funkciou. Tou je tzv. diskriminačná funkcia \delta_k. Vznikne logaritmovaním vzťahu pre Pr[Y=c_k|X=x] a vynechaním členov spoločných pre všetky triedy, takže proporcionalita medzi skupinami zostáva zachovaná, \delta_k(x) \propto Pr[Y=c_k|X=x]. Po úpravách dostávame \delta_k(x) = x\frac{\mu_k}{\sigma^2} - \frac{\mu_k^2}{2\sigma^2} + \ln\pi_k. Diskriminačná funkcia je lineárna v premennej x, čo vysvetľuje názov metódy.
Parametre skutočného rozdelenia sotvakedy poznáme. Najlepšie odhady, čo máme, sú výberový priemer a (vážený) rozptyl \begin{split} \hat\mu_k &= \frac1{n_k}\sum_{i:\ y_i=c_k}x_i\ , \\ \hat\sigma^2 &= \frac1{n-K}\sum_{k=1}^K\ \sum_{i:\ y_i=c_k}(x_i-\hat\mu_k)^2=\frac1{n-K}\sum_{k=1}^K(n_k-1)\sigma_k^2\ . \end{split}
pars <- list(list(mean=-1.25, sd=1),
list(mean=1.25, sd=1))
ggplot() +
geom_function(fun = dnorm, args = pars[[1]], color = 2, linewidth = 1.5) +
geom_function(fun = dnorm, args = pars[[2]], color = 4, linewidth = 1.5) +
geom_vline(xintercept = 0, linetype=2) +
xlab("x") + ylab("f(x)") +
xlim(-5, 5)
# data
set.seed(4)
dat <- list("k=1" = data.frame(x = rnorm(20, mean = pars[[1]]$mean, sd = pars[[1]]$sd)),
"k=2" = data.frame(x = rnorm(20, mean = pars[[2]]$mean, sd = pars[[2]]$sd))
) |>
dplyr::bind_rows(.id = "class")
# súhrny
tmp1 <- dat |>
dplyr::group_by(class) |>
dplyr::summarise(mu = mean(x),
s2 = var(x),
n = dplyr::n()
) |>
dplyr::ungroup()
tmp2 <- tmp1 |>
dplyr::summarise(mu = weighted.mean(mu, w = n),
s2 = weighted.mean(s2, w = n-1), # pooled variance
n = sum(n)
)
# priesečník
tmp2$boundary <- tmp2$mu - diff(log(tmp1$n/tmp2$n)) * tmp2$s2 / diff(tmp1$mu)
# report
dplyr::full_join(tmp1, rbind(tmp2), dplyr::join_by(mu, s2, n)) |>
pander::pander()| class | mu | s2 | n | boundary |
|---|---|---|---|---|
| k=1 | -0.8737 | 0.6421 | 20 | NA |
| k=2 | 1.556 | 0.8149 | 20 | NA |
| NA | 0.3411 | 0.7285 | 40 | 0.3411 |
# graf vyhladenej empirickej hustoty
dat |>
ggplot() + aes(x = x, fill = class) +
geom_density(alpha = 0.5, adjust = 1.5) +
geom_vline(xintercept = tmp2$boundary, linetype = 3)
# graf diskriminačných funkcií
LDfun <- Vectorize(
function(x, mu, s2, pi) x * mu/s2 - mu^2/(2*s2) + log(pi),
vectorize.args = "x")
xgrid <- c(-3,3)
LDfun(xgrid, mu = tmp1$mu, s2 = tmp2$s2, pi = tmp1$n/tmp2$n) |>
t() |> as.data.frame() |> setNames(tmp1$class) |>
cbind(x = xgrid) |>
tidyr::pivot_longer(-x, names_to = "class", values_to = "delta") |>
ggplot() + aes(x = x, y = delta, group = class, color = class) + geom_line() +
geom_vline(xintercept = tmp2$boundary, linetype = 3)
\begin{split}
L(\mathbf{u}) = \mathbf{u}^T\mathbf{B}\mathbf{u} - \lambda(\mathbf{u}^T\mathbf{V}\mathbf{u}-1) & \\
\frac{\partial L(\mathbf{u})}{\partial \mathbf{u}} = 2\mathbf{B}\mathbf{u} - 2\lambda\mathbf{V}\mathbf{u} \qquad\qquad & \\
\qquad\qquad 2\mathbf{B}\hat{\mathbf{u}} - 2\lambda\mathbf{V}\hat{\mathbf{u}} & = 0 \\
\mathbf{B}\mathbf{u} & = \lambda\mathbf{V}\mathbf{u} \\
\mathbf{V}^{-1}\mathbf{B}\mathbf{u} & = \lambda\mathbf{u}
\end{split}
takže náš hľadaný optimálny vektor \mathbf{u} je vlastný vektor matice \mathbf{V}^{-1}\mathbf{B}.
(Pre detaily riešenia problému zovšeobecnených vlastných čísel pozri Ghojogh et al., 2022.)
# odhad modelu
fit <- MASS::lda(Species ~ Petal.Width + Petal.Length, data = iris)
fitCall:
lda(Species ~ Petal.Width + Petal.Length, data = iris)
Prior probabilities of groups:
setosa versicolor virginica
0.3333333 0.3333333 0.3333333
Group means:
Petal.Width Petal.Length
setosa 0.246 1.462
versicolor 1.326 4.260
virginica 2.026 5.552
Coefficients of linear discriminants:
LD1 LD2
Petal.Width 2.402394 5.042599
Petal.Length 1.544371 -2.161222
Proportion of trace:
LD1 LD2
0.9947 0.0053 # zobrazenie tried v priestore prediktorov aj s centroidmi
centr.df <- fit |>
getElement("means") |> # treba ich ešte
tibble::as_tibble(rownames = "Species")
iris |>
ggplot() +
aes(x = Petal.Width, y = Petal.Length, color = Species) +
geom_point() +
geom_point(data = centr.df, shape = 4, size = 3, show.legend = FALSE)
# prvky transformačných vektorov "u" pre zobrazenie v priestore prediktorov
trans.df <- fit |>
coef() |> # to isté ako getElement("scaling")
# apply(2, function(x) x/sqrt(sum(x^2))) |> # konverzia na jednotkové vektory (voliteľné)
t() |>
as.data.frame() |>
# setNames(fit$scaling |> rownames()) |> # lebo apply zruší názvy riadkov
cbind(variable = fit$scaling |> colnames())
# biplot
iris |>
purrr::map_if(is.numeric, scale) |> # centrovanie a škálovanie
as.data.frame() |>
ggplot() + aes(x = Petal.Width, y = Petal.Length) +
geom_point(aes(color = Species)) +
geom_segment(aes(x = 0, y = 0, xend = Petal.Width, yend = Petal.Length),
data = trans.df,
arrow = arrow(length = unit(0.5, "char"))
) +
geom_text(aes(x = Petal.Width, y = Petal.Length, label = variable),
data = trans.df,
nudge_x = 0.5)
# biplot pomocou externého balíku
# library(devtools)
# install_github('fawda123/ggord')
# ggord::ggord(fit, iris$Species, size = 2)
# prvky transformačných vektorov "u" pre zobrazenie v priestore diskriminantov
trans.df <- fit |>
getElement("scaling") |>
as.data.frame() |>
tibble::rownames_to_column(var = "variable")
# výpočet hodnôt diskriminantov (aj pomocou predict(fit)$x)
LD <- iris[,4:3] |>
as.matrix() |>
(`%*%`)(fit$scaling) |> # transformácia dát do priestoru diskriminantov
apply(2, scale, center = TRUE, scale = FALSE) |> # centrovanie diskriminantov
cbind.data.frame(iris[5])
# biplot v priestore diskriminantov
LD |>
ggplot() + aes(LD1, LD2) +
geom_point(aes(color = Species)) +
geom_segment(aes(x = 0, y = 0, xend = LD1, yend = LD2),
data = trans.df,
arrow = arrow(length = unit(0.5, "char"))
) +
geom_text(aes(x = LD1, y = LD2, label = variable),
data = trans.df,
nudge_y = -0.5) +
stat_ellipse(aes(fill = Species), geom = "polygon", alpha = 0.2) +
ylim(-3, 6)
# hustota pravdepodobnosti podľa tried v priestore diskriminantov
LD |>
tidyr::pivot_longer(-Species, names_to = "diskriminant", values_to = "value") |>
ggplot() + aes(x = value, fill = Species) + geom_density(alpha = 0.5, adjust = 2) + facet_wrap(~diskriminant)
# empirická hustota v triedach pre prvý diskriminant vstavanou funkciou
# MASS::ldahist(iris[3:4], iris$Species)# predikcia v trénovacej vzorke (obsahuje triedy, ich pravdepodobnosti aj hodnoty diskriminantov)
# pred <- predict(fit)
# (urpbiť: vykresliť pravdepodobnosti pre jednotlivé triedy a diskrimininanty ako alternatívu ku vyhladeným frekvenčným polynómom hore)
# nové dáta a predpovede
xgrid <- expand.grid(Petal.Length = seq(1,7,0.01), Petal.Width = seq(0,2.5,0.01))
pred <- cbind.data.frame(
xgrid,
class = predict(fit, newdata = xgrid) |> getElement("class")
)
# zobrazenie oblastí klasifikácie do jednotlivých tried
pred |>
ggplot() + aes(x = Petal.Width, y = Petal.Length) + geom_tile(aes(fill = class), alpha = 0.3) +
geom_point(aes(color = Species), data = iris) +
ggtitle("Hranice tried podľa LDA")
# externým balíkom
# klaR::partimat(Species ~ Petal.Width + Petal.Length, data = iris, method = "lda")# kvadratická diskriminačná analýza
fit <- MASS::qda(Species ~ Petal.Width + Petal.Length, data = iris)
fitCall:
qda(Species ~ Petal.Width + Petal.Length, data = iris)
Prior probabilities of groups:
setosa versicolor virginica
0.3333333 0.3333333 0.3333333
Group means:
Petal.Width Petal.Length
setosa 0.246 1.462
versicolor 1.326 4.260
virginica 2.026 5.552# nové dáta a predpovede
xgrid <- expand.grid(Petal.Length = seq(1,7,0.01),
Petal.Width = seq(0,2.5,0.01)
)
pred <- cbind.data.frame(
xgrid,
class = predict(fit, newdata = xgrid) |> getElement("class")
)
# zobrazenie oblastí klasifikácie do jednotlivých tried
pred |>
ggplot() + aes(x = Petal.Width, y = Petal.Length) + geom_tile(aes(fill = class), alpha = 0.3) +
geom_point(aes(color = Species), data = iris) +
ggtitle("Hranice tried podľa QDA")
# --- balíkom klaR
# fit <- klaR::NaiveBayes(Species ~ Petal.Width + Petal.Length, data = iris, usekernel = T, adjust = 2.) # výstup je neformátovaný, list
# plot(fit)
# --- balíkom e1071
# fit <- e1071::naiveBayes(Species ~ Petal.Width + Petal.Length, data = iris)
# výstup tvorí vektor apriórnych pravdepodobností a pod názvom "Conditional probabilities" stredná hodnota a štandardná odchýlka prediktorov v jednotlivých triedach (žiaľ neoznačné).
# --- balíkom naivebayes
# odhad modelu s predpokladom normálneho (parametrického) a vyhladeného (neparametrického) rozdelenia. Argument adjust násobí prednastavenú šírku pásma jadrového vyhladenia.
fit <- list(
gauss = naivebayes::naive_bayes(Species ~ Petal.Width + Petal.Length, data = iris),
kernel = naivebayes::naive_bayes(Species ~ Petal.Width + Petal.Length, data = iris,
usekernel = TRUE, adjust = 2)
)# zobrazenie vyhladenej hustoty v jednotlivých triedach
fit$kernel |> plot(arg.num = list(lwd = 2))
xgrid <- expand.grid(Petal.Length = seq(1,7,0.02),
Petal.Width = seq(0,2.5,0.02)
)
fit |>
purrr::map_dfc(predict, newdata = xgrid, type = "class") |>
cbind(xgrid) |>
tidyr::pivot_longer(c(gauss, kernel), names_to = "distrib", values_to = "class") |>
ggplot() + aes(x = Petal.Width, y = Petal.Length) + geom_tile(aes(fill = class), alpha = 0.3) +
geom_point(aes(color = Species), data = iris) +
facet_wrap(~distrib) +
ggtitle("Hranice tried NB podľa predpokladu rozdelenia prediktorov")