Código
set.seed(123)
<- 5
theta_real <- matrix(rexp(n = 1000 * 3, rate = theta_real), ncol = 3)
muestra
<- apply(X = muestra, MARGIN = 1, FUN = sum)
suma_muestra
<- 3 / suma_muestra
theta_techo
hist(theta_techo - theta_real, breaks = 100)
Un estimador es una función que calcula una estimación o predicción de un parámetro desconocido en una distribución de probabilidad. En estadística, uno de los criterios importantes para evaluar la calidad de un estimador es si es insesgado.
Un estimador se considera insesgado si su valor esperado es igual al valor verdadero del parámetro que se está estimando. Esto se puede expresar matemáticamente como:
donde
Ejemplo 7.1 Un ejemplo de un estimador insesgado es el promedio muestral
Lo que significa que
Aunque la propiedades ser insesgado pareciera natural en los estimadores, no siempre tenemos esto.
Ejemplo 7.2 Considere una muestra
De acuerdo a los ejemplos de capítulos pasados, el estimador de máxima verosimilitud de
Podemos preguntarnos si R
para estimar el sesgo de
set.seed(123)
<- 5
theta_real <- matrix(rexp(n = 1000 * 3, rate = theta_real), ncol = 3)
muestra
<- apply(X = muestra, MARGIN = 1, FUN = sum)
suma_muestra
<- 3 / suma_muestra
theta_techo
hist(theta_techo - theta_real, breaks = 100)
El histograma muestra la diferencia entre las estimaciones
Como
1 La Gamma Inversa con parámetros
Por lo que
Si por ejemplo
mean(theta_techo - theta_real)
[1] 2.804016
Tomemos otro estimador,
Entonces
Comprobemos que efectivamente
<- 2 / suma_muestra
theta_u mean(theta_u - theta_real)
[1] 0.2026772
Con esto concluimos que el estimador de máxima verosimilitud no siempre es insesgado.
En un mundo ideal, nos gustaría tener estimadores insesgados pero que además tengan varianza pequeña, i.e.,
Ejemplo 7.3 Considere una muestra
Viendo el gráfico queda la pregunta
¿Cómo controlar sesgo y varianza?
Para esto definamos el error cuadrático medio (MSE) de
Escribiendo la defnición de esta cantidad, podemos desagregar el MSE en dos partes:
Si
Ejemplo 7.4 De acuerdo al primer ejemplo del capítulo, nos interesa comparar
Dado que
2 Si
var(theta_u) + mean(theta_u - theta_real)^2
[1] 28.17069
var(theta_techo) + mean(theta_techo - theta_real)^2
[1] 71.15413
Observación. El estimado bayesiano es
<- 4 / (2 + suma_muestra)
theta_bayes var(theta_bayes) + mean(theta_bayes - theta_real)^2
[1] 11.87193
Llegados a este punto, hemos estudiado con mucho detalle los estimadores de la media de una distribución. Hemos visto que en general estos estimadores son insesgados, suficientes, consistentes y eficientes. Sin embargo, en muchas ocasiones nos interesa encontrar un estimador insesgado de la varianza de una distribución. Para esto, consideremos una muestra
Definamos la varianza muestras como
la cual se contrapone con la varianza poblacional
Ambas son idénticas a excepción del denominador
Teorema 7.1 Si
Ejemplo 7.5 Sean
Si
<- matrix(rpois(n = 1000 * 100, lambda = 2), nrow = 100)
muestra
<- apply(muestra, 1, mean)
media <- apply(muestra, 1, var)
varianza <- apply(muestra, 1, function(x, alpha) {
ambos * mean(x) + (1 - alpha) * var(x)
alpha alpha = 0.5)
},
hist(media)
hist(varianza)
hist(ambos)
Ejemplo 7.6 En caso de distribuciones normales, ¿Cuál estimador tiene menor MSE,
Defina
Entonces
Optimizando,
se encuentra que
Ejercicio 7.1 Calcule el MSE de
La Información de Fisher es una herramienta fundamental en inferencia estadística que nos permite cuantificar la cantidad de información que una muestra proporciona acerca de un parámetro desconocido.
Consideremos una variable aleatoria
Ejemplo 7.7 Si
La siguiente función será clave para definir la información de Fisher.
Definición 7.1 Definimos la función Score como:
Definición 7.2 Si
Teorema 7.2 Bajo las condiciones anteriores, y suponiendo que las dos derivadas de
La Información de Fisher es una herramienta esencial en inferencia estadística que nos permite cuantificar la cantidad de información que una muestra proporciona acerca de un parámetro desconocido.
Dada la función score
Ejemplo 7.8 Si
La información de Fisher para esta distribución es:
Ejemplo 7.9 Para
Entonces:
La información de Fisher para esta distribución es:
Estos ejemplos fueron estimados usando solo un dato
Definición 7.3 Suponga que
Observación. La fórmula anterior, no es tan útil como quisieramos. En particular observe que
Ejemplo 7.10 Suponga que una compañía quiere conocer como se comportan sus clientes en sus tiendas. Hay dos propuestas para este modelo
Un modelo Poisson de parámetro
Un modelo donde cada cliente es una v.a. exponencial con tasa de llegada
El tiempo de llegada de cada cliente es independiente.
¿Cuál variable contiene más información de
Solución:
Caso de la variable aleatoria
Acá tenemos que:
Entonces,
Caso de la variable aleatoria
Para
Ambas variables tienen la misma información si
A partir de este ejercicio vamos a hacer un pequeño ejemplo de simulación.
Suponga que
<- 5
theta <- 20 # t = tiempo
tiempo <- tiempo * theta # n = clientes
clientes
<- rpois(n = 1000, lambda = tiempo * theta)
muestra_y <- rgamma(n = 1000, shape = clientes, rate = theta) muestra_x
Según lo estimado ambas informaciones de Fisher debería dar aproximadamente igualdad.
Para
mean(muestra_y / theta^2)
[1] 4.01772
Para
/ theta^2 clientes
[1] 4
Entonces bajo este criterio, ambas variables contienen la misma información, aunque modelen el problema desde ópticas diferentes.
El proceso
hist(muestra_y)
El proceso
hist(muestra_x)
Ejercicio 7.2 Basado en los valores de la simulación, proponga dos valores de
Teorema 7.3 Si
La igualdad se da si y solo si existen funciones
Si
Ejemplo 7.11 Sea
Vea que
y el supuesto 3 se puede verificar por la diferenciabilidad de
Así,
Por ejemplo generemos una secuencia de valores de
<- seq(1, 5, length.out = 100)
beta <- 100
n
<- lapply(
lista_muestras X = beta,
FUN = function(b) {
matrix(rexp(n = n * 500, rate = b), nrow = 500)
}
)
plot(beta, n / beta^2)
Considere el estadístico
La cota de Cramer Rao, si
por lo que
Este comportamiento podemos observarlo con nuestro ejemplo numérico.
<- sapply(
estimador1 X = lista_muestras,
FUN = function(x) {
apply(x, 1, function(xx) (n - 1) / sum(xx))
}
)
plot(beta, apply(X = estimador1, MARGIN = 2, FUN = mean))
plot(beta, apply(X = estimador1, MARGIN = 2, FUN = var))
lines(beta, beta^2 / n, col = "blue")
lines(beta, beta^2 / (n - 2), col = "red")
Ahora, estime
La cota de Cramer es
<- sapply(
estimador2 X = lista_muestras,
FUN = function(x) {
apply(x, 1, function(xx) mean(xx))
}
)
plot(1 / beta, apply(X = estimador2, MARGIN = 2, FUN = mean))
plot(beta, apply(X = estimador2, MARGIN = 2, FUN = var))
lines(beta, 1 / (n * beta^2), col = "blue")
Definición 7.4 sea
Ejemplo 7.12 Sea
Verosimilitud:
Entonces
La cota de Cramer-Rao es
Los otros candidatos para estimar
Teorema 7.4 Teorema. Bajo las condiciones anteriores y si
donde
Teorema 7.5 Recuerde que el MLE
Ejemplo 7.13 Sea
Desde la perspectiva bayesiana, el Estimador de Máxima Verosimilitud
Supongamos que la distribución a priori de