Stimatori distorti e non distorti della varianza

Dividere per n - 1.

Ci sono due formule per la varianza:

\[ S^2 = \frac{\sum_{i = 1}^n (X_i - \bar{X})^2}{n} \]

e

\[ s^2 = \frac{\sum_{i = 1}^n (X_i - \bar{X})^2}{n - 1} \] La prima formula è quella di una statistica, e si pone il problema di fornire descrizione sintetica di una proprietà del campione – in questo caso, la varianza.

La seconda formula è quella di uno stimatore, e si pone l’obiettivo di descrivere, nella maniera più fedele possibile, una proprietà della popolazione – in questo caso, la varianza – utilizzando le informazioni presenti nel campione.

La prima formula, quella con \(n\) al denominatore, fallisce nello scopo che la seconda formula si propone di raggiungere (ovvero, l’obiettivo di descrivere le proprietà della popolazione). Infatti, in media, produrrà una stima troppo piccola. Usiamo una simulazione per esaminare questa faccenda.

Iniziamo a definire le proprietà della popolazione.

mu <- 100 
sigma <- 15

Decidiamo, inoltre, di considerare campioni di ampiezza \(n\) = 30.

sample_size <- 30

Nel caso di un singolo campione, per esempio, abiamo:

one_sample <- rnorm(sample_size, mu, sigma) 
mean(one_sample)

[1] 100.8664

var(one_sample)

[1] 192.2824

Ma un singolo campione ci dice poco delle caratteristiche della “formula” che stiamo esaminando – quella che ha \(n\) al denominatore. Dato che è facile farlo con R, esaminiamo quello che succede quando consideriamo un milione di campioni:

var_distr <- replicate(
  1e6,
    var(
      rnorm(sample_size, mu, sigma) 
    ) * (sample_size - 1) / sample_size
)

La funzione rnorm() estrae un campione casuale di ampiezza sample_size (ovvero, 30) da una popolazione normale di media 100 e deviazione standard 15. La varianza della popolazione è dunque

15^2

[1] 225

La funzione var() * (sample_size - 1) / sample_size calcola la varianza delle 30 osservazioni di un campione utilizzando la prima formula (con \(n\) al denominatore).

La funzione replicate() ripete un milione di volte questi calcoli, ovvero calcola la varianza di un milione di campioni casuali di ampiezza 30 estratti da una popolazione normale di media 100 e varianza \(15^2\). Ciò significa che l’oggetto var_distr conterrà un milione di numeri: le varianze calcolate con la prima formula, per un milione di campioni casuali estratti dalla popolazione di riferimento.

Quanto bene ha funzionato la prima formula? Ovviamente, alcune volte (cioé, per alcuni campioni) meglio, altre volte peggio. Il valore più piccolo che abbiamo ottenuto è

min(var_distr)

[1] 44.00988

e il valore più grande che abbiamo ottenuto è

max(var_distr)

[1] 604.1064

Il valore medio – ovvero, il valore atteso del valore che si ottiene utilizzando la prima formula, è

mean(var_distr)

[1] 217.4179

Questo valore è chiaramente troppo piccolo. In altre parole, la prima formula, se venisse usata per stimare la varianza della popolazione produrrebbe una sottostima del valore del parametro cercato.

Si può correggere questo errore sistematico?

Certamente! E questo è lo scopo della seconda formula, quella con \(n - 1\) al denominatore. Ripetiamo la simulazione usando la seconda formula:

var_distr_c <- replicate(
  1e6,
  var(
    rnorm(sample_size, mu, sigma) 
  )
)

In questo caso, il valore atteso è

mean(var_distr_c)

[1] 225.0316

uguale al valore del parametro

15^2

[1] 225

Nella simulazione il risultato non è perfetto, ma si capisce che questa è, appunto, una simulazione: aumentando il numero delle ripetizioni si otterrebbe un valore sempre più simile al valore teorico. Ma non è necessario fare questo. La conclusione è chiara: la seconda formula ci fornisce uno stimatore corretto (ovvero, privo di errore sistematico) della varianza della popolazione. Per questa ragione dividiamo per (\(n\) - 1).