Ako je \(Y\) slučajna veličina, poznato je da je \[F_Y(Y) \stackrel{d}{=} U \sim U[0,1],\] pa odatle i \[F_Y^{-1}(U) \stackrel{d}{=} Y,\] pa je zato standaran način generisanja brojeva iz raspodele \(F_Y\) taj da se generiše niz \(u_1, \dots u_n\) iz uniformne \([0,1]\) raspodele, pa željeni uzorak \(y_i\) dobijamo sa \[y_i = F_Y^{-1}(u_i).\]

• Problem koji se često pojavljuje je da ne postoji analitički oblik inverza \(F_Y^{-1}\) i njegovo izračunavanje može biti veoma vremenski zahtevno, pa je ovaj metod generisanja neupotrebljiv za iole veće uzorke, jer je za svaki element uzorka potrebno izračunati vrednost inverza.

• Zato bismo da nadjemo način da sa manjim brojem izračunavanja inverza dobijemo metod za generisanje uzoraka.

• Moguć način da to izvedemo je interpolacija.

• Cilj interpolacije je da aproksimiramo funkciju koju teško izračunavamo nekim polinomom, čije se vrednosti poklapaju sa vrednostima funkcije u odredjenim čvornim tačkama.

• Mi ćemo koristiti samo interpolaciju Lagranževim polinomom, koji je oblika \[F_Y^{-1}(u) \approx \sum_{i=1}^NF_Y^{-1}(x_i)l_i(u),\] gde je \(N\) broj čvornih tačaka \(x_i\), а \(l_i(u) = \prod\limits_{j=1, j\neq i}^N\frac{u-x_j}{x_i-x_j}\).

• Iz izraza aproksimacije \[F_Y^{-1}(u) \approx \sum_{i=1}^NF_Y^{-1}(x_i)l_i(u),\quad l_i(u) = \prod\limits_{j=1, j\neq i}^N\frac{u-x_j}{x_i-x_j},\] vidimo da nam je potrebno samo \(N\) inverzija \(F_Y^{-1}(x_i)\), dok ostatak izraza zavisi samo od čvornih tačaka.

• Kada se ova suma sredi zapravo se dobije polinom po \(x\), pa se zato može jednostavno i veoma brzo izračunavati.

• Razvijen oblik polinoma se direktno može i vektorski izračunavati, što je od značaja kada se radi sa programskim jezicima poput R-a.

• Uzmimo da je \(Y\sim Logistic(0,1)\) i interpolacijom generišimo uzorak. Ako su \(u_1,\dots,u_n\sim U[0,1]\), uzorak \(y_i\) pravimo sa \(y_k = \sum_{i=1}^NF_Y^{-1}(x_i)l_i(u_k)\). Ocena gustine takvog uzorka je:

• 

Problem možemo da vidimo kada pogledamo konkretno kvalitet same interpolacije.

Pri krajevima pravimo velike greške, jer tu \(F_Y^{-1}\) naglo raste, a polinom je mnogo sporiji. Možemo da povećamo broj tačaka, ali ne dobijamo veliko poboljšanje.

• Možemo i da dodamo par tačaka bliže krajevima intervala, da bismo ispratili rast funkcije, ali…

• 

• Da bismo smanjili problem nestabilnosti na krajevima intervala, možda bi pomoglo da nekako razvučemo funkciju \(F_Y^{-1}(u)\) na veći domen, npr. celo \(\mathbb{R}\), da bi polinom manje morao da se adaptira ekstremnim vrednostima.

• To možemo da postignemo nekom smenom oblika \(u=h(x)\), gde je \(h:\mathbb{R}\to[0,1]\).

• Sada interpoliramo funkciju \(g(x) = F_Y^{-1}(h(x)) :\mathbb{R}\to\mathbb{R}\).

• Problem: Kako biramo smenu \(h\)?

• Da bismo koristili \(g(x)\) za generisanje slučajnih veličina, \(h(x)\) mora da ima \(U[0,1]\) raspodelu.

• Rešenje: Prirodno je da \(h(x)\) bude funkcija raspodele neke nove slučajne veličine \(X\), odakle je \(h(X)=F_X(X) \sim U[0,1].\)

• Zato je i \[F_Y^{-1}(F_X(X)) \stackrel{d}{=} Y.\]

• Uzorak iz \(Y\) generišemo sa \[y_i = F_Y^{-1}(F_X(\xi_i)),\] gde je \(\xi_1,\dots,\xi_n\) uzorak iz raspodele \(F_X\).

• Sada nećemo da interpoliramo funkciju \(F_Y^{-1}\), već \(F_Y^{-1}\circ F_X : \mathbb{R} \to \mathbb{R}\), pa dobijamo polinom \[g(x) := F_Y^{-1}(F_X(x)) \approx \sum_{i=1}^NF_Y^{-1}(F_X(x_i))l_i(x) =: g_N(x),\] gde su opet \(l_i(x) = \prod\limits_{j=1, j\neq i}^N\frac{x-x_j}{x_i-x_j}\), a \(x_i\) čvorne tačke.

• Dakle, uzorak iz \(Y\) generišemo sa \[y_i = g_N(\xi_i),\] gde je \(\xi_1,\dots,\xi_n\) uzorak iz \(F_X\).

• Mi ćemo se fokusirati samo na slučaj kada je \(X\sim\mathcal{N}(0,1)\). Ovaj izbor daje dobre rezultate za veliku klasu raspodela za \(Y\).

• Ponovo ćemo pristupiti generisanju logističke raspodele. Funkcija koju aproksimiramo je sada \(g_N(x) \approx F_Y^{-1}(\Phi(x))\).

• Uzorak generišemo tako što generišemo standardan normalan uzorak \(\xi_1,\dots,\xi_n\) i primenimo \(y_i = g_N(\xi_i)\).

• Na sledećim slajdovima vidimo rezultate aproksimacije.

Ovaj put smo koristili samo 5 čvornih tačaka!

Razlike se javljaju oko -5 i 5, ali nisu toliko ekstremne kao kod uniformne.

KS test nam potvrdjuje dobru ocenu (za uzorak obima 100 000):

set.seed(216)
ks.test(pol_norm(rnorm(1e5)), plogis)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pol_norm(rnorm(1e+05))
## D = 0.0026158, p-value = 0.5006
## alternative hypothesis: two-sided

… i gustina izgleda jako dobro

Ovako bi izgledalo sa 7 tačaka…

Razlike sad gotovo i da nema.

KS test, očekivano, opet prolazi.

set.seed(216)
ks.test(pol_norm(rnorm(1e5)), plogis)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pol_norm(rnorm(1e+05))
## D = 0.0016458, p-value = 0.9493
## alternative hypothesis: two-sided

• U prethodnim primerima, čvorne tačke polinoma koje smo uzimali su bile optimalne u nekom smislu. Proći ćemo ukratko motivaciju izbora tačaka.

• Cilj nam je da smanjimo \(L^2\) grešku interpolacije, koja je oblika \[E\left[(g(X) -g_N(X))^2\right]=\int_{-\infty}^{\infty}\left(g(x) - g_N(x)\right)^2f_X(x)\,dx.\]

• Ovo je zapravo integracija \((g(X) -g_N(X))^2\), sa težinskom funkcijom \(f_X(x)\).

• Zato možemo primeniti znanje iz numeričke matematike…

• Ideja je da ovaj integral možemo da predstavimo pomoću Gausove kvadraturne formule \[\int_{-\infty}^{\infty}\left(g(x) - g_N(x)\right)^2f_X(x)\,dx =\\ \sum_{i=1}^N(g(x_i) -g_N(x_i))^2w_i + \varepsilon_N.\]

• Stoga možemo da koristimo i teoriju vezanu za Gausove kvadraturne formule.

• Pokazuje se da optimalne tačke \(x_i\) i težine \(w_i\) možemo izračunati koristeći samo prvih \(2N\) momenata raspodele \(X\).

• Kako je \(g_N\) interpolacioni polinom za \(g\), u čvornim tačkama oni se poklapaju pa je \[\sum_{i=1}^N(g(x_i) -g_N(x_i))^2w_i + \varepsilon_N = \varepsilon_N.\]

• Stoga se \(L^2\) greška interpolacije svodi samo na grešku integracije, koju smo smanjili izborom optimalnih tačaka \(x_i\) i težina \(w_i\).

Za \(X\sim\mathcal{N}(0,1)\), izvodi se da je ova greška jednaka:

\[\varepsilon_N = \frac{N!\sqrt\pi}{2^N}\frac{\Psi^{(2N)}(\hat\xi)}{2N},\] gde je \(\Psi(x) = (g(x)-g_N(x))^2 = \left(\frac1{N!}\left.\frac{d^Ng(x)}{dx^N}\right|_{x=\hat\xi}\prod_{i=1}^N(x-x_i)\right)^2\) kvadrat standardne greške Lagranževe interpolacije.

Za dovoljno glatke \(\Psi\), greška teži nuli, kad \(N\to\infty\).

• Pokrenućemo u nastavku još nekoliko primera generisanja slučajnih veličina pomoću SCMC metode.

• Cilj ovih primera je da demonstriraju brzinu koja se dobija SCMC metodom.
  Zato su prikazani samo rezultati merenja vremena izvršavanja.

• Generisanje vršimo pomoću paketa scmc, koji razvijamo, a razvojna verzija se može instalirati sa

  devtools::install_github("blaza/scmc")

• Odabir broja čvorova je takav da ks.test potvrdi da je raspodela tačna.

library(scmc)
scmc_gamma <- univariate_sampler(function(p) qgamma(p, 5, 2),
                                 optimal_points(5, norm_moment_fun))
scmc_gamma_c <- univariate_sampler(function(p) qgamma(p, 5, 2),
                                   optimal_points(5, norm_moment_fun),
                                   compile_c = TRUE)
microbenchmark(rnorm = rnorm(1e5), scmc = abs(scmc_gamma(1e5)),
               scmc_C = abs(scmc_gamma_c(1e5)), 
               rgamma = rgamma(1e5, 5, 2))

## Unit: milliseconds
##    expr       min        lq      mean    median        uq        max neval
##   rnorm 13.515143 13.635918 14.999591 14.263949 16.268135  20.423739   100
##    scmc  4.694448  4.922791  6.020736  5.303754  7.029676   9.632422   100
##  scmc_C  3.067484  3.185539  6.983829  3.555152  6.727582 138.167565   100
##  rgamma 18.434582 18.603646 20.268926 18.766650 22.229071  26.323712   100

library(scmc)
scmc_chsq <- univariate_sampler(function(p) qchisq(p, 4, 1),
                                optimal_points(9, norm_moment_fun))
scmc_chsq_c <- univariate_sampler(function(p) qchisq(p, 4, 1),
                                  optimal_points(9, norm_moment_fun),
                                  compile_c = TRUE)
microbenchmark(rnorm = rnorm(1e5), scmc = abs(scmc_chsq(1e5)),
               scmc_C = abs(scmc_chsq_c(1e5)), 
               rchisq = rchisq(1e5, 4, 1))

## Unit: milliseconds
##    expr       min        lq      mean    median       uq       max neval
##   rnorm 13.558461 13.746552 15.773072 14.431104 17.34818  23.08362   100
##    scmc  6.854567  7.154001  8.839418  8.371104 10.58059  13.70566   100
##  scmc_C  3.241296  3.380656  5.982371  3.829946  6.05731 129.76888   100
##  rchisq 41.607670 42.031585 45.465797 42.436744 49.12297  59.42664   100

library(scmc)
scmc_t <- univariate_sampler(function(p) qt(p, df = 4),
                             optimal_points(11, norm_moment_fun))
scmc_t_c <- univariate_sampler(function(p) qt(p, df = 4),
                               optimal_points(11, norm_moment_fun),
                               compile_c = TRUE)
microbenchmark(rnorm = rnorm(1e5), scmc = scmc_t(1e5),
               scmc_C = scmc_t_c(1e5), 
               rt = rt(1e5, df = 4))

## Unit: milliseconds
##    expr       min        lq      mean    median        uq        max neval
##   rnorm 13.626112 13.755363 15.178586 14.272531 15.817555  23.339215   100
##    scmc  7.629217  8.234433 11.063055  9.034328 10.867329 139.061500   100
##  scmc_C  2.990105  3.175205  4.044539  3.571476  4.350525   8.520896   100
##      rt 39.430932 40.006968 42.488959 40.585304 43.985382  57.383725   100

• Primetimo da, ako bismo interpolirali \(F_{Z|Y=y_i}^{-1}(\Phi(\xi_i))\), morali bismo da konstruišemo zaseban polinom za svako \(y_i\), što je potpuno neupotrebljivo.

• To možemo da izbegnemo koristeći dvodimenzionu Lagranževu interpolaciju.

• Posmatramo funkciju \(g(x, y) := F_{Z|Y=y}^{-1}(\Phi(x))\), i nju interpoliramo:

• \[g(x, y) \approx \sum_{i=1}^{N_X}\sum_{j=1}^{N_Y}F_{Z|Y=y_j}^{-1}(\Phi(x_i))l_Y^{(j)}(y)l_X^{(i)}(x) =: g_N(x,y),\] gde su \(l_X^{(i)}(x)=\prod\limits_{j=1, j\neq i}^{N_X}\frac{x-x_j}{x_i-x_j}\) i \(l_Y^{(j)}(y)=\prod\limits_{i=1, j\neq i}^{N_Y}\frac{y-y_i}{y_j-y_i}\), a \(y_1,\dots, y_{N_Y}\) čvorne tačke za interpolaciju za promenljivu \(y\).

• Ovime smo pokrili ceo domen \((x,y)\) dvodimenzionim polinomom, pa \(Z\) sada generišemo na isti način kako smo i zamislili, samo što umesto \(F_{Z|Y=y_i}^{-1}(\Phi(\xi_i))\) koristimo aproksimaciju koju smo naveli, pa dobijamo \[z_i = g_N(\xi_i, y_i) \approx F_{Z|Y=y_i}^{-1}(\Phi(\xi_i)).\]

• Generisanje se sada svodi na generisanje \(\xi_i\), \(y_i\) i evaluaciju dvodimenzionog polinoma \(g_n(x,y)\), što se sve radi vrlo brzo.

• Optimalne tačke za \(Y\) možemo računati na osnovu uzoračkih momenata generisanog uzorka \(y_i\), dok za \(X\) koristimo isti metod kao i za jednodimenzione raspodele.

• Ovaj postupak se može lako uopštiti na veličine \(Z|X_1,\dots,X_d\), gde samo koristimo \((d+1)\)-dimenzioni polinom.

• Generisanje uslovnih raspodela je od velike koristi za generisanje višedimenzionih raspodela.

• Recimo da hoćemo da generišemo vektor \((X, Y)\).
  To možemo da uradimo u dva koraka:
  1. Generišemo uzorak \(x_1,\dots,x_n\) iz raspodele \(X\)
  2. Generišemo uzorak \(y_1,\dots,y_n\) iz \(Y\) koristeći uslovne raspodele \(Y|X=x_i\)

• Vidimo da nam je za generisanje potrebno znanje marginalne raspodele \(X\) i uslovne raspodele \(Y|X=x_i\).

• Svaki korak se može odraditi SCMC metodom

• Postupak se lako uopštava na više od dve dimenzije.

• Kad generišemo \((X_1,\dots,X_d)\), imamo \(d\) koraka:
  • Generišemo \(X_1\)
  • Generišemo \(X_2|X_1\)
  • Generišemo \(X_3|X_2, X_1\)
    \(\vdots\)
  • Generišemo \(X_d|X_{d-1},\dots,X_2, X_1\)

• U svakom koraku \(k\) imamo \(k\)-dimenzioni interpolacioni polinom.

• Neka je \((X_1,X_2) \sim \mathcal{N}\left(\begin{pmatrix}1 \\ 2\end{pmatrix}, \begin{pmatrix}1 & 0.3\\ 0.3 & 1\end{pmatrix}\right)\).

• Znamo da je onda \(X_1 \sim \mathcal{N}(1, 1)\).

• Takodje, \(X_2|X_1=x \sim \mathcal{N}\left(\mu_2 + \frac{\sigma_2}{\sigma_1} \rho (x-\mu_1), (1-\rho^2)\sigma_2^2\right)\),
  
  odnosno \(X_2|X_1=x \sim \mathcal{N}\left(1.7 + 0.3x, 0.91\right)\)

• Imamo sve što je potrebno za generisanje.

library(scmc)
gen_x1 <- univariate_sampler(function(p) qnorm(p, mean = 1, sd = 1),
                             optimal_points(3, norm_moment_fun))
gen_x2 <- conditional_sampler(function(x, p) qnorm(p, mean = 1.7 + 0.3*x, sd = 0.91),
                              list(optimal_points(3, sample = gen_x1(1e5)),
                                   optimal_points(3, norm_moment_fun)))

scmc_rmvnorm <- function(n) {
  x1 <- gen_x1(n)
  x2 <- gen_x2(n, list(x1))
  cbind(x1,x2)
}

Uporedimo za početak brzinu sa paketom mvtnorm

microbenchmark(mvtnorm = mvtnorm::rmvnorm(1e5, mean = c(1, 2),
                                sigma = matrix(c(1,0.3,0.3,1), ncol=2)),
               scmc = scmc_rmvnorm(1e5))

## Unit: milliseconds
##     expr      min       lq     mean   median       uq      max neval
##  mvtnorm 36.20379 44.31492 58.13147 46.29024 52.22350 205.4573   100
##     scmc 23.02030 26.90734 32.36595 29.59714 31.21924 165.3798   100

Ocena gustine raspodele jednog uzorka generisanog pomoću SCMC…

Uradimo još i jedan statistički test da potvrdimo da je raspodela korektna…

library(ks)
kde.test(scmc_rmvnorm(1e4), mvtnorm::rmvnorm(1e4, mean = c(1, 2),
                                sigma = matrix(c(1,0.3,0.3,1), ncol=2)))$pvalue

## [1] 0

P-vrednost je velika, stoga ne odbacujemo pretpostavku o jednakosti raspodele našeg uzorka i onog generisanog sa mvtnorm paketom.