AJO_1617_Resumen_estimación_parámetros

hace 206 días por sevillad1617

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 
       
[232.3365119270521,
 230.50460439912243,
 231.67242343856336,
 234.54219360684428,
 232.2568317595671]
[232.3365119270521,
 230.50460439912243,
 231.67242343856336,
 234.54219360684428,
 232.2568317595671]
mean(d) 
       
232.26251302622987
232.26251302622987
 
       

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 
       
[232.1210829121442,
 231.99340298656918,
 231.7042161971901,
 235.00915270824717,
 233.2840592795767]
[232.1210829121442,
 231.99340298656918,
 231.7042161971901,
 235.00915270824717,
 233.2840592795767]
xbarra = mean(d) sum([(x_i-xbarra)^2 for x_i in d])/(n-1) 
       
1.8561083547078603
1.8561083547078603
variance(d) 
       
1.8561083547078603
1.8561083547078603
 
       

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 
       
[232.18391510575853,
 232.51078528776378,
 235.39105146685287,
 233.247508394749,
 234.74800741915774]
[232.18391510575853,
 232.51078528776378,
 235.39105146685287,
 233.247508394749,
 234.74800741915774]
xbarra = mean(d) print "media", xbarra 
       
media 233.616253535
media 233.616253535
conf = 0.95 a_Y = r.qnorm((1-conf)/2)._sage_() b_Y = -a_Y print "cuantiles estándar", (a_Y, b_Y) 
       
cuantiles estándar (-1.95996398454005, 1.95996398454005)
cuantiles estándar (-1.95996398454005, 1.95996398454005)
a_mu = -b_Y * (sigma/sqrt(n)) + xbarra b_mu = -a_Y * (sigma/sqrt(n)) + xbarra a_mu = a_mu.n() b_mu = b_mu.n() print "intervalo", (a_mu, b_mu) 
       
intervalo (231.863208453703, 235.369298616010)
intervalo (231.863208453703, 235.369298616010)
 
       

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="red",legend_label=u"normal estándar",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..30]] 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)] print d 
       
[233.2610140007237, 232.01677133852925, 234.60504296947695,
233.64185757761024, 234.98925176996133]
[233.2610140007237, 232.01677133852925, 234.60504296947695, 233.64185757761024, 234.98925176996133]
xbarra = mean(d) s2 = variance(d) print "media", xbarra print "varianza", s2 
       
media 233.702787531
varianza 1.37764550673
media 233.702787531
varianza 1.37764550673

Los cuantiles de $t$ dan un intervalo un poco más grande que los de $N(0,1)$. El efecto es mayor cuanto menos datos tengamos.

conf = 0.95 a_Y = r.qt((1-conf)/2, n-1)._sage_() b_Y = -a_Y print "cuantiles t", (a_Y, b_Y) 
       
cuantiles t (-2.77644510519779, 2.77644510519779)
cuantiles t (-2.77644510519779, 2.77644510519779)
a_mu = -b_Y * (sqrt(s2)/sqrt(n)) + xbarra b_mu = -a_Y * (sqrt(s2)/sqrt(n)) + xbarra a_mu = a_mu.n() b_mu = b_mu.n() print "intervalo", (a_mu, b_mu) 
       
intervalo (232.245407451518, 235.160167611003)
intervalo (232.245407451518, 235.160167611003)
 
       

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..10]] animate(frames).show(delay=30) 
       

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 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 
       
[235.74317043081484, 236.1909547593392, 236.18072383406033,
233.27142955148122, 236.78334483771417]
[235.74317043081484, 236.1909547593392, 236.18072383406033, 233.27142955148122, 236.78334483771417]
s2 = variance(d) print "varianza", s2 
       
varianza 1.88093909737
varianza 1.88093909737
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 chi2", (a_Y, b_Y) 
       
cuantiles chi2 (0.48441855708793, 11.1432867818778)
cuantiles chi2 (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.6751828734867765, 15.531519755775543)
intervalo (0.6751828734867765, 15.531519755775543)
 
       

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", n_X print "var_poblacio:", sigma_X^2 print "var_muestral:", s2_X print print "Y: n", n_Y print "var_poblacio:", sigma_Y^2 print "var_muestral:", s2_Y print print "cociente real", sigma_X/sigma_Y.n() 
       
X: n 50
var_poblacio: 4
var_muestral: 3.49492243467

Y: n 10
var_poblacio: 9
var_muestral: 7.62337940747

cociente real 0.666666666666667
X: n 50
var_poblacio: 4
var_muestral: 3.49492243467

Y: n 10
var_poblacio: 9
var_muestral: 7.62337940747

cociente real 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 F", (a_Z, b_Z) 
       
cuantiles F (0.419013382865503, 3.47468082184094)
cuantiles F (0.419013382865503, 3.47468082184094)
a_cociente = s2_X/s2_Y * 1/b_Z b_cociente = s2_X/s2_Y * 1/a_Z print "cociente real", sigma_X/sigma_Y.n() print print "estimación puntual", variance(d_X)/variance(d_Y) print "intervalo", (a_cociente, b_cociente) 
       
cociente real 0.666666666666667

estimación puntual 0.458447920256
intervalo (0.13193957769415343, 1.0941128350614577)
cociente real 0.666666666666667

estimación puntual 0.458447920256
intervalo (0.13193957769415343, 1.0941128350614577)
 
       

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", n_X print "mu:",mu_X print "xb:",xbarra print "var_poblacio:", sigma_X^2 print "var_muestral:", s2_X print print "Y: n", n_Y print "mu:",mu_Y print "yb:",ybarra print "var_poblacio:", sigma_Y^2 print "var_muestral:", s2_Y print print "diferencia real", mu_X-mu_Y 
       
X: n 50
mu: 234
xb: 233.651271719
var_poblacio: 4
var_muestral: 3.6650502945

Y: n 10
mu: 157
yb: 157.475794678
var_poblacio: 9
var_muestral: 5.51991858237

diferencia real 77
X: n 50
mu: 234
xb: 233.651271719
var_poblacio: 4
var_muestral: 3.6650502945

Y: n 10
mu: 157
yb: 157.475794678
var_poblacio: 9
var_muestral: 5.51991858237

diferencia real 77
conf = 0.95 a_Z = r.qnorm((1-conf)/2)._sage_() b_Z = -a_Z print "cuantiles estándar", (a_Z, b_Z) 
       
cuantiles estándar (-1.95996398454005, 1.95996398454005)
cuantiles estándar (-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 "diferencia real", mu_X-mu_Y print print "estimación puntual", xbarra-ybarra print "intervalo", (a_diferencia, b_diferencia) 
       
diferencia real 77

estimación puntual 76.1754770414
intervalo (74.62562647170537, 77.72532761111788)
diferencia real 77

estimación puntual 76.1754770414
intervalo (74.62562647170537, 77.72532761111788)