Proiect de laborator

Mai jos regăsiți subiectele pentru proiectul de laborator pentru cursul de Statistică.

Precizări importante

1. Fiecare echipa (care conține maxim 2 persoane) va trimite prin liderul echipei un singur e-mail la adresa alexandru.amarioarei@fmi.unibuc.ro ce va conține proiectul. De asemenea am să fac o temă (un assignment) pe Teams și vă rog să încărcați proiectele și acolo (toată lumea, chiar dacă se vor repeta). Documentația proiectului trebuie scrisă prin intermediul pachetului Rmarkdown/Quarto din R. Informații introductive despre modul de folosire a acestui pachet pot fi găsite aici sau aici iar pentru mai multe detalii se poate consulta cartea R Markdown: The definitive guide. Toate simulările, figurile și codurile folosite trebuie incluse în raport. Se va folosi doar limbajul R. Documentația trimisă trebuie să conțină pe lângă fișierul generat (HTML, Microsoft Word sau \(\LaTeX\)) și fișierul .Rmd care conține codul sursă (comentat!).

2. Documentația proiectului trebuie să conțină, pe prima pagină, numele membrilor echipei, liderul echipei și grupa din care face parte fiecare. Pentru fiecare exercițiu în parte documentația trebuie să conțină:

   • Calculele matematice solicitate (pot fi tehnoredactate sau scrise de mână și scanate, însă dacă apelați la ultima variantă vă rog să vă asigurați că ați scris citeț, fără ștersături sau zone scanate deficitar)
   • Codul R comentat
   • Graficele realizate
   • Comentariile și concluziile voastre

3. Pentru obținerea punctajului maxim trebuie să rezolvați toate subiectele.

Important

Data de predare a proiectului este 4 februarie 2024 ora 22:00.

Subiect

Fie \(X\) și \(Y\) două variabile aleatoare reale independente astfel încât:

• \(X\) este repartizată \(\Gamma(\alpha,\theta)\) cu densitatea

\[ f_{X}(x)=\frac{1}{\Gamma(\alpha) \theta^\alpha} x^{\alpha-1} \exp \left(-\frac{x}{\theta}\right), \, x\in\mathbb{R}^{*}_{+} \]

unde \(\alpha>0\) și \(\theta>0\) sunt parametrii cunoscuți

• \(Y\) admite densitatea de repartiție

\[ f_{Y}(y)=\frac{y-\mu+2 b}{4 b^2} \mathbb{1}_{\{y \in[\mu-2 b, \mu]\}}+\frac{1}{2 b} \exp \left(-\frac{y-\mu}{b}\right) \mathbb{1}_{\{y>\mu\}}, \, x\in\mathbb{R} \]

unde \(\mu\in\mathbb{R}\) și \(b>0\) sunt parametrii cunoscuți.

Obiectivul aceste probleme este de a estima probabilitatea

\[ \delta=\mathbb{P}[X+Y>t] \]

pentru un \(t\) dat. În cele ce urmează putem considera că \(\alpha = 3\), \(\theta = 2\), \(\mu = -13\), \(b = 2\) și \(t = 0\) dar calculele pot fi făcute pentru cazul general.

Metoda 1 - Monte Carlo clasic (naiv)

1. Calculați exact probabilitatea \(\delta\). Comparați cu rezultatul obținut prin aplicarea funcției integrate.

2. Descrieți o metodă de simulare a unor observații din repartiția vectorului \((X,Y)\) și scrieți codul R corespunzător. Precizați eventualele calcule efectuate.

3. Verificați grafic (e.g. prin intermediul unor histograme) că metoda propusă de simulare și codul sunt corecte.

4. Descrieți estimatorul Monte Carlo naiv \(\hat{p}_n\) a lui \(\delta\) folosind un eșantion \(\left(X_k, Y_k\right)_{k \geq 1}\) i.i.d. de volum \(n\) din repartiția vectorului \((X,Y)\).

5. Scrieți un cod R care să returneze pentru o valoare a lui \(n\) (e.g. \(n = 10000\)) valoarea estimatorului \(\hat{p}_n\) și eroarea pătratică medie a lui \(\hat{p}_n\).

6. Notăm cu \(F_X\) funcția de repartiție a lui \(X\). Arătați că

\[ \delta=\mathbb{E}\left[1-F_X(t-Y)\right], \] și găsiți un estimator \(\hat{\delta}_n\) a lui \(\delta\) plecând de la un eșantion \(U_1, U_2,\ldots, U_n\sim\mathcal{U}([0,1])\) de volum \(n\).

7. Implementați un cod R care să returneze pentru o valoare a lui \(n\) (e.g. \(n = 10000\)) valoarea estimatorului \(\hat{\delta}_n\) și eroarea pătratică medie a lui \(\hat{\delta}_n\).

8. Plecând de la aceste rezultate, putem spune că estimatorul \(\hat{\delta}_n\) este mai bun decât estimatorul \(\hat{p}_n\) în termeni de varianță? Dacă da, atunci arătați că are loc inegalitatea \(\operatorname{Var}\left[\widehat{\delta}_n\right] \leq \operatorname{Var}\left[\widehat{p}_n\right]\)

Metoda 2 - variabile antitetice

Metoda bazată pe variabile antitetice are ca scop propunerea unui estimator a cărui varianță este mai redusă decât cea a estimatorului Monte Carlo naiv. Metoda ține cont de simetria problemei. Astfel, dacă \(X\) este o variabilă aleatoare (sau vector aleator), \(h:\mathbb{R}(\mathbb{R}^d)\to\mathbb{R}\) o funcție măsurabilă iar scopul este de a estima \(\mathbb{E}[h(X)]\) atunci dacă există o transformare măsurabilă \(T:\mathbb{R}(\mathbb{R}^d)\to\mathbb{R}(\mathbb{R}^d)\) pentru care \(T(X)\) are aceeași repartiție cu \(X\) (\(T(X)\) se numește variabilă antitetică pentru \(X\)), estimatorul obținut prin metoda variabilelor antitetice este:

\[ \hat{\delta}_n^{A} = \frac{1}{n}\sum_{i=1}^{n}\frac{h(X_i) + h(T(X_i))}{2} \]

unde \(X_1,\ldots,X_n\) sunt variabile aleatoare i.i.d. cu repartiția lui \(X\). Se poate arăta cu ușurință că estimatorul astfel obținut este nedeplasat și are varianța mai mică decât a estimatorului Monte Carlo naiv dacă și numai dacă \(\operatorname{Cov}(h(X),h(T(X)))\leq 0\)¹. Pentru mai multe detalii referitoare la metoda variabilelor antitetice se poate consulta (Ross 2022, capitolul 9).

9. Folosind metoda variabilelor antitetice, propuneți un estimator \(\hat{\delta}_n^{A}\) a lui \(\delta\). Arătați că \(\operatorname{Var}\left[\hat{\delta}_n^{A}\right] \leq \operatorname{Var}\left[\hat{\delta}_n\right]\).

10. Implementați în R un cod care să returneze pentru o valoare a lui \(n\) (e.g. \(n = 10000\)) valoarea estimatorului \(\hat{\delta}_n^{A}\) și eroarea pătratică medie a lui \(\hat{\delta}_n^{A}\).

Metoda 3 - variabile de control

O altă metodă de reducere a varianței unui estimator este metoda bazată pe o variabilă de control. Metoda presupune găsirea unei funcții \(g\) mai simplă astfel încât \(f\approx g\) și pentru care să putem calcula analitic \(\mathbb{E}[g(X)]\). Astfel, dacă \(g:\mathbb{R}(\mathbb{R}^d)\to\mathbb{R}\) este o funcție măsurabilă pentru care \(\mathbb{E}[g(X)]=m\) este ușor de calculat și \(\operatorname{Var}(g(X))<\infty\) atunci definim estimatorul obținut prin metoda variabilei de control ca fiind

\[ \widehat{\delta}_n(b)=\frac{1}{n} \sum_{i=1}^n\left[h\left(X_i\right)-b\left(g\left(X_i\right)-m\right)\right] \]

unde \(X_1,\ldots,X_n\) sunt variabile aleatoare i.i.d. cu repartiția lui \(X\), iar \(b\) este un parametru real. Se poate arăta cu ușurință că estimatorul obținut este nedeplasat și are varianța mai mică decât a estimatorului Monte Carlo clasic dacă variabilele \(h(X)\) și \(g(X)\) sunt necorelate. Pentru mai multe detalii referitoare la metoda variabilelor antitetice se poate consulta (Ross 2022, capitolul 9).

Notăm cu \(x_{p}\) cuantila de ordin \(p\) a repartiției variabilei aleatoare \(X\) și definim funcțiile \(h_1, h_2:\mathbb{R}\to\mathbb{R}\) prin \(h_1(y) = y\) și respectiv \(h_2(y) = \mathbb{1}_{\left\{y \geq t-x_p\right\}}\). În cele ce urmează considerăm \(p = 0.6\).

11. Între \(h_1\) și \(h_2\) ce funcție este de preferat pentru a reduce varianța estimatorului \(\hat{\delta}_n\) prin intermediul metodei variabilei de control? Notăm acest estimator cu \(\hat{\delta}_n^{c}\).

12. Implementați un cod R care să returneze pentru o valoare a lui \(n\) (e.g. \(n = 10000\)) valoarea estimatorului \(\hat{\delta}_n^{c}\) și eroarea pătratică medie a lui \(\hat{\delta}_n^{c}\).

Pentru un șir \(\left(Y_k\right)_{k \geq 1}\) de variabile aleatoare i.i.d. din repartiția \(f_Y\), ne propunem să construim cel mai bun estimator de forma

\[ \widehat{\delta}_n(\beta)=1-\frac{1}{n} \sum_{k=1}^n F\left(t-Y_k\right)-<\beta,\binom{Y_k-\mathbb{E}[Y]}{\mathbb{1}_{\left\{Y_k \geq t-x_p\right\}}-\mathbb{P}\left[Y \geq t-x_p\right]}>, \quad \text { cu } \quad \beta=\binom{\beta_1}{\beta_2} \in \mathbb{R}^2 \]

unde \(<\cdot,\cdot>\) este produsul scalar din \(\mathbb{R}^2\), i.e. \(<u,v>=u^{\intercal} v\).

13. Arătați că

\[ \operatorname{Var}\left[\widehat{\delta}_n(\beta)\right]=\operatorname{Var}\left[\widehat{\delta}_n\right]-\frac{1}{n}(2<\beta, C>-<\beta, \Sigma \beta>) \]

unde

\[ C=\binom{\operatorname{Cov}[F_X(t-Y), Y]}{\operatorname{Cov}\left[F_X(t-Y), \mathbb{1}_{\left\{Y \geq t-x_p\right\}}\right]} \quad \text { și } \quad \Sigma=\left(\begin{array}{cc} \operatorname{Var}[Y] & \operatorname{Cov}\left[Y, \mathbb{1}_{\left\{Y \geq t-x_p\right\}}\right] \\ \operatorname{Cov}\left[Y, \mathbb{1}_{\left\{Y \geq t-x_p\right\}}\right] & \operatorname{Var}\left[\mathbb{1}_{\left\{Y \geq t-x_p\right\}}\right] \end{array}\right) . \]

14. Calculați explicit \(\Sigma\) și verificați dacă este inversabilă.

15. Demonstrați că

\[ \beta^{\star}=\underset{\beta \in \mathbb{R}^2}{\operatorname{argmin}} \operatorname{Var}\left[\widehat{\delta}_n(\beta)\right]=\Sigma^{-1} C \quad \text { și } \quad \operatorname{Var}\left[\widehat{\delta}_n\left(\beta^{\star}\right)\right] \leq \operatorname{Var}\left[\widehat{\delta}_n\right] . \]

16. Implementați un cod R care să returneze pentru o valoare a lui \(n\) (e.g. \(n = 10000\)) valoarea estimatorului \(\widehat{\delta}_n(\beta^{\star})\) și eroarea pătratică medie a lui \(\widehat{\delta}_n(\beta^{\star})\).

Metoda 4

17. Fie \(h\in L^2([0,1])\). Arătați că pentru un șir \(\left(U_k\right)_{k \geq 1}\) de variabile aleatoare i.i.d repartizate \(\mathcal{U}([0,1])\) are loc relația

\[ \operatorname{Var}\left[\frac{1}{n} \sum_{k=1}^n h\left(\frac{k-1+U_k}{n}\right)\right] \leq \operatorname{Var}\left[\frac{1}{n} \sum_{k=1}^n h\left(U_k\right)\right] \]

18. Deduceți un estimator \(\widehat{\Delta}_n\) a lui \(\delta\). Cum puteți interpreta acest estimator?

19. Implementați un cod R care să returneze pentru o valoare a lui \(n\) (e.g. \(n = 10000\)) valoarea estimatorului \(\widehat{\Delta}_n\) și eroarea pătratică medie a lui \(\widehat{\Delta}_n\).

20. Comparați eficacitatea relativă a estimatorilor lui \(\delta\) propuși în punctele anterioare.

Referințe

Ross, Sheldon M. 2022. Simulation. academic press.

Note de subsol

1. Atenție la volumul eșantionului: dacă \(\bar{\delta}_n\) este estimatorul Monte Carlo clasic a lui \(\delta\) atunci \(\operatorname{Cov}[h(X), h(T(X))]<0 \quad \Longleftrightarrow \quad \operatorname{Var}\left[\widehat{\delta}_n\right]<\frac{1}{2} \operatorname{Var}\left[\bar{\delta}_n\right]=\operatorname{Var}\left[\bar{\delta}_{2 n}\right].\)