Metodi per la costruzione della distribuzione approssimata di una statistica: Monte Carlo

Alla fine di questa puntata…

• dovreste avere imparato le funzioni disponibili in R base per simulare campioni per l’utilizzo del metodo Monte Carlo
• dovreste avere imparato le funzioni disponibili in R base per estrarre campioni per l’utilizzo del metodo Bootstrap
• dovreste avere imparato le funzione disponibili in R base (funzioni della classe `apply`) per effettuare calcoli su un oggetto di tipo lista
• dovreste avere imparato le funzioni disponibili nel package purrr (installato e caricato automaticamente con tidyverse) per effettuare calcoli usando un approccio di programmazione funzionale
• dovreste essere in grado di costruire la distribuzione campionaria approssimata di una statistica con il metodo Monte Carlo
• dovreste essere in grado di costruire la distribuzione campionaria approssimata di una statistica con il metodo Bootstrap

Prima di cominciare… prerequisiti

• Sul tuo PC ipotizzo che tu abbia installato sia R che RStudio. Se ti servono istruzioni al riguardo consulta la pagina con le istruzioni disponibile a questo link.

• Sul tuo PC ipotizzo che tu abbia installato anche la collezione di package tidyverse e il package ineq

Cosa utilizzeremo in questa lezione?

Le funzioni sfruttate in queste lezioni sono disponibili:

• nell’installazione base di R
• nel package purrr (installato automaticamente quando si installa tidyverse)
• nel package ineq che contiene le principali misure di disguaglianza: sfrutteremo la funzione Gini per il calcolo del coefficiente di concentrazione di Gini

Il metodo Monte Carlo

E’ possibile ottenere la costruzione della distribuzione campionaria approssimata di una statistica con il metodo Monte Carlo utilizzando poche righe di codice e sfruttando le funzioni disponibili con l’installazione base di R.

La procedura Monte Carlo può essere sintetizzata nei seguenti passi:

1. Scelta di una pseudo-popolazione (modello generatore dei dati, DGP - Data Generating Process)
2. Estrazione di $M$ campioni casuali dalla pseudo-popolazione
3. Calcolo della statistica $T_n$ di interesse sugli $M$ campioni ottenuti al passo precedente

I valori della statistica \(T_n\) ottenuti al passo 3. della procedura sopra riportata costituiscono una distribuzione approssimata Monte Carlo della statistica.

Possiamo mostrare il funzionamento nel caso di una statistica \(T_n\) di cui conosciamo dalla teoria il funzionamento in modo da verificare empiricamente quello che abbiamo già studiato. Scegliamo a tal fine come primo esempio la statistica media campionaria \(\bar{X}_n\), simulando dati da un modello generatore normale e da un modello generatore non normale: nel primo caso sappiamo che la media campionaria segue esattamente una distribuzione normale con media \(\mu\) e varianza \(\frac{\sigma^2}{n}\) per la proprietà riproduttiva della v.c. normale, nel secondo caso una distribuzione asintotica con gli stessi parametri per il teorema del limite centrale.

Fissiamo i valori di \(M\) (numero di campioni Monte Carlo da generare), \(n\) (ampiezza campionaria), \(\mu\) e \(\sigma^2\) (parametri del modello generatore, da cui possiamo ottenere dalla teoria i corrispondenti parametri della statistica media campionaria):

# numero di campioni e taglia del singolo campione
M <- 1000      # numero di replicazioni Monte Carlo
n <- 50        # ampiezza singolocampione
#
# parametri modello generatore (pseudo-popolazione)
media <- 100
sqm <- 10

Quando si lavora con i numeri (pseudo-)casuali, è conveniente fissare il seme dell’algoritmo di generazione ai fini della riproducibilità. A parità di seme otterremo cioè sempre gli stessi numeri sia su successivi utilizzi del codice che su utilizzi dello stesso codice su diverse macchine:

# scegliamo un numero come seme dell'algoritmo di generazione dei numeri casuali
set.seed(12)

Media campionaria nel caso di pseudo-popolazione normale

La funzione replicate permette di replicare un’espressione (secondo argomento) per un numero fissato di volte (primo argomento). Stiamo sostanzialmente utilizzando un ciclo di iterazione lasciando a R la complessità della gestione dell’output, in quando non dobbiamo preoccuparci di alcuna inizializzazione, iterazione e puntamento alle posizioni opportune per memorizzare i risultati delle singole iterazioni:

campioni_mc <- replicate(M, rnorm(n = n, mean = media, sd = sqm))

La funzione replicate, per default, cerca di semplificare la struttura di output in un vettore o in un array, a seconda del tipo di funzione. Per default infatti l’argomento simplify, che regola tale semplificazione, è impostato a TRUE:

campioni_mc_array <- replicate(M, rnorm(n = n, mean = media, sd = sqm), 
                               simplify = TRUE)

Possiamo andare ad esplorare la classe della struttura di output usando la funzione class:

class(campioni_mc_array)

[1] "matrix" "array" 

ed andando a calcolare la dimensione della matrice ottenuta ci accorgiamo subito che in questo modo otteniamo una matrice dove in ciascuna colonna è presente uno degli \(M\) campioni Monte Carlo:

dim(campioni_mc_array)

[1]   50 1000

Per lavorare su uno dei singoli campioni possiamo usare la sintassi standard che ci permette di estrarre una singola colonna:

# ecco il primo degli M campioni Monte Carlo
campioni_mc_array[, 1]

 [1] 103.50159 100.93605 104.46844  87.30646 102.63608  85.23516 105.20665
 [8] 113.65104 107.85544  95.05632 111.23594 109.42525 112.04507  90.59259
[15] 108.02223 102.91878  97.61508 108.31631 107.50594  90.98721  93.43796
[22] 101.61035  77.26247 117.04986 109.59251 115.39315 117.45238 104.45735
[29] 110.55206 117.69836 109.82019  98.93695 104.18338 110.98241 101.36923
[36] 103.60483 102.66435 100.63068 107.65729  81.12922 113.13717  96.96738
[43] 108.65947  97.45926 100.41171 108.75260 110.08252  99.13132  85.82960
[50] 112.29631

Se si imposta invece l’argomento simplify a FALSE si ottiene invece un oggetto di tipo lista:

campioni_mc_lista <- replicate(M, rnorm(n = n, mean = media, sd = sqm), 
                               simplify = FALSE)
class(campioni_mc_lista)

[1] "list"

che ha lunghezza pari ad \(M\):

length(campioni_mc_lista)

[1] 1000

Ciascun elemento della lista contiene uno degli \(M\) campioni simulati dalla pseudo-popolazione:

# ecco il primo degli M campioni Monte Carlo
campioni_mc_lista[[1]]

 [1]  96.20926 100.01220 113.05492  85.67140 103.68631 112.01379  94.75258
 [8]  89.34047 110.29806 102.96797 103.17488 107.21611  91.02462 101.86665
[15] 111.89525 109.49269 100.25497 104.75443  92.17095 112.15244 106.99232
[22] 101.96798 112.61655  99.53630 105.13423  93.30923 111.51515 107.20495
[29] 102.88861 103.41417  94.56612 110.64814  97.63495  97.97171 100.85683
[36] 105.41704  88.02295  89.92695  99.06479 109.17660 105.83947 102.30404
[43] 107.44457 103.96636 101.90229 104.17690 106.89087  93.38691  95.67662
[50]  97.96788

Anche se apparentemente la struttura di tipo lista è più complessa di una struttura di tipo array, si rivela molto comoda ai fini delle elaborazioni.

Per il calcolo su strutture dati di tipo lista si può sfruttare la funzione lapply che “applica” una funzione (secondo argomento) su una lista (primo argomento). La funzione restituisce una lista che ha pari lunghezza, ciascun elemento conterrà la media del corrispondente campione Monte Carlo:

# calcolo le medie sulla lista
medie_mc <- lapply(campioni_mc_lista, mean)
# ottengo in output un oggetto di tipo lista
class(medie_mc)

[1] "list"

# di pari lunghezza
length(medie_mc)

[1] 1000

Si può sfruttare anche la variante sapply che ha la stessa sintassi ma differisce in quanto restituisce (se il tipo di output calcolato lo consente) un vettore o un array:

# calcolo le medie sulla lista chiedendo di semplificare la struttura di output
medie_mc <- sapply(campioni_mc_lista, mean)
# ottengo in output un oggetto di tipo array
class(medie_mc)

[1] "numeric"

Su questo oggetto semplificato è semplice calcolare quello che ci interessa, ovvero esplorare le caratteristiche della distribuzione Monte Carlo della statistica ottenuta:

# media della distribuzione campionaria della media
mean(medie_mc)

[1] 100.0319

# media della popolazione
media

[1] 100

#
# scarto quadratico medio della distribuzione campionaria della media
sd(medie_mc)

[1] 1.393886

# scarto quadratico medio della popolazione / radice di n
sqm / sqrt(n)

[1] 1.414214

hist(medie_mc, freq = FALSE)
x <- seq(90, 110, 0.1)
y <- dnorm(x, mean = media, sd = sqm / sqrt(n))
lines(x, y, col = "red", lwd = 2)

boxplot(medie_mc, horizontal = TRUE)

qqnorm(medie_mc)

Media campionaria nel caso di pseudo-popolazione non normale (Poisson)

campioni_mc <- replicate(M, 
                         rpois(n = n, lambda = media), 
                         simplify = FALSE)

medie_mc <- sapply(campioni_mc, mean)
mean(medie_mc)

[1] 99.96232

sd(medie_mc)

[1] 1.449544

Media campionaria nel caso di pseudo-popolazione non normale (chi-quadro)

campioni_mc <- replicate(M, 
                         rchisq(n = n, df = media), 
                         simplify = FALSE)

medie_mc <- sapply(campioni_mc, mean)
mean(medie_mc)

[1] 99.96456

sd(medie_mc)

[1] 1.955139

sqrt(200) / sqrt(50)

[1] 2

Coefficiente di concentrazione di Gini nel caso di pseudo-popolazione non normale (chi-quadro)

# DGP: Chi-quadro
library(ineq)
campioni_mc <- replicate(M, 
                         rchisq(n = n, df = 30), 
                         simplify = FALSE)

hist(campioni_mc[[1]])

ineq(campioni_mc[[1]], type = "Gini")

[1] 0.1303027

gini_mc <- sapply(campioni_mc, Gini)
hist(gini_mc)

Nota tecnica: le funzione “moderne” (map e dintorni) del package purrr

library(purrr)

out_map <- map(campioni_mc, mean)
class(out_map)

[1] "list"

out_map_vec <- map_vec(campioni_mc, mean)
class(out_map_vec)

[1] "numeric"

set.seed(12)
campioni_mc <- map(1:M, ~rnorm(n = n, mean = media, sd = sqm))
class(campioni_mc)

[1] "list"

length(campioni_mc)

[1] 1000

out_map <- map(campioni_mc, mean)
out_map_vec <- map_vec(campioni_mc, mean)

mean(out_map_vec)

[1] 100.0227

sd(out_map_vec)

[1] 1.412397