Elemente de Simulare în R

Scopul acestor note de laborator este de a prezenta succint metodele de simulare introduse la curs și de a le ilustra printr-o serie de exemple în R.

Elemente de simulare

În acest capitol vom prezenta o serie de exemple de generare a observațiilor dintr-o repartiție dată, folosind două metode generale de generare: metoda inversă și metoda respingerii. Înainte de a prezenta metodele de simulare vom ilustra cum putem estima valorile constantelor \(\pi\) și \(e\) aplicând Legea Numerelor Mari (a se vedea Legea Numerelor Mari).

Estimarea constantelor \(\pi\) și \(e\)

În exercițiul de mai jos sunt prezentate două metode de estimare a lui \(\pi\) prin simulare.

Exercițiul 1 Fie \(U_{i1}, U_{i2}, V_{i}\), \(i\in\{1,2,\dots,n\}\), variabile aleatoare independente repartizate uniform \(\mathcal{U}([0,1])\). Definim variabilele aleatoare

\[ X_i=\left\{\begin{array}{ll} 1, & U_{i1}^2+U_{i2}^2<1\\ 0, & \text{altfel} \end{array}\right.\;\;\;\text{și}\;\;\; Y_i=\sqrt{1-V_i^2},\;\; i\in\{1,2,\dots,n\} \]

Considerăm variabilele aleatoare \(\hat{\pi}_1=\frac{4}{n}\sum_{i=1}^{n}X_i\) și \(\hat{\pi}_2=\frac{4}{n}\sum_{i=1}^{n}Y_i\). Calculați media și varianța acestor variabile și stabiliți care este mai eficientă¹ în estimarea lui \(\pi\).

Soluție. Pentru \(\hat{\pi}_1\): observăm că v.a. \(X_i\) sunt v.a. de tip Bernoulli cu

\[ \begin{aligned} \mathbb{P}(X_i=1)&=\mathbb{P}\left(U_{i1}^2+U_{i2}^2<1\right)=\iint_{\{u^2+v^2<1\}\cap[0,1]^2}f_{(U_{i1},U_{i2})}(u,v)\,dudv \\ &\overset{indep.}{=}\iint_{\{u^2+v^2<1\}\cap[0,1]^2}f_{U_{i1}}(u)f_{U_{i2}}(v)\,dudv\\ &=\int_{0}^{1}\int_{0}^{\sqrt{1-u^2}}1\,dvdu = \int_{0}^{1}\sqrt{1-u^2}\,du\\ &\overset{u=\sin\alpha}{=}\int_{0}^{\frac{\pi}{2}}\cos^2\alpha\,d\alpha=\int_{0}^{\frac{\pi}{2}}\frac{\cos2\alpha+1}{2}\,d\alpha\\ &=\frac{\pi}{4}+\frac{1}{2}\int_{0}^{\frac{\pi}{2}}\cos2\alpha\,d\alpha=\frac{\pi}{4} \end{aligned} \]

O altă variantă de calcul pentru \(\mathbb{P}(X_i=1)\) era să observam că această probabilitate se exprima și ca raportul dintre aria mulțimii \(\{(u,v)\in[0,1]^2\,|,u^2+u^2<1\}\) și cea a pătratului \([0,1]^2\), deci tot \(\frac{\pi}{4}\).

Dacă \(T=\sum_{i=1}^{n}X_i\) atunci \(T\sim\mathcal{B}\left(n,\frac{\pi}{4}\right)\) de unde avem că media este \(\mathbb{E}[T]=\frac{n\pi}{4}\) iar varianța

\[ Var[T]=n\frac{\pi}{4}\left(1-\frac{\pi}{4}\right). \]

Cum \(\hat{\pi}_1=\frac{4}{n}T\) deducem că \(Var[\hat{\pi}_1]=\frac{4\pi-\pi^2}{n}\). Din Legea Numerelor Mari obținem că \(\hat{\pi}_1=\frac{4}{n}\sum_{i=1}^{n}X_i\overset{a.s.}{\to}4\mathbb{E}[X_1]=4\mathbb{P}(X_1=1)=\pi\).

Pentru \(\hat{\pi}_2\), să observăm pentru început că media lui \(Y_1\) este

\[ \mathbb{E}[Y_1] = \int_{0}^{1}\sqrt{1-u^2}\,du = \frac{\pi}{4} \]

iar varianța lui \(Y_1\) este

\[ Var[Y_1] = \mathbb{E}[Y_1^2]-\mathbb{E}^2[Y_1] = \int_{0}^{1}1-u^2\,du -\frac{\pi^2}{16} = \frac{2}{3}-\frac{\pi^2}{16}. \]

Prin aplicarea Legii Numerelor Mari rezultă că

\[ \hat{\pi}_2=\frac{4}{n}\sum_{i=1}^{n}Y_i\overset{a.s.}{\to}4\mathbb{E}[Y_1]=\pi \]

iar varianța lui \(\hat{\pi}_2\) este

\[ Var[\hat{\pi}_2] = \frac{16}{n^2}\sum_{i=1}^{n}Var[Y_i] = \frac{16}{n}\left(\frac{2}{3}-\frac{\pi^2}{16}\right). \]

Pentru a vedea care dintre cei doi estimatori este mai eficient trebuie să verificăm care are varianța mai mică. Cum \(\frac{32}{3}<12<4\pi\) rezultă că \(Var[\hat{\pi}_2]<Var[\hat{\pi}_1]\) deci al doilea estimator este mai eficient.

Figura 1 de mai jos ilustrează acest fenomen:

l <- seq(10,10000,10)
ll <- length(l)

dat <- data.frame(a = 1:length(l), pi1 = 1:length(l), pi2 = 1:length(l))

for (i in 1:length(l)){
  u1 <- runif(l[i])
  u2 <- runif(l[i])
  v <- runif(l[i])
  
  x <- ifelse(sqrt(u1^2+u2^2)<1,1,0)
  y <- sqrt(1-v^2)
  
  dat$pi1[i] <- 4*sum(x)/l[i]
  dat$pi2[i] <- 4*sum(y)/l[i]
}

plot(dat$a, dat$pi1, 
     col = myblue, 
     xlab = "Volumul esantionului",
     ylab = TeX("$\\hat{\\pi}_1$ (albastru) si $\\hat{\\pi}_2$ (rosu)"),
     bty = "n",
     pch = 16)
points(dat$a, dat$pi2,
       col = mygreen,
       pch = 16)
abline(h = pi, 
       lty = 2,
       lwd = 3,
       col = myred)

Figura 1: Ilustrarea celor doi estimatori în raport cu volumul eșantionului.

Următorul exercițiu prezintă o metodă de estimare a lui \(e\) prin simulare și are la bază articolul (Russell 1991).

Exercițiul 2 Fie \(U_1, U_2,\dots, U_n\) variabile aleatoare independente și repartizate \(\mathcal{U}([0,1])\) și \(S_n=\sum_{i=1}^{n}U_i\). Dacă variabila aleatoare \(N\) este definită prin

\[ N = \min\{k\,|\,S_k>1\} \]

atunci:

1. Arătați că dacă \(0\leq t\leq 1\) atunci \(\mathbb{P}(S_k\leq t) = \frac{t^k}{k!}\).

2. Determinați \(\mathbb{E}[N]\) și \(Var[N]\).

Soluție. Avem:

1. Pentru a calcula probabilitatea \(\mathbb{P}(S_k\leq t)\) cu \(0<t<1\) să ne reamintim că dacă \(X\) și \(Y\) sunt două variabile aleatoare independente cu densitățile \(f_X\) și \(f_Y\) atunci densitatea sumei \(Z = X+Y\) (convoluția) este dată de

\[ f_Z(z) = \int f_X(z-t)f_Y(t)\, dt. \]

Fie \(f_n\) densitatea variabilei aleatoare \(S_n\) pentru \(n\geq1\). Avem, pentru \(0<x<1\), că \(f_1(x) = 1\) și pentru a calcula densitatea \(f_{n+1}\) a variabilei aleatoare \(S_{n+1}\) să observăm că \(S_{n+1} = S_{n} + U_{n+1}\) cu \(S_n\) și \(U_{n+1}\) variabile aleatoare independente, de unde aplicând formula pentru densitatea sumei deducem că

\[ f_{n+1}(x) = \int_{0}^xf_n(t)\,dt, \quad n\geq 1. \]

Prin inducție rezultă că \(f_n(x) = \frac{x^{n-1}}{(n-1)!}\) pentru \(0<x<1\) de unde

\[ \mathbb{P}(S_n\leq t) = \int_{0}^{t}f_n(x)\,dx = \int_{0}^{t}\frac{x^{n-1}}{(n-1)!}\,dx = \frac{t^{n}}{n!}. \]

2. Pentru \(n\geq 2\) să observăm că \(\mathbb{P}(N = n) = \mathbb{P}(S_{n-1}<1\leq S_n)\) de unde

\[ \mathbb{P}(N = n) = \mathbb{P}(S_{n-1}<1) - \mathbb{P}(S_n<1) = \frac{1}{(n-1)!} - \frac{1}{n!} = \frac{n-1}{n!}. \]

Pentru medie avem

\[ \mathbb{E}[N] = \sum_{n = 2}^{\infty}n\mathbb{P}(N = n) = \sum_{n = 2}^{\infty}n\frac{n-1}{n!} = \sum_{n = 2}^{\infty}\frac{1}{(n-2)!} = e. \]

În mod similar se poate arăta că \(Var[N] = e(3-e)\).²

l <- seq(10,10000,10)
ll <- length(l)

dat <- data.frame(a = 1:length(l), e1 = 1:length(l))

Sime <- function(n){
  # genereaza observatii din repartitia N
  y <- rep(0, n)
  
  for (i in 1:n){
    k <- 1
    u <- runif(1)
    
    while(u<1){
      u <- u + runif(1)
      k <- k+1
    }
    
    y[i] <- k
  }
  
  return(y)
}

Simev <- Vectorize(Sime, vectorize.args = "n") # vectorize

for (i in 1:length(l)){
  x <- Simev(l[i])
  
  dat$e1[i] <- sum(x)/l[i]
}

plot(dat$a, dat$e1, 
     col = myblue, 
     xlab = "Volumul esantionului",
     ylab = TeX("$N$"),
     pch = 16,
     bty = "n")
abline(h = exp(1), 
       lty = 2,
       lwd = 3,
       col = myred)

Figura 2: Ilustrarea convergenței lui \(N\) către \(e\).

Generarea variabilelor aleatoare prin metoda inversă

Metoda de simulare prezentată în această secțiune se bazează pe Teorema de universalitate a repartiției uniforme. În esență, rezultatul afirmă că putem genera (cel puțin din punct de vedere teoretic) o observație dintr-o repartiție dată plecând de la o observație uniformă pe \([0,1]\). În practică, acest rezultat este ușor de aplicat atunci când putem determina funcția cuantilă sub formă compactă.

Teorema 1 (Universalitatea Repartiției Uniforme) Fie \(X\) o variabilă aleatoare reală cu funcția de repartiție \(F\), \(U\) o variabilă aleatoare repartizată uniform pe \([0,1]\) și fie funcția cuantilă (inversa generalizată) asociată lui \(F\), \(F^{-1}:(0,1)\to\mathbb{R}\) definită prin

\[ F^{-1}(u) = \inf\{x\in\mathbb{R}\,|\,F(x)\geq u\}, \quad \forall u\in(0,1). \]

Atunci

1. Variabilele aleatoare \(X\) și \(F^{-1}(U)\) sunt repartizate la fel.
2. Dacă \(F\) este continuă atunci variabila aleatoare \(F(X)\) este repartizată uniform \(\mathcal{U}(0,1)\).

Generarea variabilelor aleatoare discrete

În cazul în care \(X\) este o variabilă aleatoare discretă ce ia valori într-o mulțime finită \(X(\Omega)=\{x_1,x_2,\ldots,x_n\}\), unde \(x_1\leq\cdots\leq x_n\), și are funcția de masă \(p_k=\mathbb{P}(X = x_k)\), \(k\in\{1,\ldots,n\}\), putem aplica metoda inversă.

Funcția de repartiție este

\[ F(x) = \left\{\begin{array}{lllll} 0 &,\, x<x_1\\ p_1 &,\, x_1\leq x < x_2\\ p_1 + p_2 &,\, x_2\leq x < x_3\\ & \vdots\\ p_1 + \cdots + p_{n-1} &,\, x_{n-1}\leq x < x_n\\ 1 &,\, x\geq x_n \end{array}\right. = \sum_{i = 1}^{n}(p_1+\cdots+p_{i-1})\mathbf{1}_{[x_{i-1},x_i)}(x), \]

unde \(x_0 = -\infty\) iar \(p_0 = 0\), de unde deducem că funcția cuantilă este

\[ F^{-1}(u) = \left\{\begin{array}{lllll} x_1 &,\, 0 < u \leq p_1\\ x_2 &,\, p_1 < u \leq p_1 + p_2\\ & \vdots\\ x_n &,\, p_1+\cdots+p_{n-1}\leq x \leq 1 \end{array}\right. = \sum_{i=1}^{n}x_i\mathbf{1}_{(p_1+\cdots+p_{i-1}, p_1+\cdots+p_i]}(u). \]

În acest caz putem genera o observație cu repartiția lui \(X\) (\(\mathbb{P}\circ X^{-1}\)) plecând de la \(U(\omega)\), o observație \(\mathcal{U}(0,1)\), și considerând \(Y(\omega) = F^{-1}(U(\omega))\). În practică, este suficient să determinăm indicele unic \(i\) pentru care \(p_1+\cdots+p_{i-1}< U(\omega) \leq p_1+\cdots+p_i\) și să luăm \(Y(\omega) = x_i\).

Trebuie menționat că timpul necesar pentru generarea unei observații este proporțional cu numărul de intervale de tipul \((F(x_{i-1}), F(x_i))]\) pe care le căutăm prin urmare este indicat să ordonăm valorile variabilei aleatoare \(X\) după ordinea descrescătoare a funcției de masă, i.e. după \(p_i\).

În cazul în care \(X(\Omega)\) este numărabilă, rezultatele rămân adevărate numai că trebuie să schimbăm sumele finite cu serii.

Exercițiul 3 Fie \(X\) o variabilă aleatoare cu valori în mulțimea \(\{1,2,3,4\}\). Repartiția \(\nu\) a lui \(X\) este

\[ \mathbb{P}(X = 1) = 0.2, \quad \mathbb{P}(X = 2) = 0.5, \quad \mathbb{P}(X = 3) = 0.1, \quad \mathbb{P}(X = 4) = 0.2. \]

1. Simulați un eșantion \(u\) de volum \(10000\) de variabile aleatoare i.i.d. repartizate \(\mathcal{U}([0,1])\).

2. Plecând de la acest eșantion, construiți un eșantion de volum \(10000\), \(x\), de variabile aleatoare i.i.d. repartizate conform \(\nu\).

3. Comparați, cu ajutorul diagramei cu bare verticale (barplot), repartiția eșantionului \(x\) și cea teoretică \(\nu\).

Soluție. Avem:

1. Pentru a genera observații independente repartizate uniform vom folosi funcția runif:

n <- 10000
u <- runif(n)

2. În primul rând să observăm că pentru a genera o variabilă Bernoulli \(\mathcal{B}(1-p)\) este suficient să considerăm variabila \(Y=\mathbf{1}_{\{U>p\}}\) unde \(U\sim\mathcal{U}(0,1)\) care în R se scrie as.integer(runif(1)>p) (folosim funcția as.integer pentru a transforma un vector logic, runif(1)>p, într-un vector de numere întregi).

Astfel, pentru variabila noastră \(X\) putem scrie

x <- 1 + (u > 0.2) + (u > 0.7) + (u > 0.8)

Trebuie remarcat că funcția sample din R permite simularea repartițiilor discrete (e.g. repartiția uniformă pe o mulțime discretă). Pentru exemplul nostru putem să extragem, cu întoarcere, \(n\) numere din mulțimea \(\{1,2,3,4\}\) urmând un vector de probabilități prob:

x <- sample(1:4, n, replace = TRUE, 
           prob = c(0.2, 0.5, 0.1, 0.2))

3. Pentru început trebuie să determinăm câte observații sunt din fiecare valoare unică a lui \(x\). Putem face acest lucru folosind funcția table:

Figura 3: Compararea repartiției teoretice cu cea empirică pentru \(n=10000\) de observații.

Observăm că eșantionul este repartizat conform \(\nu\).

Următorul exercițiu generalizează rezultatul anterior.

Exercițiul 4 Definiți o funcție care să genereze un eșantion de volum \(n\) dintr-o repartiție discretă definită pe mulțimea \(x=\{x_1,\dots,x_N\}\) (vector numeric sau de caractere) cu probabilitatea \(p=\{p_1,\dots,p_N\}\) pe \(x\).

Soluție. Din punct de vedere practic, pentru a implementa metoda inversă va trebui să reordonăm valorile pe care le poate lua \(X\) în așa fel încât vectorul de probabilități \(p\) să fie descrescător.

Avem următoarea funcție:

GenerateDiscrete <- function(n = 1, x, p, err = 1e-15){
  # talia esantionului
  # x alfabetul 
  # p probabilitatile
  lp <- length(p)
  lx <- length(x)
  
  if(abs(sum(p)-1)>err | sum(p>=0)!=lp){
    
    stop("suma probabilitatilor nu este 1 sau probabilitatile sunt mai mici decat 0")
    
  }else if(lx!=lp){
    
    stop("x si p ar trebui sa aiba aceeasi marime")
    
  }else{
    out <- rep(0, n)
    
    indOrderProb <- order(p, decreasing = TRUE) # index
    pOrdered <- p[indOrderProb] # rearanjarea probabilitatilor 
    xOrdered <- x[indOrderProb] # rearanjarea valorilor lui x
    
    u <- runif(n) # generam n uniforme
    pOrderedCS <- cumsum(pOrdered)
    
    for (i in 1:n){
      # determinam indicele
      k <- min(which(u[i]<=pOrderedCS))
      out[i] <- xOrdered[k]
    }
  }
  
  return(out)
}

Pentru a testa această funcție să considerăm următoarele două exemple:

1. Ne propunem să generăm observații din \(X\sim\begin{pmatrix}1 & 2 & 3\\ 0.2 & 0.3 & 0.5\end{pmatrix}\), în acest caz: \(x=[1,2,3]\) și \(p=[0.2,0.3,0.5]\). Începem prin generarea a \(n = 10\) observații din repartiția lui \(X\)

GenerateDiscrete(10, c(1,2,3), c(0.2,0.3,0.5))

 [1] 3 3 3 1 3 3 3 3 2 2

Plecând de la un eșantion de \(n = 10000\) de observații vrem să comparăm, cu ajutorul diagramei cu bare verticale (barplot), repartiția eșantionului cu cea teoretică. Pentru aceasta considerăm următoarea funcție generică:

Compar_diag_bare <- function(pX, pT, legend.loc = "topright", 
                             inset.legend = 0, cex.legend = 1, ...){
  indX <- names(pX)
  
  barplot(rbind(pX, pT), 
        beside = TRUE, 
        space = c(0.1, 1),
        col = c(myblue, myred),
        names.arg = indX, 
        ylab = "Functia de masa",
        cex.axis = 0.7,
        legend.text = c("Repartitia Esantionului", "Repartitia teoretica"),
        args.legend = list(x = legend.loc, bty = "n", 
                           inset = inset.legend, cex = cex.legend), ...)
}

Pentru un eșantion de \(n = 10000\) de observații avem:

Figura 4: Compararea repartiției teoretice cu cea empirică pentru \(n=10000\) de observații.

2. În acest caz considerăm variabila aleatoare \(X\sim\begin{pmatrix}a & b & c & d\\ 0.15 & 0.25 & 0.15 & 0.45\end{pmatrix}\), deci \(x=[a,b,c,d]\) și \(p=[0.15,0.25,0.15,0.45]\). Mai jos generăm \(n = 15\) observații din repartiția variabilei aleatoare \(X\):

GenerateDiscrete(15, c('a','b','c','d'), c(0.15,0.25,0.15,0.45))

 [1] "d" "b" "d" "b" "d" "d" "d" "b" "c" "a" "a" "d" "d" "d" "d"

Ca și în cazul primului exemplu, vom compara repartiția teoretică cu cea a unui eșantion de \(n = 10000\) de observații:

Figura 5: Compararea repartiției teoretice cu cea empirică pentru \(n=10000\) de observații.

În subsecțiunile de mai jos vom testa funcția și pentru o parte din repartițiile discrete de bază: uniformă (discretă), binomială, geometrică, hipergeometrică sau Poisson.

Repartiția Uniformă

În cazul repartiției uniforme pe \(\{1, 2, \ldots, N\}\) (ilustrat mai jos pentru \(N = 7\)) funcția generică conduce la:

Figura 6: Compararea repartiției teoretice \(\mathcal{U}(\{1,2,\ldots,7\})\) cu cea empirică pentru un eșantion de volum \(n=10000\).

Pentru acest caz putem să particularizăm și să găsim un algoritm mai eficient. Ținând cont că funcția de repartiție pentru uniforma pe \(\{1, 2, \ldots, N\}\) este \(F(i) = \frac{i}{N}\) și folosind Teorema de universalitate deducem că

\[ F(i-1) < u \leq F(i) \iff i-1< Nu \leq i \]

cu alte cuvinte \(X\) ia valoarea \(i\) dacă \(i-1< Nu \leq i\) sau altfel spus, dat fiind că \(i\) este întreg, \(i = [Nu] + 1\) adică

\[ X = [Nu] + 1 \]

unde \([x]\) este partea întreagă a lui \(x\). Următoarea funcție implementează această relație:

sim_discrete_unif <- function(n = 100, N = 7){
     # n - nr de observatii
     # N - uniforme pe 1, 2, ..., N
     
     u <- runif(n)
     x <- floor(N*u) + 1
     
     return(x)
}

În R putem utiliza funcția sample(1:N, n, replace = TRUE) pentru a genera un eșantion de volum \(n\) din repartiția uniformă pe \(\{1, 2, \ldots, N\}\).

Repartiția Binomială

Vom aplica funcția generică pentru \(\mathcal{B}(10, 0.3)\):

Figura 7: Compararea repartiției teoretice \(\mathcal{B}(10, 0.3)\) cu cea empirică pentru un eșantion de volum \(n=10000\).

O altă variantă de simulare are la bază calculul funcției de repartiție a \(\mathcal{B}(n,p)\). Ținând cont de relația de recurență

\[ \mathbb{P}(X = i+1) = f(i+1) = \frac{n-i}{i+1}\frac{p}{1-p} f(i), \quad i\geq 1 \]

dintre funcția de masă a \(\mathcal{B}(n,p)\) în \(i+1\) și \(i\), unde \(f(0) = (1-p)^n\) atunci pentru \(u\) o observație uniformă pe \([0,1]\), o observație \(i\) din repartiția \(\mathcal{B}(n,p)\) verifică

\[ \mathbb{P}(X\leq i-1) = F(i-1) < u\leq F(i) = \mathbb{P}(X\leq i) \]

ceea ce conduce la algoritmul

sim_discrete_binom <- function(n = 100, size = 10, p = 0.5){
     # n - nr de observatii
     # size, p - parametrii binomialei
     
     x <- rep(0, n)
     r <- p/(1-p)
     
     for (j in 1:n){
          u <- runif(1)
          x[j] <- 0
          
          pb <- (1-p)^size
          Fr <- pb
          
          while(u >= Fr){
               pb <- (size - x[j])/(x[j] + 1) * r * pb
               Fr <- Fr + pb
               x[j] <- x[j] + 1
          }
     }
     
     return(x)
}

În R putem folosi funcția rbinom() pentru a genera observații repartizate binomial.

Repartiția Geometrică

Aplicăm funcția GenerateDiscrete pentru \(\mathrm{Geom}(0.3)\) și obținem:

Figura 8: Compararea repartiției teoretice \(\mathrm{Geom}(0.3)\) cu cea empirică pentru un eșantion de volum \(n=10000\).

Pentru a îmbunătății algoritmul putem să observăm că în cazul repartiției \(\mathrm{Geom}(p)\) funcția de masă este

\[ \mathbb{P}(X = i) = f(i) = p(1-p)^{i-1}, \quad i\in\{1, 2, \ldots\} \]

ceea ce implică că funcția de repartiție este

\[ F(i) = \sum_{j = 1}^{i}f(j) = 1 - (1-p)^{i-1}. \]

Dacă \(U\sim\mathcal{U}([0,1])\) atunci \(X\) ia valoarea \(i\) care verifică relația \(F(i-1)\leq U< F(i)\) de unde

\[ 1 - (1-p)^{i-1} < U \leq 1 - (1-p)^{i} \iff i< \frac{\log(1-U)}{\log(1-p)} \leq i + 1 \]

ceea ce conduce la \(X = \left[\frac{\log(1-U)}{\log(1-p)}\right]\). Mai mult, pentru \(U\sim\mathcal{U}([0,1])\) avem \(1-U\sim\mathcal{U}([0,1])\) prin urmare \(X = \left[\frac{\log(U)}{\log(1-p)}\right]\). Următoarea funcție implementează acest algoritm:

sim_discrete_geom <- function(n = 100, p = 0.5){
     # n - nr de observatii
     # p - probabilitatea de succes
     
     u <- runif(n)
     x <- floor(log(u)/log(1-p))
     
     return(x)
}

În R putem folosi funcția rgeom() pentru a genera observații repartizate geometric.

Repartiția Hipergeometrică

Funcția generică GenerateDiscrete aplicată pentru \(\mathrm{HG}(20,30,15)\) conduce la:

Figura 9: Compararea repartiției teoretice \(\mathrm{HG}(20,30,15)\) cu cea empirică pentru un eșantion de volum \(n=10000\).

Ca și în cazul repartiției binomiale și după cum vom vedea mai jos a repartiției Poisson, și în cazul repartiției hipergeometrice putem utiliza relația de recurență pentru funcția de masă a \(\mathrm{HG}(n, N, M)\) în punctele \(i+1\) și \(i\). Pentru aceasta trebuie să ținem cont că în R repartiția hipergeometrică primește ca argumente de intrare \(N\) numărul de bile albe (nu numărul total de bile cum am văzut la curs!), \(M\) numărul de bile negre și \(n\) numărul de bile extrase cu întoarcere, i.e. \(\mathrm{HG}(n, N+M, M)\). Astfel funcția de masă devine

\[ \mathbb{P}(X = i) = f(i) = \frac{\binom{M}{i}\binom{N}{n-i}}{\binom{N+M}{n}} \]

iar între \(f(i)\) și \(f(i+1)\) există relația

\[ f(i+1) = \frac{M-i}{i+1}\frac{n-i}{N-n+i+1}f(i), \quad i\geq 0 \]

cu \(f(0)=\frac{\binom{N}{n}}{\binom{N+M}{n}}\). Folosind Teorema de universalitate deducem că pentru \(u\) uniform în \([0,1]\), \(X\) ia valoarea \(i\) care verifică

\[ F(i-1) < u \leq F(i), \]

de unde găsim următorul algoritm:

sim_discrete_hyper <- function(nn = 100, n = 1, N = 10, M = 5){
     # nn - nr de observatii generate
     # n - nr de bile extrase fara intoarcere
     # N - nr de bile albe in urna
     # M - nr de bile negre in urna
     
     x <- rep(0, nn)
     
     for (j in 1:nn){
          u <- runif(1)
          
          x[j] <- 0
          pb <- dhyper(0, M, N, n)
          Fr <- pb
          
          while(u >= Fr){
               pb <- (M - x[j])/(x[j] + 1) * (n - x[j])/(N - n + x[j] + 1) * pb
               Fr <- Fr + pb
               x[j] <- x[j] + 1
          }
     }
     
     return(x)
}

În R putem folosi și funcția rhyper() pentru a genera observații repartizate hipergeometric.

Repartiția Poisson

Aplicând funcția GenerateDiscrete pentru \(\mathrm{Pois}(5)\) găsim că:

Figura 10: Compararea repartiției teoretice \(\mathrm{Pois}(5)\) cu cea empirică pentru un eșantion de volum \(n=10000\).

Reamintim că o variabilă repartizată \(\mathrm{Pois}(\lambda)\) are funcția de masă dată de

\[ \mathbb{P}(X = i) = f(i) = e^{-\lambda}\frac{\lambda^i}{i!},\quad i\geq 0. \]

O metodă mai eficientă de simulare combină rezultatul din Teorema de universalitate și relația de recurență pentru funcția de masă următoare

\[ f(i+1) = \frac{\lambda}{i+1}f(i),\quad i\geq 0, \quad f(0) = e^{-\lambda}. \]

Astfel, pentru \(u\) o observație uniformă pe \([0,1]\), o observație \(i\) din repartiția \(\mathrm{Pois}(\lambda)\) verifică

\[ \mathbb{P}(X\leq i-1) = F(i-1) < u \leq F(i) = \mathbb{P}(X\leq i) \]

ceea ce conduce la algoritmul

sim_discrete_pois1 <- function(n = 100, lambda = 1){
     # n - nr de observatii
     # lambda - media Poisson
     
     x <- rep(0, n)
     
     for (j in 1:n){
          u <- runif(1)
          x[j] <- 0
          p <- exp(-lambda)
          Fr <- p
          
          while(u >= Fr){
               p <- p*lambda/(x[j] + 1)
               Fr <- Fr + p
               x[j] <- x[j] + 1
          }
     }
     
     return(x)
}

O altă modalitate de generare a observațiilor din repartiția \(\mathrm{Pois}(\lambda)\) este dată în ?@sec-gen-inv-other-pois. Totodată, în R putem folosi și funcția rpois() pentru a genera observații repartizate Poisson.

Generarea variabilelor aleatoare continue

Exercițiul 5 Scrieți un program care să folosească metoda transformării inverse pentru a genera \(n\) observații din densitatea

\[ f(x) = \left\{\begin{array}{ll} \frac{1}{x^2}, & x\geq 1\\ 0, & \text{altfel} \end{array}\right. \]

Testați programul trasând o histogramă a \(10000\) de observații aleatoare împreună cu densitatea teoretică \(f\).

Soluție. Primul pas este să determinăm funcția de repartiție \(F\) corespunzătoare acestei densități. Pentru \(x<1\) avem că \(f(x)=0\) deci \(F(x)=0\) iar pentru \(x\geq 1\) avem

\[ F(x) = \int_{1}^{x}\frac{1}{t^2}\, dt = 1 - \frac{1}{x}. \]

Cum \(F\) este continuă putem să determinăm \(F^{-1}\) rezolvând ecuația \(F(x)=u\). Un calcul direct conduce la \(F^{-1}(u)=\frac{1}{1-u}\) și aplicând Teorema de universalitate găsim că \(X = \frac{1}{1-U}\) cu \(U\sim \mathcal{U}([0,1])\).

Astfel putem simula un eșantion de volum \(n\) din populația \(f\) construind funcția

sim_f <- function(n){
  u <- runif(n)
  return(1/(1-u))
}

Pentru a testa comparăm valorile simulate cu densitatea teoretică

Figura 11: Compararea histogramei valorilor simultate cu repartiția teoretică.

Repartiția Exponențială

Exercițiul 6 Scrieți o funcție care să folosească metoda transformării inverse pentru a genera \(n\) observații dintr-o repartiție \(\mathrm{Exp}(\lambda)\).

Soluție. Reamintim că o variabilă aleatoare \(X\sim\mathrm{Exp}(\lambda)\) dacă admite o densitate de repartiție de forma

\[ f(x) = \lambda e^{-\lambda x}\mathbb{1}_{\mathbb{R}_+}(x),\quad \forall x\in\mathbb{R}. \]

Pentru a aplica Teorema de uniersalitate va trebui să determinăm funcția cuantilă. Știm că funcția de repartiție a repartiției exponențiale este

\[ F(x) = (1 - e^{-\lambda x})\mathbb{1}_{\mathbb{R}_+}(x), \quad x\in \mathbb{R} \]

ceea ce conduce la

\[ F^{-1}(u) = -\frac{1}{\lambda}\log(1-u), \quad u\in [0,1]. \]

Astfel putem genera o variabilă aleatoare \(X\) repartizată \(\mathrm{Exp}(\lambda)\), generând \(U\sim\mathcal{U}([0,1])\) și, ținând seama că și \(1-U\sim\mathcal{U}([0,1])\), luând

\[ X = -\frac{1}{\lambda}\log(U). \]

Următorul cod implementează funcția din enunțul exercițiului:

sim_cont_exp <- function(n, lambda){
     # n - nr de observatii
     # lambda - parametrul repartitiei Exponentiale
     
     u <- runif(n)
     x <- -log(u)/lambda
     
     return(x)
}

Pentru a testa comparăm valorile simulate cu densitatea teoretică

Figura 12: Compararea histogramei valorilor simultate cu repartiția teoretică \(\mathrm{Exp(5)}\).

Să observăm că în R putem folosi funcția rexp(n, rate = ...) pentru generarea de observații repartizate exponențial.

Repartiția Gamma

Exercițiul 7 Fie \(X_1, \ldots, X_n\), \(n\) variabile aleatoare i.i.d. repartizate \(\mathcal{E}(\lambda)\). Atunci variabila aleatoare \(S_n = X_1 + X_2 +\cdots+ X_n\) este repartizată \(\Gamma(n, \lambda)\). Plecând de la acest rezultat, generați \(10000\) de observații din repartiția \(\Gamma(n, \lambda)\) (\(n = 10\)). Trasați histograma acestui eșantion și comparați cu densitatea repartiției.

Soluție. Faptul că \(X_1 + X_2 +\cdots+ X_n\sim\Gamma(n, \lambda)\) este arătat în ?@exr-stat-sim-met-part-pois-e1. Pentru generarea a \(m\) variabile aleatoare i.i.d. repartizate \(\Gamma(n, \lambda)\), simulăm \(m\times n\) variabile repartizate exponențial pe care le stocăm într-o matrice cu \(m\) linii și \(n\) coloane. Suma elementelor de pe o linie reprezintă o realizare a repartiției \(\Gamma(n, \lambda)\).

sim_cont_gamma <- function(m, n, lambda = 1){
  out <- rowSums(matrix(sim_cont_exp(m * n, lambda), m, n))
  return(out)
}

Din punct de vedere al costului de memorie, această metodă poate să nu fie optimă. Putem folosi o abordare mai simplă folosind bucla for:

sim_cont_gamma2 <- function(m, n, lambda = 1){
  out <- numeric(m)
  
  for (i in 1:m){
    out[i] <- sum(sim_cont_exp(n, lambda))
  }
  
  return(out)
}

Pentru verificare considerăm repartiția \(\Gamma(10,1)\) și avem

Figura 13: Compararea histogramei valorilor simultate cu repartiția teoretică \(\Gamma(10,1)\).

În R putem folosi funcția rgamma(n, ...) pentru generarea de observații repartizate Gamma.

Repartiția Weibull

Exercițiul 8 Scrieți o funcție care să folosească metoda transformării inverse pentru a genera \(n\) observații dintr-o repartiție Weibull \(\mathrm{Weib}(\alpha, \beta)\).

Soluție. Spunem că o variabilă aleatoare \(X\) este repartizată Weibull de parametrii \(\alpha>0\) și \(\beta>0\), notat \(X\sim\mathrm{Weib}(\alpha, \beta)\), dacă admite ca densitate de repartiție pe

\[ f(x) = \frac{\alpha}{\beta}\left(\frac{x}{\beta}\right)^{\alpha-1}e^{-\left(\frac{x}{\beta}\right)^{\alpha}}\mathbf{1}_{[0,\infty)}(x). \]

Funcția de repartiție este dată de

\[ F(x) = \int_0^x \frac{\alpha}{\beta}\left(\frac{x}{\beta}\right)^{\alpha-1} \exp \left(-\left(\frac{t}{\beta}\right)^\alpha\right) \mathrm{d} t = 1-e^{ -\left(\frac{x}{\beta}\right)^\alpha}\mathbf{1}_{[0,\infty)}(x) \]

iar un calcul simplu conduce la funcția cuantilă

\[ F^{-1}(u) = \beta(-\log (1-u))^{1 / \alpha}. \]

Astfel pentru a genera \(X\sim\mathrm{Weib}(\alpha, \beta)\), considerăm \(U\sim\mathcal{U}([0,1])\) și, ținând cont de \(1-U\sim\mathcal{U}([0,1])\), luăm

\[ X = \beta(-\log (U))^{1 / \alpha}. \]

Următoarea funcție implementează această procedură:

sim_cont_weib <- function(n = 1, alpha = 1, beta = 1){
     # n - nr de obs 
     # alpha - parametrul de forma
     # beta - parametrul de scala
     
     u <- runif(n)
     x <- beta*(-log(u))^(1/alpha)
     
     return(x)
}

Pentru a testa comparăm valorile simulate cu densitatea teoretică

Figura 14: Compararea histogramei valorilor simultate cu repartiția teoretică \(\mathrm{Weib}(2, 3)\).

Să observăm că pentru \(\alpha = 1\), densitatea de repartiție a repartiției Weibull coincide cu cea a repartiției Exponențiale. De asemenea, în R putem folosi funcția rweibull() pentru a genera observații din \(\mathrm{Weib}(\alpha, \beta)\).

Repartiția Cauchy

Exercițiul 9 Scrieți o funcție care să folosească metoda transformării inverse pentru a genera \(n\) observații dintr-o repartiție Cauchy de parametrii \(\alpha\) și \(\beta\) - \(\mathrm{C}(\alpha, \beta)\).

Soluție. Reamintim că variabila aleatoare \(X\) este repartizată Cauchy de parametrii \((\alpha, \beta)\), cu \(-\infty<\alpha<\infty\) și \(\beta>0\), și notăm \(X\sim \mathrm{C}(\alpha, \beta)\), dacă densitatea ei de repartiție este dată de

\[ f_X(x) = \frac{1}{\pi\beta} \frac{1}{1+\left(\frac{x-\alpha}{\beta}\right)^2},\quad \forall x\in\mathbb{R}. \]

Trebuie specificat că media și varianța variabilei aleatoare \(\mathrm{C}(\alpha, \beta)\) nu există. Funcția de repartiție este dată de

\[ F(x) = \int_{-\infty}^x \frac{1}{\beta \pi}\left(1+\frac{(t-\alpha)^2}{\beta^2}\right)^{-1} \mathrm{~d} t = \frac{1}{2}+\frac{1}{\pi} \arctan\left(\frac{x-\alpha}{\beta}\right) \]

iar funcția cuantilă devine

\[ F^{-1}(u) = \alpha+\beta \tan \left(\pi\left(u-\frac{1}{2}\right)\right). \]

Astfel din Teorema de uniersalitate, generând \(U\sim\mathcal{U}([0,1])\) avem că \(X\) definit de relația de mai sus este repartizat \(\mathrm{C}(\alpha, \beta)\). Următoarea funcție implementează acest procedeu:

sim_cont_cauchy <- function(n = 1, alpha = 0, beta = 1){
     # n - nr de obs 
     # alpha - parametrul de locatie
     # beta - parametrul de scala
     
     u <- runif(n)
     x <- alpha + beta*tan(pi*(u - 0.5))
     
     return(x)
}

Figura de mai jos compară histograma unui eșantion de volum \(10000\) de observații din \(\mathrm{C}(2, 3)\) cu repartiția teoretică:

Figura 15: Compararea histogramei valorilor simultate cu repartiția teoretică \(\mathrm{C}(2, 3)\).

În R putem folosi funcția rcauchy() pentru a simula observații din repartiția \(\mathrm{C}(\alpha, \beta)\). O altă modalitate de generare a observațiilor din repartiția Cauchy este dată în ?@exr-stat-sim-met-part-cauchy-e1.

Repartiția Logistică

Exercițiul 10 Scrieți o funcție care să folosească metoda transformării inverse pentru a genera \(n\) observații dintr-o repartiție Logistică de parametrii \(\alpha\) și \(\beta\), \(\mathrm{Logis}(\alpha, \beta)\).

Soluție. Spunem că o variabilă aleatoare \(X\) este repartizată Logistic de parametrii \(-\infty<\alpha<\infty\) și \(\beta>0\), notat \(X\sim\mathrm{Logis}(\alpha, \beta)\), dacă admite ca densitate de repartiție pe

\[ f(x) = \frac{\exp \left(\frac{x-\alpha}{\beta}\right)}{\beta\left(1+\exp \left(\frac{x-\alpha}{\beta}\right)\right)^2}\mathbf{1}_{\mathbb{R}}(x). \]

Funcția de repartiție este dată de

\[ F(x) = \int_{-\infty}^x \frac{\exp \left(\frac{x-\alpha}{\beta}\right)}{\beta\left(1+\exp \left(\frac{x-\alpha}{\beta}\right)\right)^2} \mathrm{~d} t = \frac{1}{1+e^{-(x-\alpha) / \beta}} \]

iar un calcul direct conduce la funcția cuantilă

\[ F^{-1}(u) = \alpha+\beta \log \left(\frac{u}{1-u}\right). \]

Astfel pentru a genera \(X\sim\mathrm{Logis}(\alpha, \beta)\), considerăm \(U\sim\mathcal{U}([0,1])\) și luăm

\[ X = \alpha+\beta \log \left(\frac{U}{1-U}\right). \]

Următoarea funcție implementează această procedură:

sim_cont_logis <- function(n = 1, alpha = 0, beta = 1){
     # n - nr de obs 
     # alpha - parametrul de locatie
     # beta - parametrul de scala
     
     u <- runif(n)
     x <- alpha + beta*log(u/(1-u))
     
     return(x)
}

Pentru a testa funcția, comparăm valorile simulate pentru un eșantion de volum \(n=10000\) de observații cu densitatea teoretică

Figura 16: Compararea histogramei valorilor simultate cu repartiția teoretică \(\mathrm{Logis}(2, 3)\).

În R putem folosi funcția rlogis() pentru a genera observații din \(\mathrm{Logis}(\alpha, \beta)\).

Repartiția Pareto

Exercițiul 11 Scrieți o funcție care să folosească metoda transformării inverse pentru a genera \(n\) observații dintr-o repartiție Pareto de parametrii \(\alpha\) și \(\beta\), \(\mathrm{Pareto}(\alpha, \beta)\).

Soluție. Spunem că o variabilă aleatoare \(X\) este repartizată Pareto de parametrii \(\alpha>0\) și \(\beta>0\), notat \(X\sim\mathrm{Laplace}(\alpha, \beta)\), dacă admite ca densitate de repartiție pe

\[ f(x) = \begin{cases}\beta \alpha^\beta x^{-\beta-1}, & \text { pentru } \alpha<x<\infty, \\ 0, & \text { altfel}\end{cases}. \]

Integrând densitatea, obținem că funcția de repartiție este dată de

\[ F(x)=\int_\alpha^x \beta \alpha^\beta t^{-\beta-1} \mathrm{~d} t=1-\left(\frac{\alpha}{x}\right)^\beta \]

ceea ce conduce la funcția cuantilă

\[ F^{-1}(u) = \alpha(1-u)^{-1 / \beta}. \]

Astfel pentru a genera \(X\sim\mathrm{Pareto}(\alpha, \beta)\), luăm \(U\sim\mathcal{U}([0,1])\) și, ținând cont că \(1-U\sim\mathcal{U}([0,1])\), rezultă

\[ X = \alpha(1-U)^{-1 / \beta}\sim\mathrm{Pareto}(\alpha, \beta). \]

Următoarea funcție implementează această procedură:

sim_cont_pareto <- function(n = 1, alpha = 0, beta = 1){
     # n - nr de obs 
     # alpha - parametrul de forma
     # beta - parametrul de scala
     
     u <- runif(n)
     x <- alpha*(u)^(-1/beta)
     
     return(x)
}

Pentru a testa funcția, comparăm valorile simulate pentru un eșantion de volum \(n=10000\) de observații cu densitatea teoretică

Figura 17: Compararea histogramei valorilor simultate cu repartiția teoretică \(\mathrm{Pareto}(2, 3)\).

În R nu avem funcții predefinite aferente repartiției Pareto (rpareto, dpareto, ppareto sau qpareto) dar, ca și în cazul repartiției Laplace, acestea pot fi construite pe baza formulelor de mai sus.

Repartiția Normală

Exercițiul 12 Scrieți o funcție care să folosească metoda transformării inverse pentru a genera \(n\) observații dintr-o repartiție Normală de parametrii \(\mu\) și \(\sigma^2\), \(\mathcal{N}(\mu, \sigma^2)\).

Vom prezenta mai jos, în ?@sec-gen-inv-other-norm, o metodă de generare a observațiilor dintr-o repartiție \(\mathcal{N}(0, 1)\) folosindu-ne de transformarea Box-Muller, unde plecând de la două observații uniforme independente se obțin două observații repartizate normal și independente. De asemenea, în Exercițiul 16, vom ilustra o altă metodă de generare a unei normale folosind metoda respingerii.

Soluție. În acest exemplu vom considera o metodă de generare bazată pe metoda inversă. Cu toate că în cazul repartiției normale nu avem o formă compactă a funcției de repartiție \(\Phi(x)\) sau a funcției cuantilă \(\Phi^{-1}(u)\) și prin urmare nu putem aplica în mod direct această metodă, în literatură au fost propuși mai mulți algoritmi de aproximare a lui \(\Phi^{-1}(u)\) (a se vedea de exemplu (Odeh, Evans, et al. 1974) sau (Wichura 1988)). Conform (Odeh, Evans, et al. 1974), dacă \(X\sim\mathcal{N}(0,1)\) atunci pentru \(10^{-20}<u<1-10^{-20}\) putem aproxima funcția cuantilă prin

\[ \Phi^{-1}(u) = y + \frac{S_4(y)}{T_4(y)} = y + \frac{a_0 + a_1 y + a_2 y^2 + a_3 y^3 + a_4 y^4}{b_0 + b_1 y + b_2 y^2 + b_3 y^3 + b_4 y^4} \]

unde \(y = \sqrt{-2\log(u)}\) iar \(S_4(y)\) și \(T_4(y)\) sunt polinoame de gradul \(4\). Conform autorilor, eroarea maximă de aproximare este de ordinul \(1.5\times 10^{-8}\). Următorul cod implementează această aproximare:

q_norm_OE <- function(p, mean = 0, sd = 1){
     # functia cuantila dupa Odeh 1974
     
     # coeficientii polinoamelor S_4 si T_4
     a0 <- -0.322232431088
     a1 <- -1.0
     a2 <- -0.342242088547
     a3 <- -0.204231210245e-1
     a4 <- -0.453642210148e-4
     b0 <- 0.993484626060e-1
     b1 <- 0.588581570495
     b2 <- 0.531103462366
     b3 <- 0.103537752850
     b4 <- 0.385607006340e-2
     
     ps <- p
     
     ps <- ifelse(ps > 0.5, 1-ps, ps)
     
     y <- sqrt(-2*log(ps))
     
     out <- y + ((((y*a4+a3)*y+a2)*y+a1)*y+a0) /
                        ((((y*b4+b3)*y+b2)*y+b1)*y+b0)
     
     out <- ifelse(ps == 0.5, 0, out)
     out <- ifelse(p < 0.5, - out, out)
     
     out <- mean + sd*out
     return(out)
}

Vom testa această funcție comparând cu rezultatul obținut prin aplicarea funcției qnorm() (care implementează algoritmul din (Wichura 1988)):

Tabelul 1: Compararea rezultatelor cu funcția qnorm

t	qnorm	qnorm_OE	error
0.05	-1.6448536	-1.6448536	0
0.10	-1.2815516	-1.2815516	0
0.15	-1.0364334	-1.0364334	0
0.20	-0.8416212	-0.8416212	0
0.25	-0.6744898	-0.6744897	0
0.30	-0.5244005	-0.5244005	0
0.35	-0.3853205	-0.3853205	0
0.40	-0.2533471	-0.2533471	0
0.45	-0.1256613	-0.1256613	0
0.50	0.0000000	0.0000000	0
0.55	0.1256613	0.1256613	0
0.60	0.2533471	0.2533471	0
0.65	0.3853205	0.3853205	0
0.70	0.5244005	0.5244005	0
0.75	0.6744898	0.6744897	0
0.80	0.8416212	0.8416212	0
0.85	1.0364334	1.0364334	0
0.90	1.2815516	1.2815516	0
0.95	1.6448536	1.6448536	0

Aplicând metoda inversă obținem următoarea funcție pentru generarea unei observații normale:

sim_cont_normal <- function(n, mean = 0, sd = 1){
     u <- runif(n)
     x <- q_norm_OE(u, mean = mean, sd = sd)
     
     return(x)
}

Verificăm algoritmul prin trasarea histogramei:

Figura 18: Ilustrarea metodei inverse pentru repartiția normală.

Generarea variabilelor aleatoare prin metoda respingerii

În această secțiune vom prezenta o serie de exerciții aplicative ale metodei de simulare bazate pe acceptare-respingere. Astfel, vom începe prin descrierea modului în care putem genera observații repartizate uniform pe o mulțime din \(\mathbb{R}^d\) în Secțiunea 1.3.1, continuând apoi cu descrierea principiului metodei de acceptare-respingere pentru cazul variabilelor aleatoare care admit densitate în Secțiunea 1.3.2 și încheind cu aplicarea metodei pentru cazul discret în Secțiunea 1.3.3.

Generarea variabilelor aleatoare repartizate uniform pe o mulțime dată

Vom începe această secțiune prin a ne reaminti definiția repartiției uniforme pe o mulțime boreliană dată:

Definiția 1 Fie \((\Omega, \mathcal{F}, \mathbb{P})\) un câmp de probabilitate și \(B\in\mathcal{B}_{\mathbb{R}^d}\) o mulțime boreliană astfel încât \(0< \lambda_d(B)<\infty\) (\(\lambda_d\) este măsura Lebesgue pe \(\mathbb{R}^d\)). Spunem că vectorul aleator \(M:\Omega\to\mathbb{R}^d\) este repartizat uniform pe \(B\), și notăm \(M\sim \mathcal{U}(B)\), dacă repartiția sa verifică

\[ \left(\mathbb{P}\circ M^{-1}\right)(A) = \mathbb{P}(M\in A) = \frac{\lambda_d(A\cap B)}{\lambda_d(B)}, \quad \forall\, A\in\mathcal{B}_{\mathbb{R}^d}. \]

Remarcă

Dacă vectorul \(M\sim \mathcal{U}(B)\) atunci el admite densitatea de repartiție

\[ f_{M}(x) = \frac{1}{\lambda_d(B)}\mathbf{1}_{B}(x). \]

Următorul rezultat stă la baza metodei de generare a observațiilor repartizate uniform pe o mulțime boreliană dată:

Propoziția 1 Fie \((\Omega, \mathcal{F}, \mathbb{P})\) un câmp de probabilitate și \((M_n)_n\) un șir de vectori aleatori, \(M:\Omega\to\mathbb{R}^d\), independenți și de repartiție comună \(\mu\). Fie \(B\in\mathcal{B}_{\mathbb{R}^d}\) o mulțime boreliană astfel încât \(\mu(B)>0\). Pentru toți \(\omega\in\Omega\) luăm

\[ T(\omega) = \inf\{n\in\mathbb{N}^{*}\,|\,M_n(\omega)\in B\} \] adoptând convenția \(\inf\emptyset = \infty\). Definim \(M_{T}:\Omega\to\mathbb{R}^d\) prin

\[ M_{T}(\omega)=\left\{\begin{array}{ll} M_{T(\omega)}(\omega), & T(\omega)<\infty\\ 0, & T(\omega) = \infty \end{array}\right. \]

Atunci

1. \(\mathbb{P}(T<\infty)=1\) și dacă \(T'=T\mathbf{1}_{T<\infty}\) atunci \(T'\sim\mathrm{Geom}(\mu(B))\)

2. \(M_{T}\) este un vector aleator de repartiție \(\nu\) definită prin

\[ \nu(A) = \mathbb{P}(M_{T}\in A) = \frac{\mu(A\cap B)}{\mu(B)}, \quad \forall\, A\in\mathcal{B}_{\mathbb{R}^d}. \]

Cu alte cuvinte \(\nu\) este repartiția condiționată \(\mu(\cdot|B)\).

Remarcă

Dacă presupunem că \(C\in\mathbb{R}^d\) este un dreptunghi astfel încât \(B\subset C\) și \(0< \lambda_d(B)\leq \lambda_d(C)<\infty\) iar \(\mu\) este repartiția uniformă pe \(C\) atunci \(\nu\) este repartiția uniformă pe \(B\).

Ca aplicație a rezultatului de mai sus considerăm următorul exercițiu:

Exercițiul 13 Considerăm pătratul \(C = [0,L]^2\) și discul \(D\) de centru \((\frac{L}{2},\frac{L}{2})\) și rază \(\frac{L}{2}\). Considerăm șirul de v.a. \(\left(Y_n\right)_{n\geq1}\) pe \(\mathbb{R}^2\) i.i.d. repartizate uniform pe pătratul \(C\).

1. Aproximați valoarea lui \(\pi\) prin ajutorul numărului de puncte \(Y_n\) care cad în interiorul discului \(D\) (Metoda respingerii)
2. Simulați \(n\) puncte uniforme pe disc.

Soluție. Pentru primul punct avem:

1. Definim v.a. \(X_n=\mathbf{1}_{\{Y_n\in D\}}\), \(n\geq1\), care formează un șir de v.a. i.i.d. de repartiție \(\mathcal{B}(\mathbb{P}(Y_n\in D))\), deoarece \(\left(Y_n\right)_{n\geq1}\) este un șir de v.a. i.i.d. repartizate uniform pe \(C\), \(\mathcal{U}(C)\). Din Legea Numerelor Mari avem că

\[ \displaystyle\frac{1}{n}\sum_{i=1}^{n}X_{i} \overset{a.s.}{\to} \mathbb{E}[X_1] = \mathbb{P}(Y_1\in D), \]

prin urmare trebuie să calculăm probabilitatea \(\mathbb{P}(Y_1\in D)\). Știm că densitatea v.a. \(Y_1\) este dată de \(f_{Y_1}(x,y)=\frac{1}{\mathcal{A}(C)}\mathbf{1}_{C}(x,y)\) de unde

\[ \begin{aligned} \mathbb{P}(Y_1\in D) &= \iint_{D}f_{Y_1}(x,y)\,dxdy = \iint \mathbf{1}_{D}(x,y)\mathbf{1}_{C}(x,y)\,dxdy\\ &= \frac{1}{\mathcal{A}(C)}\iint \mathbf{1}_{D}(x,y)\,dxdy = \frac{\mathcal{A}(D)}{\mathcal{A}(C)} = \frac{\pi \frac{L^2}{4}}{L^2} = \frac{\pi}{4}. \end{aligned} \]

Astfel, putem estima valoarea lui \(\pi\) prin \(\displaystyle\frac{4}{n}\sum_{i=1}^{n}X_{i}\) pentru valori mari ale lui \(n\).

# Estimam valoarea lui pi
L <- 3 # lungimea laturii patratului 
R <- L/2 # raza cercului inscris

n <- 2000 # numarul de puncte din patratul C

# generam puncte uniforme in C
x <- L*runif(n)
y <- L*runif(n)

# metoda respingerii 
l <- (x-R)^2+(y-R)^2 # distanta dintre centrul cercului si punct
ind <- l<=(R)^2 # indicii pentru care distanta este mai mica sau egala cu R

xc <- x[ind] # coordonatele punctelor din interiorul cercului  
yc <- y[ind] 

estimate_pi <- 4*sum(ind)/n # estimarea lui pi
err <- abs(estimate_pi-pi) # eroarea absoluta

Aplicând acest procedeu obținem că valoarea estimată a lui \(\pi\) prin generarea a \(n=\) 2000 puncte este 3.096 iar eroarea absoluta este 0.04559.

Pentru punctul doi:

2. Una dintre metodele prin care putem simula puncte repartizate uniform pe suprafața discului \(D\) este Metoda respingerii așa cum este descrisă la începutul secțiunii. Această metodă consistă în generarea de v.a. \(Y_n\) repartizate uniform pe suprafața pătratului \(C\), urmând ca apoi să testăm dacă \(Y_n\) aparține discului \(D\) (deoarece \(D\subset C\)). Dacă da, atunci le păstrăm, dacă nu, atunci mai generăm. Următoarea figură ilustrează această metodă:

Figura 19: Ilustrarea metodei de simulare a punctelor repartizate uniform pe disc prin metoda respingerii.

Funcția de mai jos permite generarea a \(n\) observații repartizate uniform pe discul \(D\) de centru \((xc, yc)\) și rază \(R\):

gen_unif_disc <- function(n = 10, R = 2, xc = 2, yc = 2){
     # n - nr observatii
     # R - raza cercului 
     # xc, yc - centrul cercului 
     # j (output) - nr de obs necesar pentru generarea unei obs
     
  out <- matrix(0, ncol = 3, nrow = n)
  
  for (i in 1:n){
    
    d <- R^2 + 1
    j <- 0
    
    while(d > R^2){
      x <- xc - R + 2*R*runif(1)
      y <- yc - R + 2*R*runif(1)
    
      d <- (x - xc)^2 + (y - yc)^2
      
      j <- j + 1
    }
    
    out[i, ] <- c(x, y, j)
  }
  
  return(out)
}

Putem observa că numărul mediu de puncte repartizate uniform pe \(C\) necesare pentru a genera un punct uniform pe \(D\) este \(4/\pi\) (raportul dintre aria pătratului și aria discului). Acest fenomen se poate observa și prin simulare. Atunci când generăm \(n=1000\) de puncte repartizate uniform pe discul \(D\) de centru \((xc,yc)=(\frac{L}{2},\frac{L}{2})=(1.5,1.5)\) și rază \(R=\frac{L}{2} = 1.5\), obținem un număr mediu de puncte de 1.279, cu o eroare de 0.0057605.

Metodă alternativă: Vom prezenta mai jos o altă metodă de simulare a punctelor distribuite uniform pe discul \(D\) de rază \(L\).

O primă idee, eronată, ar fi să generăm cuplul de v.a. \((X_1,Y_1)\) astfel încât, ținând seama de schimbarea de variabile în coordonate polare, \(X_1 = \tilde{R}\cos(\tilde{\Theta})\) și \(Y_1 = \tilde{R}\sin(\tilde{\Theta})\) unde \(\tilde{R}\sim\mathcal{U}([0,L])\), \(\tilde{\Theta}\sim\mathcal{U}([0,2\pi])\) cu \(tilde{R}\) și \(\tilde{\Theta}\) independente (ceea ce nu este adevărat în realitate). Vom vedea (printr-o ilustrație grafică) că această abordare este greșită (punctele sunt concentrate în centrul cercului). O funcție care să implementeze această abordare ar putea fi:

gen_unif_disc_wrong <- function(n = 10, R = 2, xc = 2, yc = 2){
  
  r <- R*runif(n)
  theta <- 2*pi*runif(n)
  
  x <- xc + r*cos(theta)
  y <- yc + r*sin(theta)
  
  return(cbind(x, y))
}

O abordare corectă este următoarea: căutăm să simulăm un cuplu de v.a. \((X,Y)\) care este uniform distribuit pe suprafața discului \(D\), i.e. densitatea cuplului este dată de \(f_{(X,Y)}(x,y)=\frac{1}{\pi L^2}\mathbf{1}_{D}(x,y)\). Considerând schimbarea de variabile în coordonate polare: \(x=r\cos(\theta)\) și \(y=r\sin(\theta)\), obiectivul este de a găsi densitatea variabilelor \(R\) și \(\Theta\) astfel încât \(X = R\cos(\Theta)\) și \(Y = R\sin(\Theta)\).

Fie \(g(x,y)=\left(\sqrt{x^2+y^2},\arctan(y/x)\right)=(r,\theta)\), transformarea pentru care avem \((R,\Theta)=g(X,Y)\). Știm că inversa acestei transformări este \(g^{-1}(r,\theta)=(r\cos(\theta),r\sin(\theta))\), prin urmare

\[ \begin{aligned} f_{(R,\Theta)}(r,\theta) &= f_{(X,Y)}\left(g^{-1}(r,\theta)\right)|\det(J_{g^{-1}}(r,\theta))|\\ &= \frac{1}{\pi L^2}\mathbf{1}_{D}(r\cos(\theta),r\sin(\theta))\left|\begin{array}{cc} \cos(\theta) & \sin(\theta)\\ r\sin(\theta) & -r\cos(\theta) \end{array}\right|\\ &= \frac{1}{\pi L^2} \mathbf{1}_{[0,L]}(r)\mathbf{1}_{[0,2\pi]}(\theta)r. \end{aligned} \]

Observăm că densitatea (marginală) a v.a. \(\Theta\) este

\[ \begin{aligned} f_{\Theta}(\theta) &= \int f_{(R,\Theta)}(r,\theta)\,dr = \mathbf{1}_{[0,2\pi]}(\theta)\int \frac{r}{\pi L^2} \mathbf{1}_{[0,L]}(r)\,d\theta\\ &= \frac{1}{\pi L^2} \mathbf{1}_{[0,2\pi]}(\theta) \frac{L^2}{2} = \frac{1}{2\pi} \mathbf{1}_{[0,2\pi]}(\theta), \end{aligned} \]

iar densitatea v.a. \(R\) este

\[ \begin{aligned} f_{R}(r) &= \int f_{(R,\Theta)}(r,\theta)\,d\theta = \frac{r}{\pi L^2} \mathbf{1}_{[0,L]}(r)\int_{0}^{2\pi}\,d\theta\\ &= \frac{r}{\pi L^2} \mathbf{1}_{[0,L]}(r)2\pi = \frac{2r}{L^2} \mathbf{1}_{[0,L]}(r). \end{aligned} \]

Din expresiile de mai sus putem observa că \(\Theta\) este o v.a. repartizată uniform pe \([0,2\pi]\) și putem verifica ușor că repartiția v.a. \(R\) este aceeași cu cea a v.a. \(L\sqrt{U}\) unde \(U\sim\mathcal{U}([0,1])\).

Astfel, pentru simularea unui punct \((X,Y)\) uniform pe \(D\) este suficient să simulăm o v.a. \(\Theta\) uniform pe \([0,2\pi]\) și o v.a. \(U\) uniformă pe \([0,1]\) și să luăm \(X=L\sqrt{U}\cos(\Theta)\) și \(Y=L\sqrt{U}\sin(\Theta)\). Funcția de mai jos implementează această metodă:

gen_unif_disc_good <- function(n = 10, R = 2, xc = 2, yc = 2){
  
  r <- R*sqrt(runif(n))
  theta <- 2*pi*runif(n)
  
  x <- xc + r*cos(theta)
  y <- yc + r*sin(theta)
  
  return(cbind(x, y))
}

Următoarea figură ne ilustrează cele două proceduri prezentate:

Figura 20: Ilustrarea metodei de simulare a punctelor repartizate uniform pe discul de centru (5, 9) si raza 12.

Generarea variabilelor aleatoare ce admit densitate

Metoda respingerii permite și simularea variabilelor aleatoare sau a vectorilor aleatori din repartițiile care admit densitate în raport cu măsura Lebesgue. Considerăm că vrem să generăm observații din vectorul aleator \(Z\) cu valori în \(\mathbb{R}^d\) care admite densitatea \(f\) (densitate țintă). Vom presupune că știm să simulăm un vector aleator \(X\) care admite densitatea \(g\) (densitate propusă sau propunere) și că există o constantă \(c\geq1\) care verifică relația \(f\leq cg\). Metoda respingerii propune să generăm un șir de vectori aleatori \((X_n)_{n\geq 1}\) independenți de repartiție comună care admite ca densitate pe \(g\) și \((U_n)_{n\geq 1}\) un șir de variabile aleatoare independente (și independente și de \((X_n)_{n\geq 1}\)) și repartizate \(\mathcal{U}([0,1])\). Generăm punctele \(M_i=(X_i, cg(X_i)U_i)\) și ne oprim la primul indice \(i_0\) (aleator) care verifică \(cg(X_{i_0})U_{i_0}\leq f(X_{i_0})\). Luăm \(Z = X_{i_0}\), iar acest vector aleator admite ca densitate pe \(f\) (conform Propoziția 2).

Justificarea algoritmului are la bază următoarele propoziții:

Propoziția 2 Fie \(f\) o densitate de probabilitate pe \(\mathbb{R}^d\) și considerăm hipograful lui \(f\)

\[ G=\left\{(x, y)\in\mathbb{R}^d\times\mathbb{R}\,|\,0\leq y\leq f(x)\right\}. \]

Fie \(M = (Z,Y)\) un vector aleator pe \(\mathbb{R}^d\times\mathbb{R}\) repartizat uniform pe \(G\). Atunci repartiția lui \(Z\) pe \(\mathbb{R}^d\) admite densitatea \(f\) în raport cu măsura Lebesgue \(\lambda_d\) pe \(\mathbb{R}^d\).

Propoziția 3 Fie \(X\) un vector aleator pe \(\mathbb{R}^d\) cu densitatea de repartiție \(g\) și \(U\sim\mathcal{U}([0,1])\) independent de \(X\). Considerăm punctul \(M = (X, cg(X)U)\), cu \(c>0\), și fie \(H\) hipograful lui \(cg\), i.e.

\[ H=\left\{(x, y)\in\mathbb{R}^d\times\mathbb{R}\,|\,0\leq y\leq cg(x)\right\}. \]

Atunci punctul \(M\sim\mathcal{U}(H)\).

Ca aplicație a rezultatelor de mai sus considerăm o serie de exerciții:

Exercițiul 14 (Repartiția triunghiulară) Folosind metoda respingerii, propuneți o metodă de simulare pentru observații independente din densitatea de repartiție \(f:x\mapsto (1-|x|)^+\).

Soluție. Observăm că pentru toți \(x\in\mathbb{R}\) avem

\[ f(x) \leq \mathbf{1}_{\{x\in [-1, 1]\}} = 2 g(x), \]

unde \(g(x) = \frac{1}{2}\mathbf{1}_{\{x\in [-1, 1]\}}\) este densitatea repartiției uniforme pe \([-1,1]\).

Figura 21: Ilustrarea metodei respingerii pentru repartiția \(f(x)=(1-|x|)^+\).

Pentru a simula din repartiția \(f\) procedăm astfel

1. simulăm \(X\) repartizată \(\mathcal{U}[-1,1]\) (\(X = 2V-1\) cu \(V\sim\mathcal{U}(0,1)\))

2. simulăm \(U\sim\mathcal{U}(0,1)\)

3. repetăm procedeul până când \(2Ug(X)\leq f(X)\).

Următorul cod implementează acest algoritm:

f <- function(x){
  return(ifelse(abs(x)<=1, 1-abs(x), 0))
}

g <- function(x){
  return(ifelse(abs(x)<=1, 0.5, 0))
}

sim_f <- function(n = 1000){
   x <- rep(0, n)
   
   for (i in 1:n){
      # gen obs din g 
      x[i] <- runif(1, -1, 1)
      
      # gen uniforma
      u <- runif(1)
      
      while(2*u*g(x[i]) > f(x[i])){
          x[i] <- runif(1, -1, 1)
          u <- runif(1)
      }
   }
   
   return(x)
}

Putem valida algoritmul propus trasând histograma unui eșantion de volum \(10000\) și suprapunând densitatea \(f\):

Figura 22: Histograma eșantionului generat pentru repartiția triunghilară \(f(x)=(1-|x|)^+\).

Exercițiul 15 (Repartiția semicirculară a lui Wigner³) Fie \(f\) densitatea de repartiție definită prin

\[ f(x) = \frac{2}{\pi}\sqrt{1-x^2}\mathbf{1}_{[-1,1]}(x) \]

1. Folosiți metoda respingerii pentru a genera \(10000\) de observații repartizate cu densitatea \(f\).
2. Trasați histograma acestui eșantion și comparați cu densitatea \(f\).

Soluție. Pentru primul punct avem:

1. Observăm că

\[ f(x) \leq \frac{2}{\pi}\mathbf{1}_{[-1,1]}(x) = \frac{4}{\pi}g(x) \]

unde \(g(x) = \frac{1}{2}\mathbf{1}_{[-1,1]}(x)\) este densitatea repartiției uniforme pe \([-1,1]\). Astfel metoda respingerii sugerează următorul algoritm:

f <- function(x){
  return(2/pi * sqrt(1-x^2)*(abs(x) <= 1))
}

sim_f <- function(n){
  x <- rep(0, n)
  
  for (i in 1:n){
    # generam obs din g
    x[i] <- runif(1, -1, 1)
    
    # generam uniforma
    u <- runif(1)
    
    while(u > pi * f(x[i])/2){
      x[i] <- runif(1, -1, 1)    
      u <- runif(1)
    }
  }
  
  return(x)
}

Putem îmbunătății codul de mai sus (vrem să evităm să avem și bucla for și bucla while imbricate) dacă ținem cont de faptul că probabilitatea de acceptare este \(p = \frac{\pi}{4}\) iar pentru a genera un eșantion de \(m\) observații avem nevoie, în medie, de \(\frac{m}{p}\) simulări. Avem

sim_f2 <- function(n){
  out <- c() # esantionul final
  
  # cate obs mai avem de generat
  m <- n
  
  # cat timp nu avem esantionul de talia dorita 
  # continuam procedeul 
  while (m>0){
    # simulam m/p observatii
    x <- runif((4*m)%/%pi + 1, -1, 1)
    u <- runif((4*m)%/%pi + 1)
    
    # testam care obs sunt acceptate
    y <- (u <= pi * f(x)/2)
    
    # pastram doar punctele acceptate
    out <- c(out, x[which(y)])
    m <- n - length(out)
  }
  
  return(out[1:n])
}

Putem compara cele două metode, în funcție de timpul de execuție (și observăm că în fapt a doua metodă este mai lentă):

n <- 10000

# Metoda 1
start <- proc.time()
x <- sim_f(n)
proc.time() - start

   user  system elapsed 
   0.05    0.00    0.04 

# Metoda 2
star <- proc.time()
x <- sim_f2(n)
proc.time() - start

   user  system elapsed 
   0.06    0.00    0.04 

2. Putem valida algoritmul propus prin metoda respingerii trasând histograma eșantionului:

Figura 23: Histograma eșantionului generat pentru repartiția \(f(x)\).

Exercițiul 16 (Repartiția normală - metoda respingerii 1) Plecând cu o propunere de tip \(\mathrm{Exp}(\lambda)\) vrem să generăm, cu ajutorul metodei acceptării-respingerii, un eșantion din următoarea densitate (jumătate de normală):

\[ f(x) = \left\{\begin{array}{ll} \frac{2}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}, & \mbox{dacă $x\geq0$}\\ 0, & \mbox{altfel}\\ \end{array}\right. \]

Cum putem modifica algoritmul pentru a genera un eșantion dintr-o repartiție normală standard?

Soluție. Pentru început, să observăm că densitatea țintă \(f\) corespunde lui \(|X|\), unde \(X\sim\mathcal{N}(0,1)\). Alternativ se poate interpreta drept densitatea unei normale standard \(X\sim\mathcal{N}(0,1)\) condiționată la \(X>0\).

Fie \(g\) densitatea repartiției exponențiale de parametru \(\lambda\),

\[ g(x) = \left\{\begin{array}{ll} \lambda e^{-\lambda x}, & \mbox{dacă $x\geq0$}\\ 0, & \mbox{altfel}\\ \end{array}\right. \]

Pentru a aplica algoritmul de acceptare-respingere trebuie să găsim valoarea lui \(c>0\) pentru care \(f(x)\leq c g(x)\) pentru toate valorile \(x\in \mathbb{R}\). Pentru \(x\geq0\) avem

\[ \frac{f(x)}{g(x)}=\frac{2}{\lambda\sqrt{2\pi}}e^{-\frac{x^2}{2}+\lambda x} \]

și cum funcția \(-\frac{x^2}{2}+\lambda x\) își atinge valoarea maximă în punctul \(x=\lambda\) rezultă că

\[ \frac{f(x)}{g(x)}\leq c^*, \,\,\forall x \geq0 \]

unde

\[ c^*=\sqrt{\frac{2}{\pi\lambda^2}}e^{\lambda^2/2}. \]

Astfel algoritmul devine:

• pentru \(n=1,2,\dots\)

• generează \(X_n\sim Exp(\lambda)\)

• generează \(U_n\sim\mathcal{U}[0,1]\)

• dacă \(U_n\leq\exp\left(-\frac{1}{2}(X_n-\lambda)^2\right)\) atunci

• întoarceți \(X_n\)

Avem funcția:

# generarea punctelor din densitatea f

f <- function(x) {
  return((x> 0) * 2 * dnorm(x,0,1))
}

g <- function(x) { return(dexp(x,1)) }

c_star <- sqrt(2 * exp(1) / pi)

rhalfnormal <- function(n) {
  res <- numeric(length=n)
  i <- 0
  while (i<n) {
    U <- runif(1, 0, 1)
    X <- rexp(1, 1)
    
    if (c_star * g(X) * U <= f(X)) {
      i <- i+1
      res[i] <- X
    }
  }
  return(res)
}

Pentru a testa validitatea algoritmului vom trasa histograma unui eșantion de volum \(n=10000\) peste care suprapunem densitatea de repartiție teoretică:

Figura 24: Histograma eșantionului generat pentru repartiția \(f(x)\) din enunț.

Cum densitatea \(f\) este densitatea unei normale standard \(X\sim\mathcal{N}(0,1)\) condiționată la \(X>0\) și cum densitatea normală este simetrică față de medie (0 în acest caz) algoritmul se modifică acceptând \(X_n\) și \(-X_n\) cu probabilitatea de \(0.5\).

Astfel avem funcția:

f2 <- function(x) {
  return(dnorm(x,0,1))
}

normal1 <- function(n) {
  res <- numeric(length=n)
  i <- 0
  while (i<n) {
    U <- runif(1, 0, 1)
    X <- rexp(1, 1)
    if (c_star * g(X) * U <= f(X)) {
      i <- i+1
      
      res[i] <- ifelse(runif(1) <= 0.5, X, -X);
    }
  }
  return(res)
}

Verificăm dacă algoritmul este valid prin generarea unui eșantion de volum \(n=10000\), trasarea histogramei corespunzătoare peste care suprapunem densitatea de repartiție teoretică:

Figura 25: Histograma eșantionului generat pentru repartiția normală.

Generarea variabilelor aleatoare discrete

Metoda respingerii nu se limitează doar la acele repartiții care admit densitate de repartiție, ea poate fi aplicată și în cazul repartițiilor discrete. Vom presupune, pentru simplificarea notațiilor, că repartițiile discrete pe care le considerăm iau valori într-o submulțime a numerelor naturale \(\mathbb{N}\) (în orice caz, o repartiție discretă ia valori într-o mulțime ce poate fi pusă în bijecție cu o submulțime a lui \(\mathbb{N}\)).

Astfel, ne propunem să simulăm observații dintr-o variabilă aleatoare \(Z\) a cărei repartiție este dată de \(Z(\Omega)\subset\mathbb{N}\) și de funcția de masă \(f:\mathbb{N}\to [0,1]\) cu \(f(k)=p_k=\mathbb{P}(Z=k)\). Pentru aceasta presupunem că știm să simulăm o variabilă aleatoare discretă \(X\), pentru care \(X(\Omega)\subset\mathbb{N}\) iar funcția de masă \(g:\mathbb{N}\to [0,1]\) este \(g(k)=q_k=\mathbb{P}(X=k)\) și verifică \(f\leq cg\) pentru un \(c\geq 1\) (aici considerăm că \(Z(\Omega)\subset X(\Omega)\)).

Propoziția 4 Folosind notațiile de mai sus, fie \((\Omega, \mathcal{F}, \mathbb{P})\) un câmp de probabilitate și \((X_n)_n\) un șir de variabile aleatoare independente de aceeași repartiție cu \(X\) și \((U_n)_n\) un șir de variabile aleatoare i.i.d. repartizate \(\mathcal{U}([0,1])\), cele două șiruri fiind și independente între ele. Definim variabila aleatoare \(T\) prin

\[ T(\omega) = \inf\{n\in\mathbb{N}^{*}\,|\,cg(X_n(\omega))U_n(\omega)\leq f(X_n(\omega))\} \]

adoptând convenția \(\inf\emptyset = \infty\) și \(X_{T}:\Omega\to\mathbb{R}^d\) prin

\[ X_{T}(\omega)=\left\{\begin{array}{ll} X_{T(\omega)}(\omega), & T(\omega)<\infty\\ 0, & T(\omega) = \infty \end{array}\right. \]

Atunci \(X_T\) este o variabilă aleatoare discretă cu aceeași repartiție ca \(Z\), i.e. \(\mathbb{P}(X_T=k)=f(k),\quad\forall k\in\mathbb{N}\).

Algoritmul de respingere se scrie succint astfel:

1. Generăm o observație din \(X\) cu funcția de masă \(q_k\).

2. Generăm o observație \(U\sim\mathcal{U}(0,1)\).

3. Dacă \(U < \frac{p_k}{cq_k}\) atunci \(Z = X\) altfel mergem la pasul 1 și repetăm.

Exercițiul 17 Folosind metoda respingerii, propuneți un algoritm de simulare a unui eșantion de \(n = 1000\) de observații independente din variabila aleatoare

\[ Z\sim\left(\begin{array}{cccccccccc} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10\\ 0.11 & 0.12 & 0.09 & 0.08 & 0.12 & 0.10 & 0.09 & 0.09 & 0.10 & 0.10 \end{array}\right) \]

Soluție. Am văzut deja în Exercițiul 4 o soluție pentru această problemă dar acum vom rezolva problema folosind metoda respingerii. Pentru aceasta vom considera că știm să generăm observații independente dintr-o variabilă aleatoare \(X\sim\mathcal{U}(\{1,2,\ldots,10\})\), adică \(q_k = \frac{1}{10}\) pentru \(k\in\{1,2,\ldots,10\}\). Pentru această alegere a lui \(X\), respectiv \(q_k\), putem alege constanta \(c\) astfel încât

\[ c = \max_{k}\frac{p_k}{q_k} = 10\max_{k}p_k = 10\times 0.12 = 1.2. \]

Algoritmul devine astfel:

1. Generăm \(X\sim\mathcal{U}(\{1,2,\ldots,10\})\) (de exemplu folosind relația din Secțiunea 1.2.1.1, \(X = [10U_1] + 1\) cu \(U_1\sim\mathcal{U}(0,1)\))
2. Generăm \(U\sim\mathcal{U}(0,1)\)
3. Dacă \(Ucg(X)\leq f(X)\) atunci \(Z=X\) altfel trecem la pasul 1

Relația din pasul 3 se rescrie sub forma

\[ U\times 1.2\times \frac{1}{10}\leq f(X) \iff U\leq \frac{f(X)}{0.12} \] Următoarea funcție implementează acest procedeu:

pX <- c(0.11, 0.12, 0.09, 0.08, 0.12, 
        0.1, 0.09, 0.09, 0.1, 0.1)

gen_MR_X <- function(n, pX = pX){
  lx <- length(pX)
  
  m <- 0
  x <- rep(0, n)
  
  # determinam constanta
  c1 <- lx*max(pX)
  m1 <- max(pX)
  
  i <- 0
  
  while(m <= n){
    
    y <- sample(1:lx, 1, replace = TRUE, prob = rep(1/lx, lx))
    u <- runif(1)
    i <- i+1 # nr obs gen

    if (u <= pX[y]/m1){
      m <- m + 1
      x[m] <- y
    }
  }
  
  
  return(list(x = x, 
              nr_iter = i, 
              nr_med = c1*n))
}

Să observăm că acest algoritm necesită în medie \(1.2\) iterații (în general \(c\) iterații) pentru genera o observație din \(Z\). Putem testa acest lucru pentru \(n = 100\) observații:

N1 <- rep(0, 100)

for (i in 1:100){
  x <- gen_MR_X(n = 100, pX = pX)
  N1[i] <- x$nr_iter
}

mean(N1)

[1] 120.57

Putem compara repartiția empirică cu cea teoretică:

Figura 26: Compararea repartiției teoretice cu cea empirică pentru un eșantion de volum \(n=10000\).

Exercițiul 18 (Repartiția Zipf) Folosind metoda respingerii generați \(n = 1000\) de observații din repartiția \(\mathrm{Zipf}(2)\).

Soluție. Repartiția Zipf este folosită cu precădere în lingvistică și științele sociale dar nu numai. Reamintim că o variabilă aleatoare \(Z\) repartizată \(\mathrm{Zipf}(a)\), are suportul pe mulțimea \(\mathbb{N}^*\) și funcția de masă dată de

\[ \mathbb{P}(Z = k) = \frac{1}{\zeta(a)k^a}, \quad k\in \mathbb{N}^* \]

unde \(\zeta(a) = \sum_{k\geq1}\frac{1}{k^a}\) este funcția \(\zeta\) a lui Riemann.

Vom considera cazul \(a=2\). În acest caz știm că \(\zeta(2) = \sum_{k\geq1}\frac{1}{k^2} = \frac{\pi^2}{6}\) prin urmare

\[ f(k) = p_k = \frac{6}{\pi^2 k^2}, \quad k\geq1. \]

Vom folosi ca propunere repartiția cu funcția de masă

\[ g(k) = q_k = \frac{1}{k(k+1)}, \quad k\geq1. \]

Alegem această propunere deoarece este ușor de simulat din ea, mai precis dacă luăm \(U\sim\mathcal{U}(0,1)\) avem

\[ \begin{aligned} q_k &= \frac{1}{k(k+1)} = \frac{1}{k} - \frac{1}{k+1} = \mathbb{P}\left(\frac{1}{k+1}< U\leq \frac{1}{k}\right)\\ &= \mathbb{P}\left(k\leq \frac{1}{U}< k+1\right) \end{aligned} \]

ceea ce arată că variabila aleatoare \(X=\left[\frac{1}{U}\right]\) are funcția de masă \(q_k\). În plus, putem determina cu ușurință constanta \(c\) care verifică

\[ c = \sup_{k}\frac{p_k}{q_k} = \frac{6}{\pi^2}\sup_{k}\frac{k+1}{k} = \frac{12}{\pi^2}. \]

Algoritmul devine:

1. Generăm \(X\) cu repartiția \(q_k\): luăm \(U\sim\mathcal{U}(0,1)\) și \(X=\left[\frac{1}{U}\right]\)
2. Generăm \(V\sim\mathcal{U}(0,1)\)
3. Dacă \(Vcg(X)\leq f(X)\) atunci \(Z=X\) altfel trecem la pasul 1

Pasul 3 se rescrie astfel

\[ Vcg(X)\leq f(X) \iff V\frac{12}{\pi^2}\frac{1}{X(X+1)}\leq \frac{6}{\pi^2 X^2} \]

ceea ce conduce la \(V\leq \frac{1}{2} + \frac{1}{2X}\).

Următoarea funcție implementează metoda respingerii pentru repartiția \(\mathrm{Zipf}(2)\):

gen_Zipf <- function(n){
  x <- rep(0, n)
  
  m <- 0
  
  # constanta
  c1 <- 12/pi^2
  
  iter <- 0
  
  while(m <= n){
    # nr iteratii
    iter <- iter + 1
    
    # gen Y din g = 1/k(k+1)
    
    y <- floor(1/runif(1))
    
    # gen V unif [0,1]
    
    v <- runif(1)
    
    # verificam conditia
    if (v <= 0.5+1/(2*y)){
      m <- m+1
      x[m] <- y
    }
    
  }
  
  return(list(x = x, 
              nr_iter = iter, 
              nr_med = n*c1))
}

Putem compara repartiția empirică cu cea teoretică:

Figura 27: Compararea repartiției teoretice Zipf(2) cu cea empirică pentru un eșantion de volum \(n=10000\).

Referințe

Odeh, Robert E, JO Evans, et al. 1974. “The Percentage Points of the Normal Distribution.” Journal of the Royal Statistical Society Series C 23 (1): 96–97.

Russell, K. G. 1991. “Estimating the Value of e by Simulation.” The American Statistician 45: 66–68. https://api.semanticscholar.org/CorpusID:123112105.

Wichura, Michael J. 1988. “Algorithm AS 241: The Percentage Points of the Normal Distribution.” Journal of the Royal Statistical Society. Series C (Applied Statistics) 37 (3): 477–84.

Note de subsol

1. Spunem că un estimator nedeplasat este mai eficient decât un altul dacă varianța lui este mai mică

2. Această metodă de a estima e este discutată în lucrarea: Russell, K.G. Estimating the value of \(e\) by simulation, The American Statistician, Vol. 45, Nr. 1, pp 66-68, 1991.

3. Pentru mai multe detalii se poate consulta pagina https://en.wikipedia.org/wiki/Wigner_semicircle_distribution.