AJO_Resumen_estimación_parámetros

110 days ago by sevillad1819

Estimación de los parámetros estadísticos de una distribución normal

Supongamos que la variable $X$, que toma los valores resultantes de un cierto proceso de medición, sigue una distribución normal $N(\mu,\sigma)$. Casi siempre $\mu$ es desconocido (es el valor real o ideal a medir). En algunas situaciones conoceremos $\sigma$ (por las especificaciones del aparato de medida, por ejemplo), en otras no. Queremos estimar los parámetros de $X$ a partir de un conjunto finito de valores (una muestra). En esta hoja ilustraremos varios casos por medio de ejemplos, que serán generados artificialmente a modo de ilustración. En lugar de ejecutarlos sin más, estudia las fórmulas teóricas y su contrapartida en Sage.

En los casos en que hagamos estimaciones por intervalo de confianza, por defecto trabajamos con probabilidad 95%. En concreto, esto significa que de cada 20 ejecuciones del código, en 19 de ellas el valor real está contenido en el intervalo de estimación; dicho de otra manera, 1 de cada 20 ese intervalo no contiene el valor real. Puedes experimentar tú mismo/a.

 
       

Estimación puntual de la media

La estimación puntual de $\mu$ que usamos es la media muestral, $\overline{x}$. A continuación generamos varios datos aleatoriamente a partir de una distribución normal arbitraria. Puedes ejecutar el código varias veces para ver cómo cambia la estimación. También puedes cambiar los parámetros.

mu = 234 sigma = 2 n = 5 T = RealDistribution('gaussian',sigma) d = [mu + T.get_random_element() for i in (1..n)] d 
       
[235.29162323007725,
 235.36152535527063,
 233.0016299775329,
 234.7278224295953,
 235.19126502620443]
[235.29162323007725,
 235.36152535527063,
 233.0016299775329,
 234.7278224295953,
 235.19126502620443]
mean(d) 
       
234.71477320373612
234.71477320373612
 
       

Estimación puntual de la varianza

La estimación puntual de $\sigma^2$ que usamos es la cuasivarianza muestral, $S^2$. La varianza muestral no es un buen indicador: se queda corto, en general, por un factor $\frac{n-1}{n}$ (compárense las fórmulas de la varianza y cuasivarianza muestrales).

A continuación un ejemplo de estimación puntual de la varianza, hecho "a mano" y también con la función de Sage correspondiente.

mu = 234 sigma = 2 n = 5 T = RealDistribution('gaussian',sigma) d = [mu + T.get_random_element() for i in (1..n)] d 
       
[233.68959549342722,
 232.76659985979904,
 233.1085630692332,
 235.13602108613654,
 238.7734166857128]
[233.68959549342722,
 232.76659985979904,
 233.1085630692332,
 235.13602108613654,
 238.7734166857128]
xbarra = mean(d) sum([(x_i-xbarra)^2 for x_i in d])/(n-1); variance(d) 
       
6.0185823973403725
6.0185823973403725
6.0185823973403725
6.0185823973403725
 
       

Estimación por intervalo de la media (caso de varianza conocida)

Supongamos que conocemos $\sigma^2$. Hemos visto que la variable \[ \overline{X}=\frac{X_1+\cdots+X_n}{n} \] sigue una distribución $N(\mu,\frac{\sigma}{\sqrt{n}})$. Estandarizando, la nueva variable \[ Y=\frac{\overline{X}-\mu}{\sigma/\sqrt{n}} \] sigue una distribución $N(0,1)$. Usando el valor de $\overline{X}$ dado por la muestra, del intervalo de confianza que calculemos para $Y$ despejamos $\mu$ y sacamos su intervalo (nótese que en el último paso al despejar una desigualdad del tipo $a \leq -\mu \leq b$ se cambia el orden y se obtiene $-b\leq \mu \leq -a$).

mu = 234 sigma = 2 n = 5 T = RealDistribution('gaussian',sigma) d = [mu + T.get_random_element() for i in (1..n)] d 
       
[233.72345428686285,
 232.63409099564663,
 232.0607831454154,
 235.6793816392266,
 236.59208981634254]
[233.72345428686285,
 232.63409099564663,
 232.0607831454154,
 235.6793816392266,
 236.59208981634254]
xbarra = mean(d) print "media", xbarra 
       
media 234.137959977
media 234.137959977
conf = 0.95 a_Y = r.qnorm((1-conf)/2)._sage_() b_Y = -a_Y print "cuantiles", (a_Y, b_Y) 
       
cuantiles (-1.95996398454005, 1.95996398454005)
cuantiles (-1.95996398454005, 1.95996398454005)
a_mu = -b_Y * (sigma/sqrt(n).n()) + xbarra b_mu = -a_Y * (sigma/sqrt(n).n()) + xbarra print "intervalo", (a_mu, b_mu) 
       
intervalo (232.384914895546, 235.891005057852)
intervalo (232.384914895546, 235.891005057852)
 
       

Estimación por intervalo de la media (caso de varianza desconocida)

Si no conocemos $\sigma^2$, no podemos despejar $\mu$ en la fórmula de arriba. Así que lo hacemos de manera ligeramente distinta: usando \[ Y=\frac{\overline{X}-\mu}{S/\sqrt{n}} \] donde ahora sustituimos la desviación típica $\sigma$, que nos es desconocida, por la estimación $S$, la raíz cuadrada de la nueva variable \[ S^2 = \frac{1}{n-1}\sum (X_i-\overline{X})^2 . \] Lo que cambia es que ahora la distribución de $Y$ no es la normal $N(0,1)$ sino una t de Student con $n-1$ grados de libertad. A modo de ilustración aquí tienes una comparación animada entre $N(0,1)$ y las $t$ de Student para distintos grados de libertad.

pnormal = plot(RealDistribution('gaussian', 1),xmin=-4,xmax=4,color="black",figsize=6) frames = [plot(RealDistribution('t', i),xmin=-4,xmax=4,color="blue",legend_label=str(i)+" grados de libertad")+pnormal for i in [1..20]] animate(frames).show(delay=20) 
       

De nuevo, calcularemos un intervalo de confianza para $Y$ con los valores de $\overline{X}$ y $S$ de la muestra, y despejaremos $\mu$.

mu = 234 sigma = 2 n = 5 T = RealDistribution('gaussian',sigma) d = [mu + T.get_random_element() for i in (1..n)] xbarra = mean(d) s2 = variance(d) print d print "media", xbarra print "varianza", s2 
       
[235.40691836556402, 235.54704848168086, 233.307234367067,
235.12956910872828, 230.23411808972972]
media 233.924977683
varianza 5.07058842678
[235.40691836556402, 235.54704848168086, 233.307234367067, 235.12956910872828, 230.23411808972972]
media 233.924977683
varianza 5.07058842678
conf = 0.95 a_Y = r.qt((1-conf)/2, n-1)._sage_() b_Y = -a_Y print "cuantiles", (a_Y, b_Y) 
       
cuantiles (-2.77644510519779, 2.77644510519779)
cuantiles (-2.77644510519779, 2.77644510519779)
a_mu = -b_Y * (sqrt(s2)/sqrt(n).n()) + xbarra b_mu = -a_Y * (sqrt(s2)/sqrt(n).n()) + xbarra print "intervalo", (a_mu, b_mu) 
       
intervalo (231.1290027754768, 236.72095258963117)
intervalo (231.1290027754768, 236.72095258963117)
 
       

Estimación por intervalo de la varianza

Usamos que la variable \[ Y=\frac{(n-1)S^2}{\sigma^2} \] sigue una distribución chi cuadrado ($\chi^2$) con $n-1$ grados de libertad. A continuación las gráficas de varias distribuciones chi cuadrado.

frames = [plot(RealDistribution('chisquared', i),xmin=0,xmax=20,ymax=1,aspect_ratio=10,color="blue",legend_label=str(i)+" grados de libertad") for i in [1..12]] animate(frames).show(delay=20) 
       

Del intervalo de confianza para $Y$ y el valor de $S^2$ obtenido de la muestra despejamos $\sigma^2$ y obtenemos su intervalo de confianza. A la hora de despejar usamos que una desigualdad del tipo $a \leq \frac{1}{\sigma^2} \leq b$ se transforma, por la inversa, en $\frac1b \leq \sigma^2 \leq \frac1a$.

mu = 234 sigma2 = 2 n = 5 T = RealDistribution('gaussian',sigma) d = [mu + T.get_random_element() for i in (1..n)] s2 = variance(d) print "varianza", s2 
       
varianza 1.78203929599
varianza 1.78203929599
conf = 0.95 a_Y = r.qchisq((1-conf)/2, n-1)._sage_() b_Y = r.qchisq(1-(1-conf)/2, n-1)._sage_() print "cuantiles", (a_Y, b_Y) 
       
cuantiles (0.48441855708793, 11.1432867818778)
cuantiles (0.48441855708793, 11.1432867818778)
a_sigma2 = (n-1)*s2/b_Y b_sigma2 = (n-1)*s2/a_Y print "intervalo", (a_sigma2, b_sigma2) 
       
intervalo (0.6396817495135354, 14.714872251825739)
intervalo (0.6396817495135354, 14.714872251825739)
 
       

Estimación por intervalo del cociente de varianzas

Tenemos dos variables $X,Y$ que siguen distribuciones normales, cuyas varianzas son $\sigma_X^2$ y $\sigma_Y^2$ respectivamente. Usamos que la variable \[ Z=\frac{S_X^2/\sigma_X^2}{S_Y^2/\sigma_Y^2} \] sigue una distribución F de Snedecor con $n_X-1$ y $n_Y-1$ grados de libertad (no necesitamos saber nada sobre la distribución F, salvo el comando de Sage que calcula sus cuantiles, y que se usa debajo). Del intervalo de confianza para $Z$ y los valores muestrales de $S_X^2,S_Y^2$ despejamos $\sigma_X^2/\sigma_Y^2$ y obtenemos su intervalo de confianza de manera parecida al caso de arriba.

mu_X = 234 sigma_X = 2 n_X = 50 T_X = RealDistribution('gaussian',sigma_X) d_X = [mu_X + T_X.get_random_element() for i in (1..n_X)] mu_Y = 157 sigma_Y = 3 n_Y = 10 T_Y = RealDistribution('gaussian',sigma_Y) d_Y = [mu_Y + T_Y.get_random_element() for i in (1..n_Y)] s2_X = variance(d_X) s2_Y = variance(d_Y) print "X n, var, var_muestral:", n_X, sigma_X^2, s2_X print "Y n, var, var_muestral:", n_Y, sigma_Y^2, s2_Y print "cociente", sigma_X/sigma_Y.n() 
       
X n, var, var_muestral: 50 4 3.7124010477
Y n, var, var_muestral: 10 9 9.88696768183
cociente 0.666666666666667
X n, var, var_muestral: 50 4 3.7124010477
Y n, var, var_muestral: 10 9 9.88696768183
cociente 0.666666666666667
conf = 0.95 a_Z = r.qf((1-conf)/2, n_X-1, n_Y-1)._sage_() b_Z = r.qf(1-(1-conf)/2, n_X-1, n_Y-1)._sage_() print "cuantiles", (a_Z, b_Z) 
       
cuantiles (0.419013382865503, 3.47468082184094)
cuantiles (0.419013382865503, 3.47468082184094)
a_cociente = s2_X/s2_Y * 1/b_Z b_cociente = s2_X/s2_Y * 1/a_Z print "intervalo", (a_cociente, b_cociente) print "cociente", sigma_X/sigma_Y.n() 
       
intervalo (0.10806295887439013, 0.8961152700765991)
cociente 0.666666666666667
intervalo (0.10806295887439013, 0.8961152700765991)
cociente 0.666666666666667
 
       

Estimación por intervalo de la diferencia de medias (caso de observaciones independientes y varianzas conocidas)

Tenemos dos variables $X,Y$ que siguen distribuciones normales, con medias $\mu_X$ y $\mu_Y$ y varianzas $\sigma_X^2$ y $\sigma_Y^2$ respectivamente. Como la variable $\overline{X} - \overline{Y}$ sigue una distribución normal de media $\mu_X - \mu_Y$ y varianza $\frac{\sigma_X^2}{n_X} + \frac{\sigma_Y^2}{n_Y}$, tenemos que \[ Z = \frac{(\overline{X} - \overline{Y}) - (\mu_X - \mu_Y)}{\sqrt{\sigma_X^2/n_X + \sigma_Y^2/n_Y}} \] sigue una distribución $N(0,1)$. Del intervalo de confianza para $Z$ y los distintos valores muestrales despejamos $\mu_X - \mu_Y$ y obtenemos su intervalo de confianza.

mu_X = 234 sigma_X = 2 n_X = 50 T_X = RealDistribution('gaussian',sigma_X) d_X = [mu_X + T_X.get_random_element() for i in (1..n_X)] mu_Y = 157 sigma_Y = 3 n_Y = 10 T_Y = RealDistribution('gaussian',sigma_Y) d_Y = [mu_Y + T_Y.get_random_element() for i in (1..n_Y)] xbarra = mean(d_X) ybarra = mean(d_Y) s2_X = variance(d_X) s2_Y = variance(d_Y) print "X n, mu, xbarra, var_muestral:", n_X, mu_X, xbarra, s2_X print "Y n, mu, ybarra, var_muestral:", n_Y, mu_Y, ybarra, s2_Y print "diferencia", mu_X-mu_Y 
       
X n, mu, xbarra, var_muestral: 50 234 233.690181139 2.35391748324
Y n, mu, ybarra, var_muestral: 10 157 155.791291933 10.5129035197
diferencia 77
X n, mu, xbarra, var_muestral: 50 234 233.690181139 2.35391748324
Y n, mu, ybarra, var_muestral: 10 157 155.791291933 10.5129035197
diferencia 77
conf = 0.95 a_Z = r.qnorm((1-conf)/2)._sage_() b_Z = -a_Z print "cuantiles", (a_Z, b_Z) 
       
cuantiles (-1.95996398454005, 1.95996398454005)
cuantiles (-1.95996398454005, 1.95996398454005)
denom = sqrt(s2_X/n_X + s2_Y/n_Y) a_diferencia = -b_Z*denom + (xbarra-ybarra) b_diferencia = -a_Z*denom + (xbarra-ybarra) print "intervalo", (a_diferencia, b_diferencia) print "diferencia", mu_X-mu_Y 
       
intervalo (75.84478644864954, 79.95299196311974)
diferencia 77
intervalo (75.84478644864954, 79.95299196311974)
diferencia 77