my_integral <- function(f, g, R = 1000) {
mean(f(rnorm(R)) / g(rnorm(R)))
}R per l’analisi statistica multivariata
Esame 20 Febbraio 2023
Problema 1
(3pt) Calcolare una stima via simulazione Monte Carlo del seguente integrale I_1 = \int_0^1\int_0^1 \exp{(\sqrt{x y})}\mathrm{d}x\mathrm{d}y,
(4pt) Calcolare una stima via simulazione Monte Carlo del seguente integrale I_2 = \int_0^2 \exp{(\sqrt{x})}\mathrm{d}x, valutandone in qualche maniera anche la precisione.
(3pt) In un appello passato avevo chiesto agli studenti di scrivere una funzione
my_integral(f, g, R = 1000)che calcoli un’approssimazione Monte Carlo di I_3 = \int_{-\infty}^{+\infty} \frac{f(x)}{g(x)} \phi(x)\mathrm{d}x, dove \phi(x) è la densità di una normale standard. Gli argomenti dimy_integralindicano le funzioni f(x), g(x) ed il numero di simulazioni R su cui basare la stima. Uno studente ha risposto in questo modo:
Si spieghi come mai questa risposta non è corretta e si fornisca quindi la soluzione.
Problema 2
Si consideri il dataset di R che si ottiene tramite il comando data(faithful) all’interno della libreria MASS. Se ne consulti la documentazione per ulteriori informazioni.
(1pt) Si ottenga un istogramma della variabile
waiting. Si ottenga inoltre un diagramma a dispersione per le variabilieruptions(x) ewaiting(y).(2pt) Si creino le variabili
long_waitingeshort_waiting, contenenti rispettivamente i valori della variabilewaitingrelativi ad eruzioni lunghe (eruptions > 3) ed eruzioni brevi (eruptions <= 3).(3pt) Si confrontino le funzioni di ripartizioni empiriche delle variabili
long_waitingeshort_waiting. Inoltre, si calcolino media e mediana delle variabililong_waitingeshort_waiting.(3pt) Si scriva la funzione
asym(x)che calcola il coefficiente di asimmetria secondo Bowley, definito come B = \frac{\mathcal{Q}_{0.75} - 2 \mathcal{Q}_{0.5} + \mathcal{Q}_{0.25}}{\mathcal{Q}_{0.75} - \mathcal{Q}_{0.25}}, dove \mathcal{Q}_p rappresenta il quantile p-esimo dei dati. Si calcoli infine il coefficiente di asimmetria secondo Bowley per le variabililong_waitingeshort_waiting.
Problema 3
Il test di Jarque-Bera è un test molto usato in ambito econometrico per verificare l’ipotesi di normalità. Si supponga di disporre di dati x_1,\dots,x_n, che assumiamo siano determinazioni indipendenti ed identicamente distribuite da una variabile aleatoria X. Il sistema d’ipotesi è il seguente H_0 : \{X \sim \text{N}(\mu, \sigma^2), \text{ per } \mu, \sigma^2 \text{ qualsiasi} \}, \quad H_1: \{\text{La distribuzione di } X \text{ non è normale}\}.
La statistica test è \text{JB} = \frac{n}{6}\left(\gamma^2 + \frac{(\kappa - 3)^2}{4}\right), dove \gamma e \kappa sono gli indici di asimmetria e curtosi di Pearson, ovvero:
\gamma = \frac{1}{\text{sqm}(x)^3} \frac{1}{n}\sum_{i=1}^n(x_i - \bar{x})^3, \quad \kappa = \frac{1}{\text{sqm}(x)^4} \frac{1}{n}\sum_{i=1}^n(x_i - \bar{x})^4,
per dei dati x_1,\dots,x_n aventi media \bar{x}. Per calcolare lo scarto quadratico medio \text{sqm(x)} si faccia uso della funzione sd.
Il test rifiuta l’ipotesi nulla per valori grandi di JB, ovvero quando \text{JB} > h_{1-\alpha}, dove \alpha indica il livello di significatività mentre h_p il quantile p-esimo della distribuzione di quando è vera H_0.
- (10pt) Si scriva una funzione
qJB(p, n, R = 1000)che calcoli, via simulazione, i quantili p-esimi h_p della distribuzione della statistica test quando è vera H_0. In particolare, la funzione deve:- Generare R campioni di lunghezza n da una normale standard;
- Calcolare la statistica per ogni campione simulato;
- Restituire i quantili campionari dei valori di JB ottenuti al passo precedente.
- (2pt) Calcolare h_{0.95} quando n = 50.