Pokročilé metódy štatistického modelovania

Author

Tomáš Bacigál

Published

December 4, 2023

Výstavba modelov združeného rozdelenia pravdepodobnosti a ich aplikácia

Dokument obsahuje postup konštrukcie združeného rozdelenia náhodného vektora pomocou súboru pozorovaní,

od zoznámenia sa s jednotlivými náhodnými premennými a odhadom ich rozdelenia pravdepodobnosti,
cez transformáciu do jednotkového intervalu pomocou odhadnutých jednorozmerných CDF,
odhad modelov kopuly,
spojenie kopuly a marginálnych rozdelení,
výpočet (nepodmienených i podmienených) pravdepodobností, kvantilov a strednej hodnoty,
zobrazenie rozdelenia premennej v role odozvy prostredníctvom podmienenej funkcie hustoty.

Prieskumná analýza

Technický úvod

Knižnice a funkcie

Code

library(ggplot2)

# make regular sequence over the range of vector x  (similar to modelr::seq_range)
seq_range <- function(x, length = length(x)) {
  ran <- range(x)
  seq.int(ran[1], ran[2], length.out = length)
}

Načítanie údajov

Code

dat <- mtcars |>
  subset(select = c(disp, mpg))

Jednorozmerná

empirická hustota s jadrovým vyhladením

Code

dat |>
  tidyr::pivot_longer(cols = everything(), 
                      names_to = "variable", values_to = "value") |>
  ggplot() + aes(x = value) + 
  geom_histogram(aes(y = after_stat(density)), bins = 5, color="white") +
  stat_density(geom="line", adjust = 1.) +
  facet_wrap(~ variable, scales = "free")

číselné charakteristiky

Code

dat |>
  dplyr::summarise(across(.cols = c(mpg, disp), 
                          .fns = list(mean = mean, sdev = sd)
                          )
  )

  mpg_mean mpg_sdev disp_mean disp_sdev
1 20.09062 6.026948  230.7219  123.9387

Code

# alternatívne, pri väčšom počte premenných alebo charakteristík:
#dat |>
#  tidyr::pivot_longer(cols = everything(), 
#                      names_to = "variable", values_to = "value") |>
#  dplyr::group_by(variable) |>
#  dplyr::summarise(across(.cols = value, 
#                          .fns = list(mean = mean, sdev = sd),
#                          .names = "{.fn}"
#  ))

test normality (p-hodnota)

Code

dat |> sapply(function(x) shapiro.test(x)$p.value) |> round(3)

 disp   mpg 
0.021 0.123

Závislosť

bodový graf

Code

ggplot(dat) + aes(x = disp, y = mpg) + geom_point()

korelačný koeficient

Code

sapply(c("pearson", "kendall", "spearman"), function(x) cor(dat, method = x)[1,2])

   pearson    kendall   spearman 
-0.8475514 -0.7681311 -0.9088824

test nulovosti korelačného koeficientu

Code

cor.test(~ disp + mpg, dat)     # cor.test(~ ., dat)$p.value |> round(3)


    Pearson's product-moment correlation

data:  disp and mpg
t = -8.7472, df = 30, p-value = 9.38e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9233594 -0.7081376
sample estimates:
       cor 
-0.8475514

Spolu

Code

ggpubr::ggscatterhist(dat, x = "disp", y = "mpg", 
                      margin.plot = "histogram", bins = 5,
                      margin.params = list(fill = "lightgray"),
)
# Pridaním ďalšej premennej s farebným odlíšením by sa dal pozorovať jej vplyv zvlášť na marginálne rozdelenia a na vzťah medzi disp a mpg.

Odhad modelu

Code

model <- list()  # kontajner

Empirické rozdelenie

Code

model$emp <- list(
  cdf = function(x) {
    dat |>  # nxp
      apply(MARGIN = 1, '<=', x) |>  # returns pxn matrix
      apply(2, all) |>  # that's why MARGIN = 2
      mean()
  }
)
# model$emp$cdf(c(300, 20))  # alebo copula::F.n(c(300, 20), as.matrix(dat))

Združené normálne rozdelenie

odhad parametrov metódou momentov

Code

mu <- sapply(dat, mean)
Sigma <- cov(dat)
cbind(mean = mu, Sigma) |> round(1)  # spojený výpis strednej hodnoty a kovariančnej matice

      mean    disp    mpg
disp 230.7 15360.8 -633.1
mpg   20.1  -633.1   36.3

poskladanie PDF a CDF

Code

model$norm <- list(
  pdf = function(x) mvtnorm::dmvnorm(x, mean = mu, sigma = Sigma),
  cdf = function(x) apply(rbind(x), 1, function(y)
    mvtnorm::pmvnorm(upper = y, mean = mu, sigma = Sigma)) |> as.numeric()
)
# model$norm$pdf(c(300, 20))

Jadrovo vyhladené empirické rozdelenie

Code

# Balík probhat (spolu so závislosťami barsurf, intoo a kubik) treba stiahnuť z archívu CRAN a nainštalovať zo súboru. 
model$smoo <- list(
  pdf = probhat::pdfmv.cks(as.matrix(dat), a = c(0,0)),
  cdf = probhat::cdfmv.cks(as.matrix(dat), a = c(0,0))
)
# model$smoo$cdf(c(300,20))

Rozkladom na marginálie a kopulu

Marginálne rozdelenia

Neparametrické

Neparametrickým odhadom je empirické rozdelenie. Reprezentované je tabuľkou relatívnych početností alebo distribučnou funkciou (so schodkovitým grafom)

Code

model$marg$mpg$emp$cdf <- ecdf(dat$mpg)
model$marg$disp$emp$cdf <- ecdf(dat$disp)
#plot(model$marg$mpg$emp$cdf, from=0, to=30)

Parametrické

Najprv iba ilustrujeme rôzne metódy odhadu parametrov rozdelenia.

Pri premennej mpg predpokladajme, že pochádza z normálneho rozdelenia.

Code

fit_norm <-  with(dat, {
  start <- c(mean = mean(mpg), sd = sd(mpg))
  lmoments <- lmomco::lmoms(mpg)  # alebo lmom::samlmu(mpg)
  rbind(
    mom = start,
    mle = MASS::fitdistr(mpg, densfun = dnorm, 
                         start = as.list(start))$estimate,
    lse = optim(par = start, 
                fn = function(pars) sum((model$marg$mpg$emp$cdf(mpg) - pnorm(mpg, mean=pars[1], sd=pars[2]))^2))$par,
    lmo = lmomco::lmoms(mpg) |> lmomco::parnor() |> getElement("para") 
    # lmo = lmom::samlmu(dat$mpg) |> lmom::pelnor()  # balík je od autora metódy
    )
})
fit_norm

        mean       sd
mom 20.09062 6.026948
mle 20.09062 5.932030
lse 19.02397 5.378611
lmo 20.09062 6.022948

Odhady sú viac-menej podobné, najviac sa líši odhad minimalizáciou štvorca vzdialenosti ku empirickej CDF (pre ktorú možno existuje aj optimalizovaná verzia). Vezmeme najjednoduchší odhad, teda metódou momentov.

Code

model$marg$mpg$nor <- fit_norm["mom",] |> 
  (\(par) list(
    pdf = function(x) dnorm(x, mean=par[1], sd=par[2]),
    cdf = function(x) pnorm(x, mean=par[1], sd=par[2]),
    qua = function(x) qnorm(x, mean=par[1], sd=par[2])
))()
# plot(model$marg$mpg$nor$cdf, from=0, to=40, col=4) 
# lines(model$marg$mpg$emp$cdf, from=0, to=40)

Pre premennú disp, pri ktorej bola zamietnutá hypotéza o normalite, si zvolíme (minimálne) dve konkurenčné parametrické triedy rozdelenia, medzi ktorými sa potom rozhodneme, napr. log-normálne a exponenciálne. (Pre kontrolu metódy rozkladu pomocou kopúl odhadneme aj normálne rozdelenie, aby sa takto zložené rozdelenie dalo porovnať s normálnym).

Code

model$marg$disp$lnor <- dat$disp |>  
  (\(x) c(meanlog = mean(log(x)), sdlog = sd(log(x))))() |> # odhad metódou momentov
  (\(par) list(
    pdf = function(x) dlnorm(x, meanlog=par[1], sdlog=par[2]),
    cdf = function(x) plnorm(x, meanlog=par[1], sdlog=par[2]),
    qua = function(x) qlnorm(x, meanlog=par[1], sdlog=par[2])
))()

model$marg$disp$exp <- dat$disp |>  
  (\(x) c(lambda = 1 / mean(x)))() |>    # odhad metódou momentov
  (\(par) list(
      pdf = function(x) dexp(x, rate = par),
      cdf = function(x) pexp(x, rate = par),
      qua = function(x) qexp(x, rate = par)
))()

model$marg$disp$nor <- dat$disp |>  
  (\(x) c(mean = mean(x), sd = sd(x)))() |>  
  (\(par) list(
    pdf = function(x) dnorm(x, mean=par[1], sd=par[2]),
    cdf = function(x) pnorm(x, mean=par[1], sd=par[2]),
    qua = function(x) qnorm(x, mean=par[1], sd=par[2])
))()

Zobrazenie

Code

plot(model$marg$disp$lnor$cdf, from=0, to=500, col=4, ylab="CDF") 
plot(model$marg$disp$exp$cdf, from=0, to=500, col=2, add=T) 
plot(model$marg$disp$nor$cdf, from=0, to=500, col=3, add=T)
plot(model$marg$disp$emp$cdf, add=T)
legend("bottomright", legend = c("empirical", "log-normal", "exponential", "normal"), title = "distribution", col = c(1,4,2,3), lty=1, pch=c(19, NA, NA, NA))

Test dobrej zhody

Code

#ks.test(x = dat$disp, y = model$marg$disp$lnor$cdf)
#goftest::cvm.test(x = dat$disp, null = model$marg$disp$lnor$cdf)
sapply(c("lnor", "exp", "nor"), function(x) 
  goftest::cvm.test(dat$disp, 
                    null = model$marg$disp[[x]]$cdf
                    )[c("statistic","p.value")] |> unlist()
  ) |> round(3)

                  lnor   exp   nor
statistic.omega2 0.118 0.541 0.149
p.value          0.507 0.031 0.394

Exponenciálne rozdelenie nie je vhodné na popis rozdelenia disp. Keby testom prešli obe, rozhodli by sme sa aj podľa výšky testovacej štatistiky, pretože predstavuje mieru vzdialenosti od empirickej CDF (pre KS test je to maximálny rozdiel, zatiaľčo pre CvM súčet štvorcov rozdielov).

Kopula

Code

dat01 <- list()

Neparametrická

Najprv vypočítame pseudo-pozorovania U=F(X) pomocou empirických marginálnych CDF.

Code

dat01$emp <- dat |>
  transform(disp = model$marg$disp$emp$cdf(disp),
            mpg = model$marg$mpg$emp$cdf(mpg)
            )
# Kopula CDF
model$cop$emp$cdf <- function(x) copula::C.n(x, X = dat01$emp)
# model$cop$emp$cdf(c(0.5,1))

Zobrazenie pseudo-pozorovaní, čiže (prakticky) hustoty kopuly

Code

plot(dat01$emp, xlab="F(disp)", ylab="F(mpg)", main="Empirická hustota kopule")

Parametrická

Pseudopozorovania pomocou parametrických marginálnych CDF. Na odhad kopuly môžeme použiť aj dáta z neparametrických marginálnych CDF. Parametrické však budeme potrebovať pri spätnej transformácii na výpočet kvantilov a strednej hodnoty.

Code

dat01$par <- dat |>
  transform(disp = model$marg$disp$lnor$cdf(disp),
            mpg = model$marg$mpg$nor$cdf(mpg)
  )

Kopula normálneho rozdelenia (Gaussovská)

Code

model$cop$nor <- copula::normalCopula(param = -0.5) |>  # štartovacia kopula
  copula::fitCopula(data = as.matrix(dat01$par)) |>  # odhad
  slot(name="copula") |>   # extrakcia objektu 
  (\(obj) list(
    obj = obj,  # odhadnutý objekt
    pdf = function(x) copula::dCopula(x, copula = obj),
    cdf = function(x) copula::pCopula(x, copula = obj)
  ))()
# model$cop$nor$obj@parameters   # odhadnutý parameter kopule
# model$cop$nor$pdf(c(0.5,0.5))  # P(X<medianX, Y<medianY)

Preklopená (rotovaná) Gumbelova kopula

Code

model$cop$rgum <- copula::gumbelCopula(param = 2) |>  # štartovacia kopula
  copula::rotCopula(flip = c(TRUE, FALSE)) |>  # preklopenie cez u1=0.5
  copula::fitCopula(data = as.matrix(dat01$par)) |>  # odhad
  slot(name="copula") |>   # extrakcia kopule 
  (\(obj) list(
    obj = obj,  # fitted object
    pdf = function(x) copula::dCopula(x, copula = obj),
    cdf = function(x) copula::pCopula(x, copula = obj)
  ))()

Preklopená (rotovaná) Claytonova kopula (MLE zlyhala, preto inverziou tau)

Code

model$cop$rcla <- copula::claytonCopula(param = 1) |>  # štartovacia kopula
  copula::rotCopula(flip = c(FALSE, TRUE)) |>  # preklopenie cez u2=0.5
  copula::fitCopula(data = as.matrix(dat01$par), method="itau") |>  # odhad
  slot(name="copula") |>   # extrakcia kopule 
  (\(obj) list(
    obj = obj,  # fitted object
    pdf = function(x) copula::dCopula(x, copula = obj),
    cdf = function(x) copula::pCopula(x, copula = obj)
  ))()

simulácie z modelu

Code

model$cop$nor$obj |>
  copula::rCopula(n=100) |>
  plot(xlab="F(disp)", ylab="F(mpg)", main="Kopula normálneho rozdelenia")
model$cop$rgum$obj |>
  copula::rCopula(n=100) |>
  plot(xlab="F(disp)", ylab="F(mpg)", main="Preklopená Gumbelova kopula")
model$cop$rcla$obj |>
  copula::rCopula(n=100) |>
  plot(xlab="F(disp)", ylab="F(mpg)", main="Preklopená Claytonova kopula")

Nasleduje test dobrej zhody (Goodness-of-fit test). Metóda resamplingu simulation = "mult" je rýchla, ale vhodná skôr pre väčšie súbory údajov. (Pre časovú náročnosť metódy “pb” je výsledok nasledujúceho kódu uvedený iba v komentári.)

Code

# pre jednu kopulu:
# copula::gofCopula(model$cop$rcla$obj, as.matrix(dat), estim.method = "itau", simulation = "pb", ties=TRUE, verbose = T)[c("statistic","p.value")] |> sapply(unname)
# spolu:
mapply(FUN = function(kop, met) 
  copula::gofCopula(model$cop[[kop]]$obj, 
                    as.matrix(dat),
                    N = 10,
                    estim.method = met,
                    ties = TRUE,   #set to FALSE or NA, if no ties are present
                    verbose = FALSE  # whether to show progress bar
                    )[c("statistic","p.value")] |> sapply(unname),  # format
  kop = c("nor", "rgum", "rcla"),
  met = c("mpl", "mpl", "itau")
) |> round(3)
# results for bootstrap method, N=1000 and met = c("mpl", "mpl", "itau"):
#              nor   rgum   rcla
# statistic 0.0391 0.0275 0.0214
# p.value    0.153  0.214  0.504

Všetky kopuly sa zdajú byť vhodné na popis závislostnej štruktúry v náhodnom vektore (disp,mpg). Najlepšou vyzerá byť preklopená Claytonova archimedovská kopula.

Spojenie C(F₁,F₂)

Code

model$mix$obj <- copula::mvdc(
  copula = model$cop$rcla$obj,
  margins = c("lnorm", "norm"),
  paramMargins = list(
    c(meanlog = log(mean(dat$disp)), sdlog=log(sd(dat$disp))), 
    fit_norm["mom",]
    )
)

Takéto zadávanie marg. je trochu nepohodlné, ešteže vytvorenie PDF/CDF z mvdc je ľahšie.

Code

# copula::pMvdc(c(300,20), mvdc = model$mix$obj)   # Error
# copula::pMvdc(cbind(c(300,301),c(20,21)), mvdc = model$mix$obj)   # OK for nrow>1
model$mix$pdf <- function(x) copula::dMvdc(x, model$mix$obj)
model$mix$cdf <- function(x) copula::pMvdc(x, model$mix$obj)

My ho však nepoužijeme (pre problémy s výpočtom jedinej hodnoty i zobrazením hustoty). Namiesto toho, použijeme marginálne PDF/CDF (keď ich už máme explicitne definované) na zostavenie združeného rozdelenia podľa známych vzťahov:
\begin{split} f(x_1,x_2) &= c(F_1(x_1),F_2(x_2)) f_1(x_1) f_2(x_2) \\ F(x_1,x_2) &= C(F_1(x_1),F_2(x_2)) \end{split} Takto dokážeme využiť aj kopulu odhadnutú mimo balíka copula (napr. z balíkov VineCopula, HAC, acopula)

Code

model$mix <- with(model, {  
  F1 <- marg$disp$lnor$cdf
  f1 <- marg$disp$lnor$pdf
  F2 <- marg$mpg$nor$cdf
  f2 <- marg$mpg$nor$pdf
  list(obj = mix$obj,
       pdf = function(x) {
         x <- rbind(x)
         cop$rcla$pdf(cbind(F1(x[,1]), F2(x[,2]))) * f1(x[,1]) * f2(x[,2])
       },
       cdf = function(x) {
         x <- rbind(x)
         cop$rcla$cdf(cbind(F1(x[,1]), F2(x[,2])))
       }
  )
})
# model$mix$cdf(c(300,20))

Kontrolne poskladáme aj združené normálne rozdelenie

Code

model$mix.nor <- with(model, {
  F1 <- marg$disp$nor$cdf
  f1 <- marg$disp$nor$pdf
  F2 <- marg$mpg$nor$cdf
  f2 <- marg$mpg$nor$pdf
  list(
    obj = copula::mvdc(
      copula = cop$nor$obj,
      margins = c("norm", "norm"),
      paramMargins = list(
        c(mean = mean(dat$disp), sd = sd(dat$disp)), 
        fit_norm["mom",]
      )
    ),
    pdf = function(x) {
      x <- rbind(x) |> unname()
      cop$nor$pdf(cbind(F1(x[,1]), F2(x[,2]))) * f1(x[,1]) * f2(x[,2])
    },
    cdf = function(x) {
      x <- rbind(x) |> unname()
      cop$nor$cdf(cbind(F1(x[,1]), F2(x[,2])))
    }
  )
})
# model$mix.nor$cdf(c(300,20))

Zobrazenie hustoty

Empirické rozdelenie

Nespojité empirické rozdelenie

Code

pocetnosti2D <- table(cut(dat$disp, 5), cut(dat$mpg, 5))
plot3D::hist3D(z = pocetnosti2D,
               col = plot3D::ramp.col(c("lightblue","brown")), 
               border = "black", 
               xlab = "disp", ylab = "mpg", zlab = "freq", main = "3D histogram")
plot3D::image2D(z = pocetnosti2D, 
                col = plot3D::ramp.col(c("lightblue","brown")), 
                xlab = "disp", ylab = "mpg", main = "Mozaikový graf")

Jadrovo vyhladené empirické rozdelenie

Code

plot(model$smoo$pdf,,TRUE)
rect(xleft=0,ybottom=0,xright=300,ytop=20, lty="dashed")
plot(model$smoo$pdf,TRUE)

Zložené rozdelenie

Marginálie

Code

plot(model$marg$disp$lnor$pdf, from = 0, to = max(dat$disp), 
     xlab="disp", ylab="hustota")
lines(density(dat$disp), lty="dashed")
legend("topright", c("log-normal","non-param."), lty=c(1,2))
plot(model$marg$mpg$nor$pdf, from = 0, to = max(dat$mpg), 
     xlab="mpg", ylab="")
lines(density(dat$mpg), lty="dashed")
legend("topright", c("normal","non-par."), lty=c(1,2))

Kopula

Code

copula::contourplot2(model$cop$rcla$obj, FUN = copula::dCopula, 
                     region=F, labels=F, n.grid=100, cuts = 100, 
                     main = "hustota kopuly")
copula::wireframe2(model$cop$rcla$obj, FUN = copula::dCopula, region=F,
                   zlab = "hustota")

Kopula + marginálie

Code

#z neznámeho dôvodu kreslí nezmysly:
#copula::contourplot2(model$mix$obj, FUN = copula::dMvdc,  
#                     xlim = c(0,300), ylim = c(0,40),
#                     region=F, labels=F, n.grid=100, cuts = 100)
# graf pomocou základnej funkcie contour()
#local({
#  n <- 100
#  X <- seq_range(dat$disp,length = n)
#  Y <- seq_range(dat$mpg,length = n)
#  Z <- outer(X, Y, FUN = function(x,y) model$mix$pdf(cbind(x,y)))
#  contour(X, Y, Z, drawlabels = FALSE, 
#          xlab = "disp", ylab="mpg", main="hustota rozdelenia")
#})
tmp <- expand.grid(disp = seq_range(dat$disp,length = 100),
           mpg = seq_range(dat$mpg,length = 100)) |> 
  transform(density = model$mix$pdf(cbind(disp,mpg)))
# graf pomocou základného balíku lattice:
# lattice::contourplot(density ~ disp + mpg, tmp, labels=FALSE,
#                     main="hustota rozdelenia (kopula + marginálie)")
ggplot(tmp) + aes(x = disp, y = mpg) + geom_contour(aes(z = density)) +
  ggtitle("hustota rozdelenia (kopula + marginálie)")

Nepodmienené pravdepodobnosti

Pr(X \leq 300, Y \leq 20)

Code

model[c("emp", "norm", "smoo", "mix", "mix.nor")] |> 
  sapply(function(x) x$cdf(c(300,20))) |> 
  round(4)

    emp    norm    smoo     mix mix.nor 
 0.2188  0.2244  0.2617  0.2599  0.2137

Rozdiel medzi “nor” a “mix.nor” (oboje predstavujú združené normálne rozdelenie) je spôsobený odlišnými metódami odhadu.

Pr(0 < X \leq 300, 0 < Y \leq 20)

Code

with(model$emp, {
  x <- c(0,300)
  y <- c(0,20)
  cdf(c(x[2],y[2])) - cdf(c(x[1],y[2])) - cdf(c(x[2],y[1])) + cdf(c(x[1],y[1]))
})

[1] 0.21875

Vrstevnicový graf

Code

expand.grid(disp = seq(0, max(dat$disp)*1.3, length.out = 100),
            mpg = seq(0, max(dat$mpg)*1.3, length.out = 100)
            ) |> 
  transform(CDF = model$mix$cdf(cbind(disp,mpg))) |>
  ggplot() + aes(x = disp, y = mpg) + 
#  metR::geom_contour2(aes(z = CDF, label = stat(level))) +  # s popisom vrstevníc
 geom_contour(aes(z = CDF, colour = after_stat(level))) +
  geom_point(x = 300, y = 20, color = "red") + 
  ggtitle(" CDF rozdelenia (kopula + marginálie)")

Code

# pomocou lattice:
#  lattice::contourplot(x = CDF ~ disp + mpg, at=c(0.1, 0.22, 0.5,0.7,0.8,0.9,0.95),
#                       main=" CDF (kopula + marginálie)")

Podmieňovanie

Pravdepodobnosť Pr(Y < y | X = x) = ?

Jadrovo vyhladené empirické rozdelenie

podmienená CDF:

Code

model$smoo$cdf_disp300 <- dat |> 
  as.matrix() |>
  probhat::cdfc.cks(a = c(0,0), conditions = c(disp=300))  # F(mpg|disp=300)

F(mpg=20|disp=300)

Code

model$smoo$cdf_disp300(20)

[1] 0.8210729

znázornenie:

Code

plot(model$smoo$cdf_disp300, 
     main = "Pr(mpg < x | disp = 300)", xlab="mpg", ylab="CDFc")
abline(v=20, h = model$smoo$cdf_disp300(20), lty="dashed")

Normálne rozdelenie

podmienená PDF

Code

model$norm$pdfc <- function(x, cond) {
  ind <- match(names(cond), names(dat))
  cond <- unname(cond)
  mat <- matrix(nrow=length(x), ncol=2)
  mat[,3-ind] <- x
  mat[,ind] <- cond
  model$norm$pdf(mat)/dnorm(cond, mean=mu[ind], sd=sqrt(Sigma[ind,ind]))
} 
# f(mpg=20|disp=300)
# model$norm$pdfc(20, cond = c(disp = 300))
# pre porovnanie:
# probhat::pdfc.cks(as.matrix(dat), a = c(0,0), conditions = c(disp=300))(20)
data.frame(mpg = seq_range(dat$mpg, 100)) |>
  transform(PDFc = model$norm$pdfc(mpg, cond=c(disp=300))) |>
  plot(type="l", main = "f(mpg | disp = 300)")
abline(v = 20, lty="dashed")

podmienená CDF

Code

model$norm$cdfc <- Vectorize(function(x, cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  lo <- mu[ind] - sqrt(Sigma[ind,ind])*4
  integrate(f = model$norm$pdfc, lower = lo, upper = x, cond = cond)$value  # if lower=-Inf, then unexpected behavior occurs, e.g. model$norm$cdfc(70, cond = c(disp = 300)) gives zero
}, vectorize.args = "x")

F(mpg=20|disp=300)

Code

model$norm$cdfc(20, cond = c(disp = 300))

[1] 0.8063011

Kopula + marginálie

podmienená PDF

Code

model$mix$pdfc <- function(x, cond) {
  condname <- names(cond)
  ind <- match(condname, names(dat))
  cond <- unname(cond)
  mat <- matrix(nrow=length(x), ncol=2)
  mat[,3-ind] <- x
  mat[,ind] <- cond
  model$mix$pdf(mat)/model$marg[[condname]][[c("lnor","nor")[ind]]]$pdf(cond)
}
#model$mix$pdfc(20, cond = c(disp = 300))

podmienená CDF

Code

model$mix$cdfc <- Vectorize(function(x, cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  lo <- mu[ind] - sqrt(Sigma[ind,ind])*4
  integrate(f = model$mix$pdfc, lower = lo, upper = x, cond = cond)$value  # if lower=-Inf, then unexpected behavior occurs, e.g. model$norm$cdfc(70, cond = c(disp = 300)) gives zero
}, vectorize.args = "x")

F(mpg=20|disp=300)

Code

model$mix$cdfc(20, cond = c(disp = 300))

[1] 0.9584723

Kvantil Pr(Y < ? | X = x) = p

Jadrovo vyhladené empirické rozdelenie

podmienená kvantilová funkcia:

Code

model$smoo$icdf_disp300 <- dat |> 
  as.matrix() |>
  probhat::qfc.cks(a = c(0,0), conditions = c(disp=300))  # F(mpg|disp=300)

kvantil

Code

model$smoo$icdf_disp300(0.8210729)

[1] 19.99966

Normálne rozdelenie

podmienená kvantilová funkcia F^{-1}(p|X=x)

Code

model$norm$icdfc <- function(p, cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  int <- mu[ind] + c(-1,1)*sqrt(Sigma[ind,ind])*3
  sapply(p, function(p) uniroot(function(x) model$norm$cdfc(x, cond) - p, interval = int, extendInt = "yes")$root)
}

kontrolný výpočet a podmienený medián

Code

c(0.8063011, 0.5) |> 
  (\(x) data.frame(prob = x, quant = model$norm$icdfc(x, c(disp=300))))()

       prob   quant
1 0.8063011 20.0000
2 0.5000000 17.2353

Kopula + marginálie

Code

model$mix$icdfc <- function(p, cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  int <- mu[ind] + c(-1,1)*sqrt(Sigma[ind,ind])*4
  sapply(p, function(p) uniroot(function(x) model$mix$cdfc(x, cond) - p, interval = int, extendInt = "yes")$root)
}

kontrolný výpočet a podmienené kvantily (medián, hranice predikčného intervalu)

Code

c(0.9584249, 0.5, 0.1, 0.9) |> 
  (\(x) data.frame(prob = x, quant = model$mix$icdfc(x, c(disp=300))))()

       prob    quant
1 0.9584249 19.99879
2 0.5000000 15.77874
3 0.1000000 11.05249
4 0.9000000 18.97367

Stredná hodnota E[Y | X = x]

E[mpg|disp=300]

Jadrovo vyhladené empirické rozdelenie

funkcia

Code

model$smoo$expc <- function(cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  up <- mu[ind] + 4*sqrt(Sigma[ind,ind])
  integrate(f = function(x) x * probhat::pdfc.cks(as.matrix(dat), a = c(0,0), conditions = cond)(x),
            lower = 0, 
            upper = up)$value
}

výpočet

Code

model$smoo$expc(c(disp=300))

[1] 16.4832

Normálne rozdelenie

funkcia

Code

model$norm$expc <- function(cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  int <- mu[ind] + c(-1,1)*sqrt(Sigma[ind,ind])*4
  integrate(f = function(x) x * model$norm$pdfc(x, cond),
            lower = int[1], 
            upper = int[2])$value
}

výpočet

Code

model$norm$expc(c(disp=300))

[1] 17.23532

Code

# kontrola: lm(mpg ~ disp, dat) |> predict(data.frame(disp=300))

Kopula + marginálie

funkcia

Code

model$mix$expc <- function(cond) {
  ind <- 3 - match(names(cond), names(dat))  # 3 = p+1
  int <- mu[ind] + c(-1,1)*sqrt(Sigma[ind,ind])*4
  integrate(f = function(x) x * model$mix$pdfc(x, cond),
            lower = int[1], 
            upper = int[2])$value
}

výpočet

Code

model$mix$expc(c(disp=300))

[1] 15.32598

Zobrazenie podmienenej hustoty

Code

tmp <- c(mean=model$mix$expc(c(disp=300)), 
         model$mix$icdfc(c(med=0.5, q10=0.1, q90=0.9), c(disp=300))
         ) |> 
  round(2) |> print()

 mean   med   q10   q90 
15.33 15.78 11.05 18.97

Code

data.frame(mpg = seq(0, max(dat$mpg), length.out = 100)) |>
  transform(hustota = model$mix$pdfc(mpg, cond = c(disp = 300))) |>
  ggplot() + aes(x = mpg) + geom_line(aes(y = hustota)) +
  scale_x_continuous(breaks = c("0"=0, "10"=10, tmp, "20"=20, "30"=30), 
                     minor_breaks = NULL) +
  theme(axis.text.x = element_text(angle=45)) +
  ggtitle("Podmienená hustota mpg pre disp = 300")

Vektorizovana funkcia podmienenej hustoty, kde prvá súradnica x je podmieňujúca a druhá súradnica je podmieňovaná premenná:

Code

model$mix$pdf2.1 <- function(x) {
  x <- rbind(x)
  model$mix$pdf(x)/model$marg$disp$lnor$pdf(x[,1])
}

Vektorizovaná funkcia podmienenej strednej hodnoty, kde argumentom je samozrejme podmieňujúca premenná.

Code

model$mix$exp2.1 <- Vectorize(function(x) {
  int <- mu[2] + c(-1,1)*sqrt(Sigma[2,2])*4
  integrate(f = function(y) y * model$mix$pdfc(y, cond=c(disp=x)),
            lower = int[1], 
            upper = int[2])$value
})

Podmienená pravdepodobnosť pre zvolené úrovne disp a podmienená stredná hodnota (prerušovaná krivka):

Code

expand.grid(disp = c(100, 200, 300, 400),
            mpg = seq(5, 35, 0.1)) |>
  transform(density = model$mix$pdf2.1(cbind(disp, mpg))) |>
  ggplot() + aes(x = mpg, y = disp) +
  geom_point(data = dat) +
  geom_line(data = data.frame(disp = seq(80,470,10)) |>
              transform(mpg = model$mix$exp2.1(disp)), 
            linetype = 2) +
  ggridges::geom_ridgeline(aes(height = density, group = disp), 
                           stat="identity", scale=400, alpha=0.5) +
  coord_flip() +
  ggtitle("Podmienená stredná hodnota a hustota dojazdu pre zvolené úrovne objemu")

Výstavba modelov združeného rozdelenia pravdepodobnosti a ich aplikácia

Prieskumná analýza

Technický úvod

Jednorozmerná

Závislosť

Spolu

Odhad modelu

Empirické rozdelenie

Združené normálne rozdelenie

Jadrovo vyhladené empirické rozdelenie

Rozkladom na marginálie a kopulu

Marginálne rozdelenia

Neparametrické

Parametrické

Kopula

Neparametrická

Parametrická

Spojenie C(F1,F2)

Zobrazenie hustoty

Empirické rozdelenie

Zložené rozdelenie

Marginálie

Kopula

Kopula + marginálie

Nepodmienené pravdepodobnosti

Podmieňovanie

Pravdepodobnosť Pr(Y < y | X = x) = ?

Jadrovo vyhladené empirické rozdelenie

Normálne rozdelenie

Kopula + marginálie

Kvantil Pr(Y < ? | X = x) = p

Jadrovo vyhladené empirické rozdelenie

Normálne rozdelenie

Kopula + marginálie

Stredná hodnota E[Y | X = x]

Jadrovo vyhladené empirické rozdelenie

Normálne rozdelenie

Kopula + marginálie

Zobrazenie podmienenej hustoty

Spojenie C(F₁,F₂)