3  Estimación puntual

3.1 La función de distribución empírica y el método de los momentos

Sea la variable aleatoria \(X\) con función de distribución \(F\). Consideramos una muestra aleatoria simple de tamaño \(n\) de \(X\), es decir, \(X_1 ,\ldots, X_n\) v.a.i.i.d. con distribución dada por \(F\) . Sea \(x_1 ,\ldots, x_n\) una realización de esa m.a.s. Se llama función de distribución empírica a la función \[ \begin{align} F_{n}(x)=\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}\mathbf{I}_{(-\infty,x]}(x_{i})¸ \end{align} \] que a cada número real x le asigna la proporción de valores observados que son menores o iguales que x.

Es inmediato comprobar que la función \(F_n\) así definida es una función de distribución:

  1. \(F_n(x) \in [0, 1]\) para todo \(x \in \mathbb{R}\).
  2. \(F_n\) es continua por la derecha.
  3. \(F_n\) es no decreciente.
  4. \(\lim _{x\to -\infty }F_{n}(x)=0\)
  5. \(\lim _{x\to \infty }F_{n}(x)=1\)

Concretamente, \(F_n\) es la función de distribución de una variable aleatoria discreta (que podemos llamar \(X_e\) ) que pone masa \(\frac{1}{n}\) en cada uno de los n puntos \(x_i\) observados:

\(x_i\) \(1\) \(2\) \(\ldots\) \(n\)
\(p_i = P(X_e = x_i)\) \(1/n\) \(1/n\) \(\ldots\) \(1/n\)

A la distribución de \(X_e\) e se le llama distribución empírica asociada al conjunto de valores \({x_1 ,\ldots, x_n}\).

Obsérvese que si fijamos el valor de \(x\) y dejamos variar la muestra, lo que obtenemos es una variable aleatoria. En efecto, se tiene entonces que

\[ \begin{align} F_{n}(x)=\dfrac{1}{n}\displaystyle\sum_{i=1}^{n}\mathbf{I}_{(-\infty,x]}(x_{i})¸ \end{align} \]

donde

\[ \begin{align} \mathbf{I}_{(-\infty,x]}(X_{i})= \left\{ \begin{array}{lcc} 1 & si & X_{i}\leq x \\ \\ 0 & si &X_{i}> x\\ \end{array} \right. \end{align} \] y, por lo tanto, cada término \(\mathbf{I}_{(-\infty,x]}(X_{i})\) es una variable aleatoria de Bernoulli con probabilidad de éxito

\[ \begin{align} p&=P(\mathbf{I}_{(-\infty,x]}(X_{i})=1)\\ &=P(X_{i}\leq x)\\ &=F(x) \end{align} \] De ahí se deduce que \(F_n\) es una variable aleatoria y que \(nF_n(x)\) tiene distribución binomial con parámetros \(n\) y \(p = F(x)\).

3.2 Teorema de Glivenko-Cantelli

El siguiente teorema recoge algunas de las propiedades de la función de distribución empírica.

Theorem 3.1 Sea \(\{X_n\}\) para \(n\geq 1\) , sucesión de variables aleatorias independientes e idénticamente distribuidas definidas en el espacio de probabilidad \((\Omega, \mathcal{A}, P)\) con función de distribución común \(F\) . Se denota por \(F_n\) la función de distribución empírica obtenida de las \(n\) primeras variables aleatorias \(X_1 ,\ldots, X_n\) . Sea \(x\in\mathbb{R}\). Se verifica lo siguiente:

  1. \(P(nF_n(x)=j)=P(F_n(x)=\frac{j}{n})=\binom{n}{j}(F(x))^j(1-F(x))^{n-j}\), \(j=1,\ldots,n\)\
  2. \(E(F_n(x))=F(x)\); \(Var(F_n(x))=\frac{1}{n}F(x)(1-F(x))\).
  3. \(\lim _{n\to \infty }F_{n}(x)=F(x)\)
  4. \(\lim _{n\to \infty }\frac{F_{n}(x)-F(x)}{\sqrt{\frac{F(x)(1-F(x))}{n}}}=Z\), donde \(Z\) es una variable aleatoria con distribución normal estándar y la convergencia es convergencia en distribución.

El siguiente teorema refuerza el resultado (c) anterior, puesto que afirma que la convergencia de \(F_n(x)\) a \(F(x)\) se da uniformemente.

Theorem 3.2 Sea \(\{X_n\}\) para \(n\geq 1\) , sucesión de variables aleatorias independientes e idénticamente distribuidas definidas en el espacio de probabilidad \((\Omega, \mathcal{A}, P)\) con función de distribución común \(F\) . Se denota por \(F_n\) la función de distribución empírica obtenida de las \(n\) primeras variables aleatorias \(X_1 ,\ldots, X_n\) . Sea \(x\in\mathbb{R}\). Entonces \[Sup_{x\in\mathbb{R}} |F_{n}(x)-F(x)|\xrightarrow{c.s}0\]

Note 3.1

Obsérvese que según el apartado (c) del teorema Theorem 3.1, las distribuciones empíricas asociadas a muestras de tamaño n convergen débilmente a la distribución de probabilidad teórica identificada por \(F\), para casi todas las muestras de tamaño infinito que se extraigan de \(F\) . Ésta es una de las consecuencias más importantes del citado teorema: la distribución empírica converge débilmente con probabilidad 1 a la poblacional cuando el tamaño de la muestra tiende a infinito: \[F_{n}(x)\xrightarrow{c.s}F(x)\] Esto garantiza la posibilidad de realizar inferencia estadística:

  1. Los aspectos probabilísticos de una característica \(X\), medida en una población, se resumen de forma estilizada en una distribución de probabilidad \(F\).
  2. La distribución de probabilidad \(F\), puede ser aproximada mediante las distribuciones empíricas \(F_n\) obtenidas por muestreo de la población en estudio.
  3. El teorema de Glivenko-Cantelli afirma que esas aproximaciones son uniformes en x.
  4. Por esta razón el teorema de Glivenko-Cantelli se llama a veces Teorema Fundamental de la Estadística Matemática.

Podemos ver a continuación cómo, a medida que aumentamos el tamaño de la muestra (n=10,30,100,1000), la función de distribución empírica se ajusta cada vez mejor a la distribución teórica normal estándar N(0,1), tal como afirma el teorema.

# Cargar librería para gráficos
library(ggplot2)

# Definir función para graficar ECDF vs distribución teórica
comparar_ecdf_teorica <- function(n, distribucion = "normal") {
  set.seed(123)  # Para reproducibilidad

  # Muestra de tamaño n desde N(0,1)
  muestra <- rnorm(n)
  
  # Dominio común
  x_vals <- seq(-3, 3, length.out = 1000)

  # Distribución teórica
  F_teorica <- pnorm(x_vals)

  # Distribución empírica
  ecdf_muestra <- ecdf(muestra)
  F_empirica <- ecdf_muestra(x_vals)

  # Construir data frame para ggplot
  df <- data.frame(
    x = rep(x_vals, 2),
    F = c(F_empirica, F_teorica),
    Tipo = rep(c("Empírica", "Teórica"), each = length(x_vals))
  )

  # Graficar
  ggplot(df, aes(x = x, y = F, color = Tipo, linetype = Tipo)) +
    geom_line(size = 1) +
    labs(
      title = paste("ECDF vs F(x) — Tamaño de muestra n =", n),
      x = "x", y = "Probabilidad acumulada"
    ) +
    theme_minimal() +
    scale_color_manual(values = c("Empírica" = "red", "Teórica" = "black")) +
    scale_linetype_manual(values = c("Empírica" = "dashed", "Teórica" = "solid"))
}

# Generar gráficos para diferentes tamaños de muestra
comparar_ecdf_teorica(10)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

comparar_ecdf_teorica(30)

comparar_ecdf_teorica(100)

comparar_ecdf_teorica(1000)

El Teorema Fundamental de la Estadística Matemática: da una fundamentación de la inferencia estadística, cuyo objetivo principal consiste en extraer información sobre \(F\) a partir de las observaciones muestrales.

¿Por qué esto es importante?

Porque sin conocer \(F(x)\) explícitamente, podemos estimarla a partir de los datos.
Esto es la base de:

  • los histogramas acumulados,
  • las pruebas no paramétricas,
  • los intervalos de confianza empíricos, y
  • toda la inferencia estadística basada en datos reales.

🎯 Ejemplo: Estimación del percentil 90 del ingreso mensual

📌 Contexto Supón que quieres estimar el ingreso mensual por debajo del cual se encuentran el 90% de las personas en una ciudad.
No conoces la distribución real del ingreso $ F(x) $, pero tienes una muestra de datos.


Paso a paso

  1. Simulamos una población Vamos a suponer que el ingreso sigue una distribución log-normal:
set.seed(123)
poblacion <- rlnorm(1e6, meanlog = 10, sdlog = 0.5)
  1. Tomamos una muestra aleatoria
muestra <- sample(poblacion, size = 20000, replace = FALSE)
  1. Estimamos la distribución empírica
Fn <- ecdf(muestra)
  1. Estimamos el percentil 90 (cuantil 0.9)
cuantil_90_empirico <- quantile(muestra, probs = 0.9)
  1. Comparamos con el valor verdadero en la población
cuantil_90_real <- quantile(poblacion, probs = 0.9)
  1. Resultado
cat("Cuantil 90 estimado (empírico):", round(cuantil_90_empirico, 2), "\n")
Cuantil 90 estimado (empírico): 41926.83 
cat("Cuantil 90 real (poblacional):", round(cuantil_90_real, 2), "\n")
Cuantil 90 real (poblacional): 41805.56 

3.3 La función de distribución empírica y el método de los momentos

3.3.1 Principio de sustitución

En esta sección presentamos una consecuencia importante de la convergencia de \(F_n\) a \(F\) , la definición de estimadores mediante el principio de sustitución.

  1. La convergencia de \(F_n\) a \(F\) permite construir versiones factibles de características poblacionales desconocidas.
  2. Supongamos que estudiamos una característica \(X\) en una población y que el resultado de la observación de \(X\) puede ser modelado como una variable aleatoria con distribución desconocida, digamos \(F\).
  3. Muchas de las preguntas relevantes acerca de la característica \(X\) podrían ser contestadas si su función de distribución \(F\) fuese conocida.
Important 3.1

Preguntas sobre \(X\)

  • el valor esperado,
  • el número de modas de la distribución o
  • la probabilidad de que \(X\) sea negativa

Para fijar ideas podemos pensar que nos interesa conocer cantidades numéricas (parámetros) que dependen únicamente de la función de distribución desconocida \(F\):

\[ \begin{align} \theta=\psi(F) \end{align} \] El teorema de Glivenko-Cantelli nos dice que \(F_n\) se acerca a \(F\), a medida que el tamaño muestral crece. Así, podemos esperar que también se verifique que \[ \begin{align} \hat{\theta}_n=\psi(F_n)\rightarrow\theta=\psi(F) \end{align} \]

Es decir, esperamos que las cantidades numéricas calculadas para la distribución empírica (estimadores) se aproximen a las cantidades desconocidas a medida que el tamaño muestral crezca.

Esta forma de obtener estimadores de parámetros poblacionales desconocidos se denomina principio de sustitución (plug-in principle en inglés). Es un procedimiento muy general de obtención de estimadores.

Sea \(X\sim U(0,\theta)\). Se toma una m.a.s. de \(X\) de tamaño n para estimar \(\theta\). Un estimador razonable de \(\theta\) es el máximo de las observaciones, que es estadístico minimal suficiente para \(\theta\): \[ \begin{align} \hat{\theta}_2 = \max_i {X_i}. \end{align} \] El siguiente código muestra:

  • Para cada tamaño de muestra n, simula valores de \(X_i\)∼U(0,θ),

  • Calcula \[\hat{\theta} = \max_i {X_i}\]

  • Compara con el valor real de \(\theta=10\).

  • Muestra cómo, al aumentar n, el estimador se acerca a \(\theta\).

# Simulación para estimar theta en una uniforme (0, theta)
set.seed(42)
theta_real <- 10

# Tamaños de muestra
n_vals <- c(5, 10, 30, 100)

# Simular y comparar
estimadores <- sapply(n_vals, function(n) {
  muestra <- runif(n, min = 0, max = theta_real)
  max(muestra)
})

# Mostrar resultados
data.frame(
  Tamaño_muestra = n_vals,
  Estimador_maximo = round(estimadores, 3),
  Error = round(theta_real - estimadores, 3)
)
  Tamaño_muestra Estimador_maximo Error
1              5            9.371 0.629
2             10            9.347 0.653
3             30            9.889 0.111
4            100            9.828 0.172

Convergencia del estimador plug-in en la distribución uniforme

En este ejemplo, estimamos el parámetro \(\theta\) de una distribución \(X \sim \mathcal{U}(0, \theta)\) usando el estimador \(\hat{\theta}_n = \max(X_i)\). Este es un estimador tipo plug-in: se usa la distribución empírica para estimar una característica de la distribución teórica.

A continuación, simulamos cómo este estimador converge al verdadero valor \(\theta = 10\) a medida que el tamaño de muestra \(n\) crece.

# Chunk de R — solo código
#| label: fig-convergencia-max
#| fig-cap: 'Convergencia del estimador plug-in $ \hat{\theta}_n = \max X_i $ hacia el valor real $ \theta = 10 $'
#| fig-align: center
#| message: false
#| warning: false

set.seed(42)
theta_real <- 10

# Vector de tamaños de muestra crecientes
n_seq <- seq(5, 500, by = 5)

# Calcular el estimador para cada n
estimadores <- sapply(n_seq, function(n) {
  muestra <- runif(n, min = 0, max = theta_real)
  max(muestra)
})

# Crear data frame para graficar
df_estimacion <- data.frame(
  n = n_seq,
  estimador = estimadores
)

# Graficar
library(ggplot2)
ggplot(df_estimacion, aes(x = n, y = estimador)) +
  geom_line(color = "steelblue", size = 1) +
  geom_hline(yintercept = theta_real, color = "red", linetype = "dashed") +
  labs(
    title = expression("Convergencia de " * hat(theta)[n] * " al valor real " * theta),
    x = "Tamaño de muestra (n)",
    y = expression(hat(theta)[n])
  ) +
  theme_minimal()

3.4 Método de momentos

Una aplicación del principio de sustitución es la definición de los estimadores basados en momentos. El momento no centrado de orden \(k\) de una variable aleatoria \(X\) con distribución \(F\) se define como \[ \begin{align} \mu_k=E_F(X^k)=\int x^kdF(x) \end{align} \] Si \(X_e\) es una variable aleatoria con función de distribución igual a \(F_n\) , la función de distribución empírica de una m.a.s. de tamaño \(n\) de \(X\), se tiene que sus (a los que llamaremos \(m_{k,n}\)) son de la forma \[\begin{align} m_{k,n}=E_{F_n}(X_e^k)=\int x^kdF_n(x)=\frac{1}{n}\sum_{i=1}^{n}X_i^k, \end{align}\] y se denominan momentos muestrales no centrados de orden \(k\). Por ejemplo, \(µ_1\) es la esperanza poblacional y \(m_{1,n}\) la media muestral.

La siguiente proposición garantiza que los momentos muestrales convergen a los poblacionales.

Proposition 3.1 Sea \(X\) variable aleatoria con \(E(X^{2k}) < \infty\). Entonces se verifica que \(m_{k,n} \rightarrow \mu_k\) casi seguro. Además, \[ \begin{align} \frac{\sqrt{n}(m_{k,n}-\mu_k)}{\sqrt{\mu_{2k}-\mu_k^2}}\xrightarrow{d}Z, \end{align} \] con \(Z\sim N(0,1)\).

Proof. Si \(Y_i=X_i^k\) entonces \(m_{k,n}=E_{k,n}(Y_i)=E_{n}(X_i^k)=\frac{1}{n}\sum_{i=1}^{n}X_i^k=\frac{1}{n}\sum_{i=1}^{n}Y_i=\bar{Y}_n\).\

Aplicando la ley fuerte de los grandes números se tiene que \[\begin{align} \lim _{n\to \infty }\frac{S_n-E(S_n)}{n}&=\lim _{n\to \infty }\frac{\sum_{i=1}^{n}X_i^k-E(\sum_{i=1}^{n}X_i^k)}{n}\nonumber\\ &=\lim _{n\to \infty }\left[\frac{\sum_{i=1}^{n}X_i^k}{n}-\frac{E(\sum_{i=1}^{n}X_i^k)}{n}\right]\nonumber\\ &=\lim _{n\to \infty }\left[\bar{Y}_n-\frac{\sum_{i=1}^{n}E(X_i^k)}{n}\right]\nonumber\\ &=\lim _{n\to \infty }\left[\bar{Y}_n-\frac{nE(X^k)}{n}\right]\mbox{Por ser las $X_i$ una mas de X}\nonumber\\ &=\lim _{n\to \infty }\left[\bar{Y}_n-\bar{Y}\right]=0 \mbox{ Aplicando L.F.G.N} \end{align}\]

Por lo anterior, se tiene que \(m_{k,n}=\bar{Y}_n\xrightarrow{c.s} \mu_k=\bar{Y}=E_F(X^k)\). Por otro lado, veamos \[\begin{align} \frac{S_n-E(S_n)}{\sqrt{Var(S_n)}}=\frac{m_{k,n}-E(X^k)}{\sqrt{\frac{Var(X^k)}{n}}}\xrightarrow{d}Z \end{align}\] donde \(S_n=\sum_{i=1}^{n}X_i^k\). Sabemos que \(E(X^k)=\mu_k\) y que \(Var(X^k)=E[(X^k)^2]-[E(X^k)]^2=\mu_{2k}-\mu_k^2\), de ahí se sigue el resultado.

Muchas características poblacionales de interés se pueden expresar como función de los momentos no centrados de órdenes \(1,\ldots, k\): \(\theta=h(\mu_1, \ldots, \mu_k)\). Por ejemplo, la varianza de \(X\) se expresa como \(\sigma^2 = h(\mu_1, \mu_2) = \mu_2-\mu_1^2\).

El estimador de \(\theta\) basado en el principio de sustitución se conoce como estimador de los momentos de \(\theta\) y será \[\begin{align} \hat{\theta}_n = h(m_{1,n} ,\ldots, m_{k,n}). \end{align}\] Obsérvese que el estimador de los momentos de \(\theta\) puede no ser único, porque diferentes funciones \(h\) pueden conducir al mismo valor \(\theta\).

Proposition 3.2 Consideremos la variable aleatoria \(X\) con \(E(X_{2k})< \infty\). Sea \(\theta=h(\mu_{n} ,\ldots, \mu_{n})\). Si \(h\) es continua en \((\mu_{n} ,\ldots, \mu_{n})\), entonces \(\hat{\theta}_n=h(m_{1,n} ,\ldots, m_{k,n})\) converge a \(\theta\) casi seguro. Además, si \(h\) es derivable en \((\mu_{n} ,\ldots, \mu_{n})\), entonces la distribución límite de \(\hat{\theta}_n\) es normal: \[\begin{align} \sqrt{n}(\hat{\theta}_n-\theta)\xrightarrow{d} N(0, \sigma_{h,\theta}^2) \end{align}\]

Example 3.1 Sea \(X\sim U(0,\theta)\). Se toma una m.a.s. de \(X\) de tamaño n para estimar \(\theta\). Un estimador de momentos \(\hat{\theta}_M\) de \(\theta\) viene dado por la siguiente relación \[\begin{align} E(X)=\frac{\theta}{2}\Longrightarrow m_{1,n}=\frac{\hat{\theta}_M}{2}\Longrightarrow 2m_{1,n}=\hat{\theta}_M\Longrightarrow 2\bar{X}=\hat{\theta}_M \end{align}\]

Example 3.2  

Para la variable aleatoria \(X\) con varianza finita, un estimador para \(\theta=Var(X)\) es \[\begin{align} \hat{\theta}=h(m_{1,n}, m_{2,n})&=m_{2,n}-m_{1,n}^2\nonumber\\ &=\frac{1}{n}\sum_{i=1}^{n}x_i^2-\bar{x}^2\nonumber\\ &=\frac{\sum_{i=1}^{n}x_i^2-n\bar{x}^2}{n}\nonumber\\ &=\frac{\sum_{i=1}^{n}(x_i-\bar{x})^2}{n}\nonumber\\ &=\frac{(n-1)S_n^2}{n}\nonumber\\ \end{align}\]

Example 3.3  

Si \(X\sim Exp(\lambda)\) con \(E(X)=\frac{1}{\lambda}\), entonces \(m_{1,n}=\frac{1}{\hat{\lambda}_M}\) \(\Longrightarrow \hat{\lambda}_M=\frac{1}{m_{1,n}}\Longrightarrow \hat{\lambda}_M=\frac{1}{\bar{X}}\).

Si \(X\sim B(n,p)\), con \(E(X)=np\) y \(Var(X)=npq\), entonces \(m_{1,n}=n\hat{p}\) \(\Longrightarrow\) \(\frac{m_{1,n}}{n}=\hat{p}\) \(\Longrightarrow\) \(\frac{\bar{X}}{n}=\hat{p}\) \(\Longrightarrow\) \(\hat{Var(X)}=n\hat{p}(1-\hat{p})\).

3.5 Estimadores de máxima verosimilitud

3.6 Cálculo del estimador máximo verosímil

Sea \(\underset{\sim}{X} =(X_1,\ldots, X_n)\) una muestra aleatoria simple de una variable aleatoria \(X\) con función de densidad (o de masa de probabilidad) \(f(\underset{\sim}{x}|\theta)\), con \(\theta = (\theta_1,\ldots,\theta_k) \in \Theta \subseteq \mathbb{R}^k\) . Sea \(\mathcal{X}\) el espacio muestral, es decir, el conjunto de todos los posibles valores de \(\underset{\sim}{X}\) . Hemos definido la para \(\underset{\sim}{x} =(x_1,\ldots, x_n )\in \mathcal{X}\) como \[ \begin{align} L(\cdot|\underset{\sim}{x}):&\Theta\rightarrow \mathbb{R}^+\nonumber\\ &\theta\rightarrow L(\theta|\underset{\sim}{x})=f(\underset{\sim}{x}|\theta)=\prod_{i=1}^{n}f(x_i|\theta) \end{align} \]

3.6.1 Preguntas sobre \(X\)

Para cada muestra \(\underset{\sim}{x} \in \underset{\sim}{X}\) , el estimador de máxima verosimilitud \(\hat{\theta}\) de \(\theta\) es el valor de \(\Theta\) que hace máxima la verosimilitud \(L(\cdot|\underset{\sim}{x})\): \[ \begin{align} L(\hat{\theta}| \underset{\sim}{x} ) = \max_{\theta \in \Theta} L(\theta|\underset{\sim}{x}). \end{align} \] Intuitivamente \(\hat{\theta}\) es el valor del parámetro que hace más verosímil la muestra observada. Veremos más adelante que los estimadores de máxima verosimilitud son muy buenos estimadores y que en general tienen propiedades de optimalidad. Además, en muchas ocasiones el estimador máximo verosímil es el que el sentido común nos llevaría a proponer.

Example 3.4 Si \(X\sim Exp(\lambda)\Longrightarrow f(x|\lambda)=\lambda\exp\{-\lambda x\}I_{[0,\infty)}(x)\), \(\lambda>0\)

Se toma una muestra de tamaño \(n=1\) y se observa que \(x=3\). Estudiamos la función de verosimilitud \(L(\lambda|3)=\lambda\exp\{-3\lambda\}\) y buscamos su máximo para \(\lambda>0\).

Buscamos los valores de \(\lambda\) que hacen la derivada cero de \(L(\lambda|3)\): \[ \begin{align} \label{maxver_exp} L^{'}(\lambda|3)&=\exp\{-3\lambda\}(1-3\lambda) \end{align} \tag{3.1}\] Igualando a cero la expresión (Equation 3.1) se tiene que \(\lambda=\frac{1}{3}\). Como \(L^{'}(\lambda|3)\geq 0\) y \[ \begin{align} \lim _{\lambda\to 0}L(\lambda|3)=\lim _{\lambda\to \infty}L(\lambda|3)=0 \end{align} \] Se sigue que el punto crítico de \(L(\lambda|3)\) es un máximo. Así, \(\hat{\lambda}=\frac{1}{3}\)

Example 3.5 Suponga que deseamos estimar \(\theta\), la proporción de personas con tuberulosis en una gran población homogénea. Para hacer esto, seleccionamos de manera aleatoria \(n\) personas para hacer pruebas y encontrar \(x\) de estos que tienen la enfermedad.

Ya que la población es grande y homogénea, asumimos que los \(n\) individuos observados son independiente y que cada uno tiene probabilidad \(\theta\) de tener tuberculosis.

Si \(E\) es el evento tiene tuberculosis \[ \begin{align} P(E;\theta)&=P(\mbox{x entre n tienen tuberculosis})\nonumber\\ &=\binom{n}{x}\theta^x(1-\theta)^{n-x} \end{align} \tag{3.2}\] Observe que \(\binom{n}{x}\) es un factor constante no tendrá efecto sobre la maximización de la expresión (Equation 3.2) sobre \(\theta\), ver (kalbfleisch et al. (1985)). La función de verosimilitud de \(\theta\) es definida como sigue: \[ \begin{align} L(\theta)=cP(E;\theta) \end{align} \tag{3.3}\] Acá \(c\) es alguna constante positiva con respecto a \(\theta\); esto es, \(c\) no es función de \(\theta\), sin embargo esta debe ser función de los datos. Escogemos \(c\) para obtener una expresión simple para \(L(\theta)\), y resultados subsecuentes no dependerán de la escogencia específica hecha.

Usualmente \(P(E;\theta)\) y \(L(\theta)\) son productos de términos y será más conveniente trabajar con logaritmos. La función logverosimilitud es el logaritmo de \(L\):

\[ \begin{align} l(\theta)=log L(\theta) \end{align} \tag{3.4}\] Observe que, por (Equation 3.3) \[ \begin{align} l(\theta)=c^{'}+log P(E;\theta) \end{align} \tag{3.5}\] donde \(c^{'}=log c\) no es una función de \(\theta\).

Por la expresión (Equation 3.2) se tiene que

\[ \begin{align} c=\frac{1}{\binom{n}{x}}, \end{align} \] entonces \[ \begin{align} L(\theta)=\theta^x(1-\theta)^{n-x}, \mbox{para } 0\leq\theta\leq1 \end{align} \tag{3.6}\] La función de log verosimilitud es ahora \[ \begin{align} l(\theta)=xlog(\theta)+(n-x)log(1-\theta), \mbox{para } 0\leq\theta\leq1 \end{align} \tag{3.7}\] El estimador de máxima verosimilitud (MDL) \(\hat{\theta}\) es el valor de \(\theta\) el cual maximiza \(l(\theta)\).

Tomando la derivada respecto a \(\theta\) de (Equation 3.4) se tiene que \[ \begin{align} S(\theta)&=\frac{dl(\theta)}{d\theta}\nonumber\\ &=\frac{x}{\theta}-\frac{n-x}{1-\theta} \end{align} \tag{3.8}\] Haciendo cero a \(S(\theta)\) tiene una única solución \(\theta=\frac{x}{n}\), para \(1\leq x\leq n-1\). Bajo estas mismas condiciones tomamos la segunda derivada y la multiplicamos por \(-1\) en la siguiente expresión

\[ \begin{align} \mathscr{I}(\theta)&=-\frac{dS(\theta)}{d\theta}\nonumber\\ &=\frac{x}{\theta^2}+\frac{n-x}{(1-\theta)^2} \end{align} \tag{3.9}\]

Ya que \(\mathscr{I}(\theta)>0\) en \(\theta=\frac{x}{n}\), la función de verosimilitud tiene un máximo relativo en \(\theta=\frac{x}{n}\). Más aún, \(L(\theta)=0\) para \(\theta=0\) y para \(\theta=1\), hemos encontrado con esto un máximo global en \(\hat{\theta}=\frac{x}{n}\).

A las funciones \(S(\theta)\) y \(\mathscr{I}(\theta)\) definidas en (Equation 3.8}) y (Equation 3.9), se les llama función Score y función de información respectivamente.

3.6.2 Principio de invarianza del estimador máximo verosímil

Sea \(X_1,\ldots, X_n\) muestra aleatoria simple de \(X\sim f(x|\theta)\) y sea \(\hat{\theta}\) el estimador máximo verosímil de \(\theta\). Si estamos interesados en estimar una función \(\tau(\theta)\) del parámetro, podemos hacerlo mediante \(\tau(\hat{\theta})\). Éste es el resultado que garantiza el siguiente teorema y se conoce como principio de invariancia.

Theorem 3.3 Si \(\hat{\theta}\) es el estimador de máxima verosimilitud de \(\theta\), entonces para cualquier función \(\tau\) el estimador de máxima verosimilitud de \(\tau(\theta)\) es \(\tau(\hat{\theta})\).

Example 3.6 Sea \(X_1,\ldots,X_n\) una m.a.s de \(X\sim N(\mu,\sigma^2)\). Podemos probar que el estimador de máxima verosimilitud para \(\mu\), \(\hat{\mu}=\bar{X}\). ¿Cuál es el estimador de máxima verosimilitud de \(\theta_1=3\mu\), \(\theta_2=\mu^2\) y \(\theta_3=1/\mu\)?

Por el principio de invarianza tenemos que \(\hat{\theta_1}=3\bar{X}\), \(\hat{\theta_2}=\bar{X}^2\) y \(\hat{\theta_3}=\frac{1}{\bar{X}}\).

3.7 Error cuadrático medio

Una vez se han presentado diferentes métodos de estimación surge la necesidad de desarrollar criterios para evaluarlos y compararlos de acuerdo a estos criterios. En este tema estudiaremos medidas de la calidad de un estimador. Lo haremos primero para muestras finitas para pasar después a proponer medidas asintóticas de calidad.

Se define el error cuadrático medio (ECM) de un estimador \(W\) de un parámetro \(\theta\) como \[ \begin{align} E_\theta((W-\theta)^2) \end{align} \] Ésta es una medida intuitiva del comportamiento de un estimador: cuanto menor sea el error cuadrático medio mejor será el estadístico \(W\). De hecho, para cualquier función \(\phi\) creciente con \(\phi(0) = 0\), \(E_\theta(\phi(|W- \theta|))\) es una medida razonable de lo alejadas que estarán, en promedio, las estimaciones de \(\theta\) que proporcione W.

En general, se prefiere el error cuadrático medio a otras medidas por ser más tratable analíticamente. Además el error cuadrático medio puede descomponerse \[\begin{align} E_\theta((W-\theta)^2)&=E_\theta((W-E_\theta(W))^2)+E_\theta((E_\theta(W)-\theta)^2)\nonumber\\ &=Var_\theta(W)+(B_\theta(W))^2 \end{align}\] El término \(B_\theta(W) = E_\theta(W)-\theta\) se llama (en inglés bias) de \(W\) cuando se estima \(\theta\) y es una medida de la desviación sistemática que se tiene cuando se estima \(\theta\) por \(W\). Si un estimador tiene sesgo nulo para cualquier valor del parámetro se dice que es un estimador insesgado. En tal caso, \(E_\theta((W-\theta)^2)=Var_\theta(W)\).

3.7.1 Observaciones sobre el ECM

  1. Así, el error cuadrático medio de un estimador es la suma de su varianza (una medida de su dispersión) más el cuadrado de su sesgo (medida de la desviación sistemática o de la exactitud del estimador).

  2. Es una medida conjunta de precisión y exactitud del estimador.

  3. Por lo tanto, parece sensato buscar estimadores que tengan error cuadrático medio pequeño, porque de esta manera controlaremos tanto la dispersión como la exactitud de las estimaciones.

Figure 3.1: Comparación entre precisión y exactitud
(a) Exactitud vs precisión
# install.packages("ggplot2")  # si hace falta
library(ggplot2)

theta    <- 0      # valor "verdadero"
thetahat <- 1.2    # estimación

x  <- seq(theta - 3, thetahat + 3, length.out = 500)
df <- data.frame(x = x,
                 y = dnorm(x, mean = thetahat, sd = 0.9))

ggplot(df, aes(x, y)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0, linewidth = 0.8) +
  # marquitas en el eje para theta y theta-hat
  geom_segment(aes(x = theta,    xend = theta,    y = 0, yend = max(y)*0.03),
               inherit.aes = FALSE) +
  geom_segment(aes(x = thetahat, xend = thetahat, y = 0, yend = max(y)*0.03),
               inherit.aes = FALSE) +
  annotate("text", x = theta,    y = -max(df$y)*0.05, label = expression(theta)) +
  annotate("text", x = thetahat, y = -max(df$y)*0.05, label = expression(hat(theta))) +
  labs(title = "Distribución de estimaciones", x = NULL, y = NULL) +
  coord_cartesian(ylim = c(-max(df$y)*0.08, max(df$y)*1.05), clip = "off") +
  theme_classic(base_size = 14) +
  theme(
    axis.text.y  = element_blank(),
    axis.ticks.y = element_blank(),
    plot.margin  = margin(10, 20, 25, 20) # deja espacio para las etiquetas bajo el eje
  )
Warning in geom_segment(aes(x = theta, xend = theta, y = 0, yend = max(y) * : All aesthetics have length 1, but the data has 500 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = thetahat, xend = thetahat, y = 0, yend = max(y) * : All aesthetics have length 1, but the data has 500 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression

Estimador de \(\theta\) estará indicado por \(\hat{\theta}\).

Quisieramos que \(E(\hat{\theta})=\theta\).

La distribución muestral para un estimador puntual sesgado positivamente, para el que \(E(\hat{\theta})>\theta\), se muestra en la siguiente figura

library(ggplot2)

# --- Parámetros editables ---
theta       <- 0
E_thetahat  <- 0.5      # sesgo: centro de la densidad
thetahat    <- 1.2
sd0         <- 0.7      # dispersión
# ----------------------------

# Curva de densidad (normal centrada en E(\hat{theta}))
x_min <- theta - 1.2
x_max <- thetahat + 1.2
x  <- seq(x_min, x_max, length.out = 800)
y  <- dnorm(x, mean = E_thetahat, sd = sd0)
df <- data.frame(x, y)
ymax <- max(y)

# Alturas EXACTAS de las barras (tocar la curva)
y_theta     <- dnorm(theta,      mean = E_thetahat, sd = sd0)  # altura en theta
y_Ethetahat <- dnorm(E_thetahat, mean = E_thetahat, sd = sd0)  # pico

ggplot(df, aes(x, y)) +
  geom_line(linewidth = 1.1, lineend = "round") +
  # ejes "a mano"
  annotate("segment", x = x_min, xend = x_max, y = 0, yend = 0, linewidth = 0.8) +
  annotate("segment", x = x_min, xend = x_min, y = 0, yend = 1.05*ymax, linewidth = 0.8) +
  # barra corta en theta (hasta la curva)
  geom_segment(aes(x = theta, xend = theta, y = 0, yend = y_theta), inherit.aes = FALSE) +
  # barra larga en E(thetahat) (hasta el pico)
  geom_segment(aes(x = E_thetahat, xend = E_thetahat, y = 0, yend = y_Ethetahat),
               inherit.aes = FALSE) +
  # etiquetas bajo el eje x
  annotate("text", x = theta,      y = -0.06*ymax, label = expression(theta)) +
  annotate("text", x = E_thetahat, y = -0.06*ymax, label = expression(E(hat(theta)))) +
  annotate("text", x = thetahat,   y = -0.06*ymax, label = expression(hat(theta))) +
  # etiqueta del eje y
  annotate("text",
           x = x_min - 0.03*(x_max - x_min), y = 0.9*ymax,
           label = expression(f(hat(theta))), angle = 90) +
  labs(x = NULL, y = NULL) +
  coord_cartesian(xlim = c(x_min, x_max),
                  ylim = c(-0.10*ymax, 1.08*ymax), clip = "off") +
  theme_void(base_size = 14) +
  theme(plot.margin = margin(10, 20, 35, 45))
Warning in geom_segment(aes(x = theta, xend = theta, y = 0, yend = y_theta), : All aesthetics have length 1, but the data has 800 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = E_thetahat, xend = E_thetahat, y = 0, yend = y_Ethetahat), : All aesthetics have length 1, but the data has 800 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression

library(ggplot2)
library(patchwork)

# ===== Panel (a): arco =====
# soporte con ~5% a la izquierda y 95% a la derecha
xL <- -0.10
xR <-  2.10
theta     <- (xL + xR)/2    # centro del arco: ahí debe ir la barra
thetahat1 <- xR

# definir el arco (semicircunferencia escalada)
c_a  <- (xL + xR)/2
R_a  <- (xR - xL)/2
x_a  <- seq(xL, xR, length.out = 700)
hsca <- 0.95
y_a  <- hsca * sqrt(pmax(0, R_a^2 - (x_a - c_a)^2))
dfa  <- data.frame(x = x_a, y = y_a)
ymax_a <- max(dfa$y)

# altura de la curva en theta
y_theta <- approx(x_a, y_a, xout = theta)$y

p_a <-
  ggplot(dfa, aes(x, y)) +
  geom_line(linewidth = 1.2, lineend = "round") +
  # ejes
  annotate("segment", x = xL-0.3, xend = xR+0.3, y = 0, yend = 0, linewidth = 0.8) +
  annotate("segment", x = 0, xend = 0, y = 0, yend = 1.05*ymax_a, linewidth = 0.8) +
  # barra en el centro (theta)
  geom_segment(aes(x = theta, xend = theta, y = 0, yend = y_theta), inherit.aes = FALSE) +
  # etiquetas
  annotate("text", x = theta,     y = -0.08*ymax_a, label = expression(theta)) +
  annotate("text", x = thetahat1, y = -0.08*ymax_a, label = expression(hat(theta)[1])) +
  annotate("text", x = -0.15, y = 0.9*ymax_a,
           label = expression(f(hat(theta)[1])), angle = 90) +
  annotate("text", x = (xL + xR)/2, y = -0.22*ymax_a, label = "(a)") +
  coord_cartesian(xlim = c(xL-0.2, xR+0.2),
                  ylim = c(-0.25*ymax_a, 1.1*ymax_a), clip = "off") +
  theme_void(base_size = 14) +
  theme(plot.margin = margin(10, 20, 40, 45))

# ===== Panel (b): normal (sin cambios) =====
theta_b    <- 1.00
thetahat2  <- 1.60
sd_b       <- 0.25
x_b <- seq(0, theta_b + 3.2*sd_b, length.out = 700)
y_b <- dnorm(x_b, mean = theta_b, sd = sd_b)
dfb <- data.frame(x = x_b, y = y_b)
ymax_b <- max(dfb$y)

p_b <-
  ggplot(dfb, aes(x, y)) +
  geom_line(linewidth = 1.2, lineend = "round") +
  annotate("segment", x = 0, xend = max(x_b)+0.3, y = 0, yend = 0, linewidth = 0.8) +
  annotate("segment", x = 0, xend = 0, y = 0, yend = 1.05*ymax_b, linewidth = 0.8) +
  geom_segment(aes(x = theta_b, xend = theta_b, y = 0, yend = dnorm(theta_b, theta_b, sd_b)),
               inherit.aes = FALSE) +
  annotate("text", x = theta_b,   y = -0.08*ymax_b, label = expression(theta)) +
  annotate("text", x = thetahat2, y = -0.08*ymax_b, label = expression(hat(theta)[2])) +
  annotate("text", x = -0.15, y = 0.9*ymax_b,
           label = expression(f(hat(theta)[2])), angle = 90) +
  annotate("text", x = (min(x_b)+max(x_b))/2, y = -0.22*ymax_b, label = "(b)") +
  coord_cartesian(xlim = c(0, max(x_b)+0.2),
                  ylim = c(-0.25*ymax_b, 1.1*ymax_b), clip = "off") +
  theme_void(base_size = 14) +
  theme(plot.margin = margin(10, 20, 40, 45))

# ===== Combinar =====
p_a | p_b
Warning in geom_segment(aes(x = theta, xend = theta, y = 0, yend = y_theta), : All aesthetics have length 1, but the data has 700 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in geom_segment(aes(x = theta_b, xend = theta_b, y = 0, yend = dnorm(theta_b, : All aesthetics have length 1, but the data has 700 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression
Warning in is.na(x): is.na() aplicado a un objeto que no es (lista o vector) de
tipo 'expression

La figura (b) es la distribución deseada porque una varianza pequeña garantiza que, en un muestreo repetido, una fracción más alta de valores \(\hat{\theta}_2\) estará ``cerca’’ de \(\theta\).

Por consiguiente, además de preferir un estimador insesgado, necesitamos que la varianza de la distribución del estimador \(V(\hat{\theta})\) sea lo más pequeña posible. Dados dos estimadores insesgados de un parámetro \(\theta\) seleccionamos el estimador con la menor varianza, mientras todos los demás parece igual.

Example 3.7 Suponga que \(Y_1,Y_2,Y_3\) denotan una muestra aleatoria ind de una distribuci'on exponencial con funci'on de densidad \[f(y)=\left\{\begin{array}{ll}\left(\frac{1}{\theta}\right)e^{-y/\theta},& y>0,\\0,&\mbox{ en cualquier otro punto}\end{array}\right.\] Considere los siguientes estimadores de \(\theta\): \[\hat{\theta}_1=Y_1,\,\,\,\,\,\,\,\,\,\,\hat{\theta}_2=\frac{Y_1+Y_2}{2},\,\,\,\,\,\,\,\,\,\,\hat{\theta}_3=\frac{Y_1+2Y_2}{3},\,\,\,\,\,\,\,\,\,\,\]

\[\hat{\theta}_4=\min(Y_1,\,Y_2,\,Y_3),\,\,\,\,\,\,\,\,\,\,\hat{\theta}_5=\overline{Y}\]

  1. ¿Cuáles son insesgados?

Sol: Determinemos si \(\hat{\theta}_4=\min(Y_1,\,Y_2,\,Y_3)\) es insesgado. Para esto, determinemos la distribución de \(\hat{\theta}_4\). Supongamos que \(Y_i\sim Exp(\lambda_i)\) con \(i=1,2,3\). \[ \begin{align} P(\hat{\theta}_4>a)&=P(\min(Y_1,\,Y_2,\,Y_3)>a)\nonumber\\ &=P(Y_1>a)P(Y_2>a)P(Y_3>a)\nonumber\\ &=\exp\{-a\lambda_1\}\exp\{-a\lambda_2\}\exp\{-a\lambda_3\}\nonumber\\ &\exp\left\{-a\left(\sum_{i=1}^{3}\lambda_i\right)\right\} \end{align} \] Así que \[ \begin{align} F_{\hat{\theta}_4}(a)&=1-P(\hat{\theta}_4>a)\nonumber\\ &=1-\exp\left\{-a\left(\sum_{i=1}^{3}\lambda_i\right)\right\} \end{align} \] De ahí que \(\hat{\theta}_4\sim Exp\left(\sum_{i=1}^{3}\lambda_i\right)\) entonces, \(E(\hat{\theta}_4)=\frac{1}{\sum_{i=1}^{3}\lambda_i}\). Pero por hipótesis \(\lambda_i=\frac{1}{\theta}\), para todo \(i=1,2,3\). Luego,

\[ \begin{align} E(\hat{\theta}_4)&=\frac{1}{\sum_{i=1}^{3}\lambda_i}\nonumber\\ &=\frac{1}{3\lambda} \mbox{ por ser las $Y_i$ son iid}\nonumber\\ &=\frac{1}{\frac{3}{\theta}}\nonumber\\ &=\frac{\theta}{3} \end{align} \] Luego \(\hat{\theta}_4\) no es insesgado, es sesgado.

3.8 Eficiencia relativa

Definition 3.1 Un estimador \(W\) de \(\theta\) se denomina inadmisible si existe otro estimador \(V\) de \(\theta\) tal que \[ \begin{align} E_\theta((V-\theta)^2)\leq E_\theta((W-\theta)^2) \mbox{ para todo $\theta \in \Theta$} \end{align} \]

3.9 Mejor estimador insesgado.

Definition 3.2 Sean \(\hat{\theta}_1\) y \(\hat{\theta}_2\) dos estimadores insesgados de un parámetro \(\theta\), con varianzas \(V (\hat{\theta}_1)\) y \(V (\hat{\theta}_2)\), respectivamente. Entonces la eficiencia de \(\hat{\theta}_1\) con respecto a \(\hat{\theta}_2\) denotada como \(eff(\hat{\theta}_1,\,\hat{\theta}_2)\) se defiene por la expresión

\[eff (\hat{\theta}_1,\, \hat{\theta}_2)=\frac{V(\hat{\theta}_1)}{V(\hat{\theta}_2)}.\]

Example 3.8 Diremos que \(\hat{\theta}_1\) es más eficiente que \(\hat{\theta}_2\) si \[eff(\hat{\theta}_1,\,\hat{\theta}_2)<1\mbox{ si y sólo si }V(\hat{\theta}_2)>V(\hat{\theta}_1)\] luego \(\hat{\theta}_1\) es un mejor estimador insesgado que \(\hat{\theta}_2\).

Example 3.9 Si \(eff(\hat{\theta}_1,\,\hat{\theta}_2)=1.8\) entonces \(V(\hat{\theta}_1)=1.8V(\hat{\theta}_2)\) entonces \(\hat{\theta}_2\) se prefiere a \(\hat{\theta}_1\).

Si \(eff(\hat{\theta}_1,\,\hat{\theta}_2)=0.73\) entonces \(V(\hat{\theta}_1)=0.73V(\hat{\theta}_2)\) entonces \(\hat{\theta}_1\) se prefiere a \(\hat{\theta}_2\).

Theorem 3.4 Teorema de Lehmann-Scheffé Si \(\hat{\theta}\) es un estimador insesgado para \(\theta\) y si \(U\) es un estadístico suficiente para \(\theta\), entonces hay una función de \(U\) que también es un estimador insesgado para \(\theta\) y tiene una varianza no mayor que \(\hat{\theta}\).

En símbolos:

Si \(E(\hat{\theta})=\theta\) y \(U\) es suficiente para \(\theta\) entonces existe una \(f(U)\) tal que \[E(f(U))=\theta\] y \[Var(\hat{\theta})\geq Var(f(U))\]

Example 3.10 Contexto del teorema

Un estimador insesgado de \(\theta\) es una variable aleatoria \(\hat{\theta}\) que cumple:

\[ E(\hat{\theta}) = \theta. \]

Un estadístico suficiente \(U\) para \(\theta\) contiene toda la información de la muestra acerca de \(\theta\). Es decir, dados \(U\), la muestra completa no aporta nada más sobre \(\theta\).

El Teorema de Lehmann–Scheffé dice que si combinas estos dos elementos:

  • Un estimador insesgado cualquiera \(\hat{\theta}\).
  • Un estadístico suficiente \(U\).

Entonces puedes construir una función de \(U\), digamos \(f(U)\), que:

  • también es insesgado,
  • tiene varianza menor o igual a la de \(\hat{\theta}\).

De hecho, el resultado formal se apoya en el Teorema de Rao–Blackwell:

\[ f(U) = E(\hat{\theta} \mid U). \]

Esta transformación se llama Rao-Blackwellización.
Y el teorema de Lehmann–Scheffé añade que, si \(U\) además de suficiente es completo, entonces ese \(f(U)\) es el único estimador insesgado de varianza mínima (UMVUE: Uniformly Minimum Variance Unbiased Estimator).


2. Intuición

  • Empiezas con cualquier estimador insesgado \(\hat{\theta}\).
  • Si “condicionas” en la estadística suficiente \(U\), eliminas ruido innecesario que no aporta información sobre \(\theta\).
  • Eso reduce (o al menos no aumenta) la varianza del estimador.
  • El resultado es un estimador más eficiente.

3. Ejemplo clásico

Contexto Supongamos que:

\[ X_1, X_2, \dots, X_n \overset{iid}{\sim} \text{Bernoulli}(p), \]

donde el parámetro de interés es \(p\).


Paso 1: un estimador insesgado cualquiera

Consideremos el estimador basado sólo en el primer dato:

\[ \hat{p} = X_1. \]

Claramente:

\[ E(\hat{p}) = E(X_1) = p. \]

Así que es insesgado, pero muy ineficiente (usa un solo dato).


Paso 2: estadístico suficiente

El número total de éxitos:

\[ U = \sum_{i=1}^n X_i \]

es un estadístico suficiente para \(p\) (por factorización de la verosimilitud).


Paso 3: aplicar Rao–Blackwell

Construimos:

\[ f(U) = E(\hat{p} \mid U). \]

Como \(\hat{p} = X_1\), necesitamos \(E(X_1 \mid U)\).
Por simetría, dado que hay \(U\) éxitos en total entre \(n\) ensayos, cada \(X_i\) tiene la misma probabilidad de ser 1. Entonces:

\[ E(X_1 \mid U) = \frac{U}{n}. \]

Así:

\[ f(U) = \frac{U}{n}. \]


Paso 4: verificar propiedades

  • Insesgado:

\[ E\left(\frac{U}{n}\right) = \frac{1}{n}E(U) = \frac{1}{n}(np) = p. \]

  • Varianza:

\[ Var(X_1) = p(1-p), \]

\[ Var\left(\frac{U}{n}\right) = \frac{1}{n^2} Var(U) = \frac{1}{n^2}(np(1-p)) = \frac{p(1-p)}{n}. \]

Y efectivamente:

\[ \frac{p(1-p)}{n} \leq p(1-p). \]


4. Conclusión

El estimador inicial \(\hat{p} = X_1\) era insesgado pero ineficiente.

Rao–Blackwellizando con el suficiente \(U\), obtuvimos:

\[ f(U) = \frac{U}{n}, \]

que es la media muestral.

Por Lehmann–Scheffé, como \(U\) es suficiente y completo, \(\bar{X} = U/n\) es el UMVUE de \(p\).

3.10 Teorema de Cramér-Rao. Información de Fisher

Definition 3.3 Diremos que un estimador \(W^{*}\) es el mejor estimador insesgado de \(\tau(\theta)\), o el UMVUE (Uniformly Minimum Variance Unbiased Estimator) (estimador insesgado de \(\tau(\theta)\) uniformemente de mínima varianza) , si \(E_\theta(W^{*}) = \tau(\theta)\) para todo \(\theta\in\Theta\) y si para cualquier otro estimador \(W\) , tal que \(E_\theta(W) = \tau(\theta)\) para todo \(\theta\in\Theta\), se tiene que \(V_\theta(W^{*}) \leq V_\theta(W)\), para todo \(\theta\in\Theta\).

La búsqueda del UMVUE no debe consistir en repasar todos los estimadores insesgados posibles. El siguiente resultado aborda el problema de un modo diferente: establece una cota inferior para la varianza de todos los estimadores insesgados de un parámetro. Así, si encontramos un estimador insesgado cuya varianza iguale esa cota podremos concluir que ese estimador es el UMVUE.

Theorem 3.5 Sea \(\underset{\sim}{X} =(X_1,\ldots, X_n)\) una variable aleatoria n-dimensional con función de densidad conjunta \(f ( \underset{\sim}{x}|\theta)\), \(\theta\in\Theta\subseteq\mathbb{R}\). Sea \(W(\underset{\sim}{X} )\) un estimador insesgado para \(\tau(\theta)\), es decir, \(E_\theta(W(\underset{\sim}{X})) = \tau(\theta)\) para todo \(\theta\), donde \(\tau\) es una función de \(\theta\) que cumple:

H1: \(\tau(\theta)\) es diferenciable en \(\theta\).

H2: Se supone además que la verosimilitud conjunta \(f(\underset{\sim}{x}|\theta)\) verifica que para cualquier función \(h(\underset{\sim}{x})\) tal que \(E_\theta|h(\underset{\sim}{X})| < \infty\) se tiene que

\[ \begin{align} &\frac{d}{d\theta}\int\dotsc\int h(\underset{\sim}{x})f (\underset{\sim}{x}|\theta)dx_1\ldots dx_n\nonumber\\ &=\int\dotsc\int h(\underset{\sim}{x})\left[\frac{\partial}{\partial\theta}f (\underset{\sim}{x}|\theta)\right]dx_1\ldots dx_n \end{align} \]

Entonces, \[\begin{align} V_\theta(W(\underset{\sim}{X} ))\geq\frac{\left(\frac{d}{d\theta}\tau(\theta)\right)^2}{E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f(\underset{\sim}{X}|\theta)\right)^2\right]} \end{align}\] A la cantidad del lado derecho de la desigualdad anterior se la denomina Cota de Cramér-Rao.

Note 3.2

Nota: El teorema de Cramér-Rao es igualmente válido en el caso discreto. En este caso la hipótesis H2 afirma que pueden intercambiarse el sumatorio y la diferenciación.

  1. Un estimador insesgado para \(\tau(\theta)\) se denomina si su varianza es la mínima posible, es decir, si es igual a la cota de Cramér-Rao.

  2. La eficiencia de un estimador insesgado se define como el cociente entre la cota de Cramér-Rao y su varianza.

  3. Es un valor menor o igual que 1 si se dan las hipótesis del teorema de Cramér-Rao.

En la demostración del teorema de Cramér-Rao se ha probado que

\[ \begin{align} E_\theta(S(\theta))&=E_\theta(S(\theta|\underset{\sim}{X}))\nonumber\\ &=E_\theta\left(\frac{\partial}{\partial \theta}\log L(\theta|\underset{\sim}{X})\right)=0 \end{align} \] Obsérvese que para obtener el estimador máximo verosímil de \(\theta\) lo que se hace es resolver la ecuación \[\begin{align} S(\theta|\underset{\sim}{X})=0 \end{align}\] lo que equivale a buscar el valor de \(\theta\) para el cual el valor de \(S(\theta|\underset{\sim}{X})\) coincide con su valor esperado.

A la cantidad que aparece en el denominador de la cota de Cramér-Rao se le denomina cantidad de información de Fisher que sobre \(\theta\) contiene el vector \(\underset{\sim}{X}\):

\[ \begin{align} \mathscr{I}_{E}(\theta)=\mathscr{I}_{\underset{\sim}{X}}(\theta)&=E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{\underset{\sim}{X}}(\underset{\sim}{X}|\theta)\right)^2\right]\\ &=V\left(\frac{\partial}{\partial \theta}\log f_{\underset{\sim}{X}}(\underset{\sim}{X}|\theta)\right)\\ &=V(S(\theta|\underset{\sim}{X})) \end{align} \]

Se denomina cantidad de información de Fisher que sobre \(\theta\) contiene la variable \(X_i\) a \[ \begin{align} &=E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X_i}(X|\theta)\right)^2\right]\\ &=V\left(\frac{\partial}{\partial \theta}\log f_{X_i}(X|\theta)\right)\\ &=V(S(\theta|X_i)) \end{align} \]

Cuando \(\underset{\sim}{X}=(X_1,\ldots, X_n)\) es una muestra aleatoria simple de \(X\) se verifica que la información de Fisher contenida en la muestra es la suma de las informaciones contenidas en cada una de las observaciones y, dado que éstas son idénticamente distribuidas, se tiene que

\[ \begin{align} \mathscr{I}_{E}(\theta)=\mathscr{I}_{\underset{\sim}{X}}(\theta)=n\mathscr{I}_{X}(\theta) \end{align} \] Este resultado es consecuencia del siguiente corolario del teorema de Cramér-Rao:

Corollary 3.1 Bajo las hipótesis del teorema de Cramér-Rao, si \(\underset{\sim}{X}=(X_1,\ldots, X_n)\) es una muestra aleatoria simple de \(X\) con distribución dada por \(f(x|\theta)\) entonces \[ \begin{align*} \mathscr{I}_{E}(\theta)=\mathscr{I}_{\underset{\sim}{X}}(\theta)&=E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{\underset{\sim}{X}}(\underset{\sim}{X}|\theta)\right)^2\right]\\ &=nE_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X}(X|\theta)\right)^2\right]\\ &=n\mathscr{I}_{X}(\theta) \end{align*} \]

Proof. Por independencia, la verosimilitud de \(\underset{\sim}{X}\) es el producto de verosimilitudes, luego \[\begin{align} \frac{\partial}{\partial \theta}\log f_{\underset{\sim}{X}}(\underset{\sim}{x}|\theta)&=\frac{\partial}{\partial \theta}\sum_{i=1}^{n}\log f_{X}(x_i|\theta)\nonumber\\ &=\sum_{i=1}^{n}\frac{\partial}{\partial \theta}\log f_{X}(x_i|\theta) \end{align}\] Por lo tanto,

\[ \begin{align} E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{\underset{\sim}{X}}(\underset{\sim}{x}|\theta)\right)^2\right]&=E_\theta\left[\left(\sum_{i=1}^{n}\frac{\partial}{\partial \theta}\log f_{X_i}(x_i|\theta)\right)^2\right]\nonumber\\ =&\sum_{i=1}^{n}E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X_i}(x_i|\theta)\right)^2\right]\nonumber\\ +&\sum_{i\neq j}E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X_i}(x_i|\theta)\frac{\partial}{\partial \theta}\log f_{X_j}(x_j|\theta)\right)\right] \end{align} \tag{3.10}\]

La segunga igualdad de Equation 3.10 se tiene por que \(\left(\sum_{i=1}^{n}a_i\right)^2=\sum_{i=1}^{n}a_i^2+\sum_{i\neq j}a_ia_j\). Pero la segunda sumatoria es igual a cero debido a la independencia entre \(X_i\) y \(X_j\) y dado que las funciones score tienen esperanza 0, según se vio en la demostración del teorema de Cramér-Rao. \[ \begin{align} E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{\underset{\sim}{X}}(\underset{\sim}{x}|\theta)\right)^2\right] =&\sum_{i=1}^{n}E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X_i}(x_i|\theta)\right)^2\right]\nonumber\\ +&\sum_{i\neq j}E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X_i}(x_i|\theta)\right)\left(\frac{\partial}{\partial \theta}\log f_{X_j}(x_j|\theta)\right)\right]\nonumber\\ =&\sum_{i=1}^{n}E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X_i}(x_i|\theta)\right)^2\right]\nonumber\\ \end{align} \tag{3.11}\]

Lemma 3.1 Si la función de verosimilitud satisface:

H3: Se supone que la verosimilitud conjunta \(f (\underset{\sim}{x}|\theta)\) verifica que para cualquier función \(h(\underset{\sim}{x})\) tal que \(E_\theta|h(\underset{\sim}{X})| < \infty\) se tiene que \[ \begin{align} &\frac{\partial^2}{\partial\theta^2}\int\dotsc\int h(\underset{\sim}{x})f (\underset{\sim}{x}|\theta)dx_1\ldots dx_n\nonumber\\ &=\int\dotsc\int h(\underset{\sim}{x})\left[\frac{\partial^2}{\partial\theta^2}f (\underset{\sim}{x}|\theta)\right]dx_1\ldots dx_n \end{align} \]

Entonces \[ \begin{align} \mathscr{I}_{X}(\theta)=E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X}(x|\theta)\right)^2\right]=-E_\theta\left[\frac{\partial^2}{\partial \theta^2}\log f_{X}(x|\theta)\right] \end{align} \]

Proof. \[ \begin{align} \frac{\partial^2}{\partial \theta^2}\log f_{X}(x|\theta)&=\frac{\partial}{\partial \theta}\left[\frac{1}{f_{X}(x|\theta)}\left(\frac{\partial}{\partial \theta}f_{X}(x|\theta)\right)\right]\nonumber\\ &=-\frac{1}{f^2_{X}(x|\theta)}\left(\frac{\partial}{\partial \theta}f_{X}(x|\theta)\right)^2\nonumber\\ &+\frac{1}{f_{X}(x|\theta)}\frac{\partial^2}{\partial \theta^2}f_{X}(x|\theta) \end{align} \] Por otro lado, \[ \begin{align} E_\theta\left[\frac{1}{f_{X}(x|\theta)}\frac{\partial^2}{\partial \theta^2}f_{X}(x|\theta)\right]&=\int\frac{\partial^2}{\partial \theta^2}f_{X}(x|\theta)dx\nonumber\\ &\overset{\text{H3}}{=}\frac{d^2}{d \theta^2}\int f_{X}(x|\theta)dx\nonumber\\ &=\frac{d^2}{d \theta^2}1 =0 \end{align} \] Así pues, \[ \begin{align} E_\theta\left[\frac{\partial^2}{\partial \theta^2}\log f_{X}(x|\theta)\right] &=-E_\theta\left[\frac{1}{f^2_{X}(x|\theta)}\left(\frac{\partial}{\partial \theta}f_{X}(x|\theta)\right)^2\right]\nonumber\\ &=-E_\theta\left[\left(\frac{\partial}{\partial \theta}\log f_{X}(x|\theta)\right)^2\right]\nonumber\\ &=-\mathscr{I}_{X}(\theta) \end{align} \]

  1. En general, el teorema de Cramér-Rao no es aplicable si el soporte de \(f(x|\theta)\) depende del parámetro \(\theta\) debido a que la derivada y la integral no son intercambiables si los límites de integración dependen de \(\theta\).

  2. Aunque el teorema de Cramér-Rao pueda ser aplicado y la cota de Cramér-Rao sea efectiva, no hay garantías de que esta cota sea alcanzada por algún estimador insesgado del parámetro.

Corollary 3.2 Sea \(X_1 , \ldots, X_n\) una muestra aleatoria simple de \(X\) con distribución dada por \(f(x|\theta)\), \(\theta\in\mathbb{R}\), donde \(f\) satisface las hipótesis del teorema de Cramér-Rao. Sea \(L(\theta|\underset{\sim}{x}) = \prod_{i=1}^{n}f(x_i|\theta)\) la función de verosimilitud. Sea \(W (\underset{\sim}{x}) = W(X_1 , \ldots, X_n)\) un estimador insesgado de \(\tau(\theta)\). Entonces \(W(\underset{\sim}{x})\) alcanza la cota de Cramér-Rao si y sólo si existe una función \(a(\theta)\) tal que se tiene la igualdad \[ \begin{align} a(\theta)(W(\underset{\sim}{x})-\tau(\theta))=\frac{\partial}{\partial\theta}\log L(\theta|\underset{\sim}{x}) \mbox{ para todo $\theta$} \end{align} \]

Además, el enunciado del corolario Corollary 3.2 ocurre sí y sólo sí existen funciones \(h(\theta)\), \(k(\theta)\) y \(u(\underset{\sim}{x})\) tales que

\[ \begin{align} L(\theta|\underset{\sim}{x})=u(\underset{\sim}{x})h(\theta)\exp\{W(\underset{\sim}{x})k(\theta)\} \end{align} \] es decir, si y sólo si la distribución de partida pertenece a la familia exponencial.

3.11 Evaluación de estimadores

3.11.1 Mejor estimador insesgado.

Teorema de Rao-Blackwell. Teorema de Lehmann-Scheffé

Si buscamos estimadores insesgados con varianzas pequeñas, podemos restringir nuestra búsqueda a estimadores que sean funciones de estadísticos suficientes.

Theorem 3.6 Si \(Y_1\) y \(Y_2\) son dos variables aleatorias, entonces

\(E(Y_1) = E[E(Y_1 \mid Y_2)]\),

donde en el lado derecho de la ecuación el valor esperado interior es con respecto a la distribución condicional de \(Y_1\) dada \(Y_2\), y el valor esperado exterior es con respecto a la distribución de \(Y_2\).


Theorem 3.7 Si \(Y_1\) y \(Y_2\) representan variables aleatorias, entonces

\(V(Y_1) = E[V(Y_1 \mid Y_2)] + V[E(Y_1 \mid Y_2)]\).

Theorem 3.8 El Teorema de Rao-Blackwell Sea \(\hat{\theta}\) un estimador insesgado para \(\theta\) tal que \(V (\hat{\theta}) <\infty\) . Si \(U\) es un estadístico suficiente para \(\theta\), definamos \(\hat{\theta}^* = E(\hat{\theta} |U)\). Entonces, para toda \(\theta\), \[E(\hat{\theta}^*) =\theta\mbox{ y }V(\hat{\theta}^*) \leq V(\hat{\theta})\].

Proof. Como \(U\) es suficiente para \(\theta\), la distribución condicional de cualquier estadístico (incluyendo \(\hat{\theta}\)), dada \(U\) no depende de \(\theta\). Entonces, \(\hat{\theta}^*=E(\hat{\theta}|U)\) no es función de \(\theta\) y es por tanto un estadístico.

Como \(\hat{\theta}\) es un estimador insesgado para \(\theta\) entonces de los Teoremas Theorem 3.6 y Theorem 3.7, tenemos \[E(\hat{\theta}^*)=E[E(\hat{\theta}|U)]=E(\hat{\theta})=\theta.\] Entonces \(\hat{\theta}^*\) es un estimador insesgado para \(\theta\). El Teorema Theorem 3.7 implica que \[\begin{eqnarray*} V(\hat{\theta})&=&V[E(\hat{\theta}|U)]+E[V(\hat{\theta}|U)]\\ &=&V(\hat{\theta}^*)+E[V(\hat{\theta}|U)]. \end{eqnarray*}\] Como \(V(\hat{\theta} |U = u) \geq 0\) para toda \(u\), se deduce que \(E[V(\hat{\theta} |U)]\geq 0\) y por lo tanto que \[\begin{eqnarray*} V(\hat{\theta})&=&V(\hat{\theta}^*)+E[V(\hat{\theta}|U)]\\ &\geq&V(\hat{\theta}^*)+0 \end{eqnarray*}\] luego \(V (\hat{\theta}) \geq V (\hat{\theta}^*)\).

El Teorema de Rao-Blackwell implica que un estimador insesgado para \(\theta\) con una varianza pequeña es, o puede llegar a ser, una función de un estadístico suficiente. Si \(U\) es un estadístico suficiente para \(\theta\), \(E(\hat{\theta}|U)=\hat{\theta}^*\) entonces \(E(\hat{\theta}^*)=\theta\) y \(V(\hat{\theta}^*)\leq V(\hat{\theta})\). \(\theta^*\) Estimador insesgado, \(\hat{\theta}^*=h(U)\) entonces \(E(h(U)|U)=h(U)=\hat{\theta}^*\).

En general, \(E(h(U)|U)=h(U)\), vemos que usando de nuevo el teorema de Rao-Blackwell nuestro nuevo estimador es simplemente \(h(U)=\hat{\theta}^*\). No ganamos nada después de la primera aplicación.

Debido a que numerosos estadísticos son suficientes para un parámetro \(\theta\) asociado con una distribución ¿qué estadístico suficiente debemos usar cuando aplicamos este teorema?

El criterio de factorización de manera típica identifica un estadístico \(U\) que mejor resume la información de los datos acerca del parámetro \(\theta\). Tales estadísticos reciben el nombre de estadísticos suficinetes mínimos.

Si aplicamos el Teoremo Theorem 3.8 usando \(U\), no sólo obtenemos un estimador con varianza más pequeña si no un estimador insesgado para \(\theta\) con varianza mínima, este se le llama Estimador Insesgado de Varianza Mínima EIVM (MVUE).

Note 3.3

El Teorema de Rao–Blackwell es uno de esos resultados hermosos que nos dicen: “Siempre puedes mejorar tus estimadores, si aprovechas la información disponible en un estadístico suficiente”.

Supongamos que queremos estimar la media \(\theta = \mu\) de una distribución Bernoulli(\(p\)).

Estimador inicial ingenuo:
Tomemos una muestra de tamaño \(n\) y digamos que solo miramos la primera observación \(X_1\).
Como \(E(X_1) = p\), este es un estimador insesgado, pero claramente tiene varianza muy grande.

Estadístico suficiente:
La suma \(U = \sum_{i=1}^n X_i\) es suficiente para \(p\).

Estimador Rao–Blackwellizado:
Calculamos:

\[\hat{p}^* = E(X_1 \mid U) = \frac{U}{n}.\]

Es decir, hemos pasado de un estimador ingenuo (\(X_1\)) a la media muestral, que sabemos es mucho más eficiente.

## Teorema de Rao–Blackwell: ejemplo Bernoulli ----

set.seed(123)

# Parámetros verdaderos
p     <- 0.6
n     <- 20         # tamaño de muestra
Nsim  <- 10000      # número de simulaciones

# Vectores para guardar los estimadores
est_naive <- numeric(Nsim)  # usa solo la primera observación X1
est_RB    <- numeric(Nsim)  # E[X1 | U] = U/n  (media muestral)

for (s in 1:Nsim) {
  X <- rbinom(n, size = 1, prob = p)
  U <- sum(X)

  est_naive[s] <- X[1]
  est_RB[s]    <- U / n
}

# Resumen empírico
mean_naive <- mean(est_naive)
var_naive  <- var(est_naive)

mean_RB <- mean(est_RB)
var_RB  <- var(est_RB)

# Valores teóricos
var_naive_theo <- p * (1 - p)            # Var(X1)
var_RB_theo    <- p * (1 - p) / n        # Var(U/n)

cat("==== Resultados (empíricos vs teóricos) ====\n")
==== Resultados (empíricos vs teóricos) ====
cat(sprintf("E[X1]            empírico = %.4f  (teórico = %.4f)\n", mean_naive, p))
E[X1]            empírico = 0.6006  (teórico = 0.6000)
cat(sprintf("Var(X1)          empírico = %.4f  (teórico = %.4f)\n\n", var_naive, var_naive_theo))
Var(X1)          empírico = 0.2399  (teórico = 0.2400)
cat(sprintf("E[U/n]           empírico = %.4f  (teórico = %.4f)\n", mean_RB, p))
E[U/n]           empírico = 0.6004  (teórico = 0.6000)
cat(sprintf("Var(U/n)         empírico = %.4f  (teórico = %.4f)\n\n", var_RB, var_RB_theo))
Var(U/n)         empírico = 0.0119  (teórico = 0.0120)
cat(sprintf("Reducción de varianza (empírica): Var(U/n) / Var(X1) = %.4f\n", var_RB / var_naive))
Reducción de varianza (empírica): Var(U/n) / Var(X1) = 0.0494
cat(sprintf("Reducción de varianza (teórica) : Var(U/n) / Var(X1) = %.4f\n", var_RB_theo / var_naive_theo))
Reducción de varianza (teórica) : Var(U/n) / Var(X1) = 0.0500

Si empezamos con un estimador insesgado para \(\theta\) y el estadístico suficiente obtenido por medio del criterio de factorización, la aplicación del Teorema de Rao-Blackwell en general lleva a un MVUE para el parámetro.

Sea \(\hat{\theta}\) el parámetro insesgado para \(\theta\). \[E(\hat{\theta})=\theta\] \(U\) igual al obtenido por el método de factorización \[\hat{\theta}^*=E(\hat{\theta}|U)\] donde \(\hat{\theta}^*\) es el MVUE (Aplicar Blackwell).

La consecuencia fundamental de este teorema es que en la búsqueda del estimador UMVUE, basta con restringirnos a aquellos estimadores insesgados que son función de un estadístico suficiente: si trabajamos con un estadístico insesgado que no es función de uno suficiente, tomando esperanzas condicionadas podemos conseguir otro que sea al menos tan bueno como el anterior y sea función del estadístico suficiente. Este proceso se llama a veces Rao-Blackwellización.

3.11.2 Otra versión del teorema de Rao-Blackwell

Theorem 3.9 Otra versión del teorema de Rao-Blackwell Sea una m.a.s. de \(X\), con densidad (o masa de probabilidad) \(f(x|\theta)\). Sea \(T(\underset{\sim}{X})\) un estadístico suficiente para \(\theta\) y sea \(W(\underset{\sim}{X})\) un estimador insesgado de \(\tau(\theta)\). Definimos \[\begin{align} W_T=E_\theta(W|T) \end{align}\] Entonces,

  • \(W_T\) es función unicamente de \(T(\underset{\sim}{X})\) (Es decir, no depende de \(\theta\) y no depende de la muestra \(\underset{\sim}{X}\) sólo a través del valor de \(T(\underset{\sim}{X})\))
  • \(E_\theta(W_T)=\tau(\theta)\).
  • \(V_\theta(W_T)\leq V_\theta(W)\) .

La consecuencia fundamental de este teorema es que en la búsqueda del estimador UMVUE, basta con restringirnos a aquellos estimadores insesgados que son función de un estadístico suficiente: si trabajamos con un estadístico insesgado que no es función de uno suficiente, tomando esperanzas condicionadas podemos conseguir otro que es al menos tan bueno como el anterior y es función del estadístico suficiente. Este proceso se llama a veces Rao-Blackwellización.

3.12 Comportamiento asintótico

3.12.1 Consistencia

Definition 3.4 Se dice que el estimador \(\hat{\theta}_n\) es un estimador consistente de \(\theta\) si, para cualquier número positivo \(\varepsilon\), \[\lim_{n\rightarrow\infty}P(|\hat{\theta}_n-\theta|\leq\varepsilon)=1\] o bien, de forma equivalente, \[\lim_{n\rightarrow\infty}P(|\hat{\theta}_n-\theta|>\varepsilon)=0\]

La consistencia es la misma convergencia en probabilidad!!

Theorem 3.10 Un estimador insesgado \(\hat{\theta}_n\) para \(\theta\) es un estimador consistente de \(\theta\) si \[\lim_{n\rightarrow\infty}V(\hat{\theta}_n)=0.\]

Proof. \(E(\hat{\theta}_n)=\theta\). Sea \(\sigma_{\hat{\theta}_n}=\sqrt{V(\hat{\theta}_n)}\) aplicando el Teorema Tchebysheff para la v.a \(\hat{\theta}_n\), obtenemos \[P(|\hat{\theta}_n-\theta|>k\sigma_{\hat{\theta}_n} )\leq \frac{1}{k^2}\] Para cualquier número positivo \(\varepsilon\), y \(n\) cualquier tamaño muestral \[\varepsilon=k\sigma_{\hat{\theta}_n}\Rightarrow k=\frac{\varepsilon}{\sigma_{\hat{\theta}_n}}\] luego \(k>0\).

La aplicación del teorema de Chebyshev para esta \(n\) fija y esta selección de \(k\) muestra que \[\begin{eqnarray*} P(|\hat{\theta}_n-\theta|>\varepsilon)&=&P\left(|\hat{\theta}_n-\theta|>\left[\frac{\varepsilon}{\sigma_{\hat{\theta}_n}}\right]\sigma_{\hat{\theta}_n}\right)\\ &\leq&\frac{1}{(\varepsilon/\sigma_{\hat{\theta}_n})^2}\\ &=&\frac{V(\hat{\theta}_n)}{\varepsilon^2} \end{eqnarray*}\] Entonces, para cualquier \(n\) fija \[0\leq P(|\hat{\theta}_n-\theta|>\varepsilon)\leq \frac{V(\hat{\theta}_n)}{\varepsilon^2}\] Si \(\lim_{n\rightarrow\infty}V(\hat{\theta}_n)=0\) entonces \[\lim_{n\rightarrow\infty}0\leq \lim_{n\rightarrow\infty}P(|\hat{\theta}_n-\theta|>\varepsilon)\leq\lim_{n\rightarrow\infty} \frac{V(\hat{\theta}_n)}{\varepsilon^2}=0\] Luego \[\lim_{n\rightarrow\infty}P(|\hat{\theta}_n-\theta|>\varepsilon)=0\] luego \(\hat{\theta}_n\) es un estimador consistente para \(\theta\).

Theorem 3.11 \(\displaystyle{\hat{\theta}_n \stackrel{P}{\rightarrow} \theta}\) y \(\displaystyle{\hat{\theta}'_n \stackrel{P}{\rightarrow} \theta'}\) entonces

  1. \(\displaystyle{\hat{\theta}_n +\hat{\theta}'_n \stackrel{P}{\longrightarrow} \theta+\theta'}\).
  2. \(\displaystyle{\hat{\theta}_n \times\hat{\theta}'_n \stackrel{P}{\longrightarrow} \theta\times\theta'}\).
  3. Si \(\theta'\neq 0\), \(\displaystyle{\frac{\hat{\theta}_n}{\hat{\theta}'_n}\stackrel{P}{\longrightarrow} \frac{\theta}{\theta'}}\).
  4. Si \(g (\cdot)\) es una funci'on de valor real que es continua en \(\theta\), entonces \(g (\hat{\theta}_n )\) converge en probabilidad en \(g (\theta)\).
Tip 3.1
  1. El teorema de Rao-Blackwell establece que basta con buscar el estimador UMVUE entre aquellos estimadores que son función de un estadístico suficiente.
  2. Bajo ciertas condiciones (existencia de estadísticos suficientes y completos y de estimadores insesgados), esta combinación de los conceptos de estadístico completo y de estadístico suficiente garantiza la existencia de estimadores UMVUE de una función \(\tau(\theta)\) del parámetro y da un método para construirlos. El siguiente teorema establece este resultado. Podemos decir que este teorema resuelve teóricamente el problema de la estimación puntual, entendida ésta como la búsqueda del UMVUE.

Definition 3.6 (estadístico completo) Sea \(f_T(t|\theta)\) la función de densidad (o de masa de probabilidad) de un estadístico \(T\) . Diremos que la familia de distribuciones \({f_T(t|\theta) : \theta\in \Theta}\) es completa si se da la implicación siguiente: \[\begin{align} E_\theta(g(T))=0 \mbox{para todo $\theta$ entonces} P_\theta(g(T) = 0) = 1 \mbox{para todo $\theta$}. \end{align}\]

En ese caso diremos que \(T\) es un estadístico completo.

## Completitud para T ~ Binomial(n, p)
set.seed(1)

n     <- 10
gridp <- seq(0.1, 0.9, by = 0.1)
N     <- 100000  # simulaciones por p

# Dos funciones no triviales de T (no dependen de p)
g1 <- function(T) T - n/2
g2 <- function(T) (T - n/2) * (T - n/3)

# Función trivial (cero) para contraste
g0 <- function(T) 0 * T

sim_mean_g <- function(gfun) {
  sapply(gridp, function(p) {
    T <- rbinom(N, size = n, prob = p)
    mean(gfun(T))
  })
}

# Medias empíricas
m0 <- sim_mean_g(g0)
m1 <- sim_mean_g(g1)
m2 <- sim_mean_g(g2)

# Valores teóricos (para comparar)
# E[T] = n p; Var[T] = n p (1-p); E[T^2] = Var[T] + (E[T])^2
theo_g1 <- n * (gridp - 1/2)  # E[g1(T)] = n(p - 1/2)

ET   <- n * gridp
VarT <- n * gridp * (1 - gridp)
ET2  <- VarT + ET^2
theo_g2 <- ET2 - (n/2 + n/3) * ET + (n^2)/(6)  # expandir (T - n/2)(T - n/3)

# Mostrar resultados
tab <- data.frame(
  p = gridp,
  E_g0_emp = round(m0, 5),
  E_g1_emp = round(m1, 5),
  E_g1_teo = round(theo_g1, 5),
  E_g2_emp = round(m2, 5),
  E_g2_teo = round(theo_g2, 5)
)
print(tab, row.names = FALSE)
   p E_g0_emp E_g1_emp E_g1_teo E_g2_emp E_g2_teo
 0.1        0 -4.00184       -4 10.24421 10.23333
 0.2        0 -2.99891       -3  5.61325  5.60000
 0.3        0 -2.00216       -2  2.77471  2.76667
 0.4        0 -1.00631       -1  1.74363  1.73333
 0.5        0 -0.00727        0  2.50485  2.50000
 0.6        0  1.00466        1  5.09209  5.06667
 0.7        0  1.99571        2  9.43229  9.43333
 0.8        0  3.00090        3 15.55946 15.60000
 0.9        0  3.99831        4 23.58862 23.56667
Tip 3.2: Qué observar
  • Para \(g_0(T)\equiv 0\), la media empírica es \(\approx 0\) para todos los \(p\).

  • Para \(g_1(T) = T - \tfrac{n}{2}\), la esperanza teórica es \(n(p - \tfrac{1}{2})\):
    solo vale \(0\) cuando \(p = \tfrac{1}{2}\); no es \(0\) para los demás \(p\).

  • Para \(g_2(T) = (T - \tfrac{n}{2})(T - \tfrac{n}{3})\), la esperanza es un polinomio en \(p\) que no es idénticamente cero; de nuevo, no se anula para todos los \(p\).


Esto ilustra la completitud: si alguna \(g(T)\) (que no depende de \(p\)) tuviera
\(E_p[g(T)] = 0\) para todo \(p \in (0,1)\), entonces esa \(g\) debe ser (casi seguro) la función cero.

➡️ De ahí que, con suficiencia + completitud, el estimador insesgado derivado es el único UMVUE (Teorema de Lehmann–Scheffé).

La definición de completitud refuerza la de suficiencia en el sentido de que si un estadístico es suficiente y completo entonces, es suficiente minimal (el recíproco no es cierto)

Theorem 3.12 (Teorema de Lehmann-Scheffé) Si \(T(\underset{\sim}{X})\) es un estadístico suficiente y completo para \(\theta\) y \(W(\underset{\sim}{X})\) es un estimador insesgado cualquiera de \(\tau(\theta)\), entonces \[\begin{align} W_T(\underset{\sim}{X})=E_\theta(W|T) \end{align}\] es el mejor estimador insesgado (UMVUE) de \(\tau(\theta)\). Si, además, \(V(W_T) < \infty\) para todo \(\theta\), entonces \(W_T\) es único.

  1. La demostración del teorema de Lehmann-Scheffé se basa en el hecho de que, si existen estimadores insesgados, esencialmente sólo existe uno que sea función del estadístico suficiente y completo, pues condicionando cualquiera de los insesgados al estadístico suficiente y completo se obtiene siempre el mismo resultado.
  2. El teorema de Rao-Blackwell garantiza que al tomar esperanzas condicionadas se ha reducido la varianza, llegando así al UMVUE.
  3. La principal conclusión del teorema de Lehmann-Scheffé es que si existe un estimador insesgado de \(\tau(\theta)\) que sea función de un estadístico suficiente y completo, entonces es el único UMVUE de \(\tau(\theta)\).
Tip 3.3: Recreando Lehmann–Scheffé en R

¡Vamos a “recrear” Lehmann–Scheffé en R con un ejemplo clásico y totalmente replicable!

Usaremos \(X_1,\dots,X_n \sim \text{i.i.d. Poisson}(\lambda)\).
Entonces

\[ T=\sum_{i=1}^n X_i \sim \text{Poisson}(n\lambda) \]

es suficiente y completo para \(\lambda\).

Veremos dos blancos \(\tau(\theta)\):

  • \(\tau(\lambda)=\lambda\)
  • \(\tau(\lambda)=e^{-\lambda}\)

En cada caso partimos de un estimador insesgado cualquiera \(W(\mathbf X)\) y lo “Rao–Blackwellizamos” con \(T\):

\[ W_T(\mathbf X)=\mathbb{E}_\theta[W \mid T]. \]

Por Lehmann–Scheffé, \(W_T\) es el UMVUE.

Idea teórica rápida

  • Para \(\tau(\lambda)=\lambda\):
    Un estimador insesgado simple es \(W=X_1\) (pues \(\mathbb{E}[X_1]=\lambda\)).
    Condicionando en \(T\):

\[ W_T = \mathbb{E}[X_1 \mid T] = \tfrac{T}{n} \quad \Longrightarrow \quad \text{UMVUE}. \]

  • Para \(\tau(\lambda)=e^{-\lambda}\):
    Un insesgado simple es \(W=\mathbf 1\{X_1=0\}\) (pues \(\Pr(X_1=0)=e^{-\lambda}\)).
    Dado \(T=t\),

\[ (X_1,\dots,X_n)\mid T=t \sim \text{Multinomial}\Big(t;\tfrac{1}{n},\dots,\tfrac{1}{n}\Big), \]

así que

\[ \Pr(X_1=0 \mid T=t)=\Big(1-\tfrac{1}{n}\Big)^t. \]

Por ende,

\[ W_T = \mathbb{E}[\,\mathbf 1\{X_1=0\}\mid T] = \Big(1-\tfrac{1}{n}\Big)^T, \]

que es el UMVUE de \(e^{-\lambda}\).

## Lehmann–Scheffé con Poisson ----
## X_i ~ Poisson(lambda).  T = sum X_i ~ Poisson(n*lambda) es suficiente y completo.

set.seed(123)

# Parámetros
lambda <- 2
n      <- 10
Nsim   <- 20000

# Contenedores
W1_naive  <- numeric(Nsim)  # para tau1(lambda)=lambda, W = X1
W1_RB     <- numeric(Nsim)  # W_T = E[X1|T] = T/n

W2_naive  <- numeric(Nsim)  # para tau2(lambda)=exp(-lambda), W = 1{X1=0}
W2_RB     <- numeric(Nsim)  # W_T = E[1{X1=0} | T] = ((n-1)/n)^T

for (s in 1:Nsim) {
  X <- rpois(n, lambda)
  Tsum <- sum(X)
  
  # Caso 1: tau(lambda) = lambda
  W1_naive[s] <- X[1]
  W1_RB[s]    <- Tsum / n
  
  # Caso 2: tau(lambda) = exp(-lambda)
  W2_naive[s] <- as.numeric(X[1] == 0)
  W2_RB[s]    <- (1 - 1/n)^Tsum
}

## Resúmenes empíricos
res <- data.frame(
  Estimador = c("W = X1", "W_T = T/n", "W = 1{X1=0}", "W_T = ((n-1)/n)^T"),
  Objetivo  = c("lambda", "lambda", "exp(-lambda)", "exp(-lambda)"),
  Media     = c(mean(W1_naive), mean(W1_RB), mean(W2_naive), mean(W2_RB)),
  Varianza  = c(var(W1_naive), var(W1_RB), var(W2_naive), var(W2_RB))
)

print(res, row.names = FALSE)
         Estimador     Objetivo    Media    Varianza
            W = X1       lambda 2.015250 2.036219248
         W_T = T/n       lambda 1.999100 0.197314056
       W = 1{X1=0} exp(-lambda) 0.131500 0.114213461
 W_T = ((n-1)/n)^T exp(-lambda) 0.135274 0.003988959
## Valores teóricos para comparación
# Caso 1 (lambda):
var_W1_naive_theo <- lambda            # Var(X1)
var_W1_RB_theo    <- lambda / n        # Var(T/n)

# Caso 2 (exp(-lambda)):
# Var(1{X1=0}) = e^{-lambda}(1 - e^{-lambda})
var_W2_naive_theo <- exp(-lambda) * (1 - exp(-lambda))

# Var(((n-1)/n)^T) = E[a^{2T}] - (E[a^T])^2, con a = (n-1)/n y T~Poisson(n*lambda)
a <- 1 - 1/n
E_aT     <- exp(n * lambda * (a - 1))        # = exp(n*lambda*(-1/n)) = exp(-lambda)
E_a2T    <- exp(n * lambda * (a^2 - 1))      # = exp(-2*lambda + lambda/n)
var_W2_RB_theo <- E_a2T - E_aT^2             # = e^{-2lambda}(e^{lambda/n} - 1)

cat("\n--- Teórico ---\n")

--- Teórico ---
cat(sprintf("Var(W=X1)                 = %.5f\n", var_W1_naive_theo))
Var(W=X1)                 = 2.00000
cat(sprintf("Var(W_T=T/n)              = %.5f\n", var_W1_RB_theo))
Var(W_T=T/n)              = 0.20000
cat(sprintf("Var(W=1{X1=0})            = %.5f\n", var_W2_naive_theo))
Var(W=1{X1=0})            = 0.11702
cat(sprintf("Var(W_T=((n-1)/n)^T)      = %.5f\n\n", var_W2_RB_theo))
Var(W_T=((n-1)/n)^T)      = 0.00406
cat("Relaciones esperadas (Rao–Blackwell):\n")
Relaciones esperadas (Rao–Blackwell):
cat(sprintf("Var(T/n)  <= Var(X1):           %.5f <= %.5f\n", var_W1_RB_theo, var_W1_naive_theo))
Var(T/n)  <= Var(X1):           0.20000 <= 2.00000
cat(sprintf("Var(((n-1)/n)^T) <= Var(1{X1=0}): %.5f <= %.5f\n", var_W2_RB_theo, var_W2_naive_theo))
Var(((n-1)/n)^T) <= Var(1{X1=0}): 0.00406 <= 0.11702

Qué comprueba el script

  • Insesgadez: las medias empíricas de \(W\) y \(W_T\) aproximan \(\tau(\lambda)\) (ya sea \(\lambda\) o \(e^{-\lambda}\)).

  • Mejor varianza:
    \[\mathrm{Var}(W_T) \leq \mathrm{Var}(W)\]
    en ambos casos (se ve empíricamente y con fórmulas teóricas).

  • Forma cerrada del UMVUE:

    • Para \(\tau(\lambda)=\lambda\): \(W_T = T/n\).
    • Para \(\tau(\lambda)=e^{-\lambda}\): \(W_T = \Big(\tfrac{n-1}{n}\Big)^T\).

Y, por Lehmann–Scheffé, como \(T\) es suficiente y completo y \(\mathrm{Var}(W_T)<\infty\), este UMVUE es único.

# --- Preparación ----
# install.packages("ggplot2") # si no lo tienes
library(ggplot2)

explore_LS <- function(lambda = 2, n = 10, Nsim = 20000, seed = 123,
                       show_plots = TRUE, save_plots = FALSE, prefix = "LS") {
  set.seed(seed)

  # ---- Simulación ----
  # Generamos Nsim muestras, cada una de tamaño n (Poisson(lambda))
  # Para eficiencia, simulamos Nsim*n y damos forma de matriz.
  Xmat <- matrix(rpois(Nsim * n, lambda), nrow = Nsim, ncol = n)
  Tsum <- rowSums(Xmat)

  # Caso 1: tau(lambda) = lambda
  W1  <- Xmat[, 1]          # estimador ingenuo: X1
  WT1 <- Tsum / n           # Rao-Blackwell: E[X1 | T] = T/n  (UMVUE)

  # Caso 2: tau(lambda) = exp(-lambda)
  W2  <- as.numeric(Xmat[, 1] == 0)        # ingenuo: 1{X1=0}
  WT2 <- (1 - 1/n)^Tsum                    # RB: ((n-1)/n)^T  (UMVUE)

  # ---- Resúmenes empíricos ----
  tab <- data.frame(
    Estimador = c("W = X1", "W_T = T/n", "W = 1{X1=0}", "W_T = ((n-1)/n)^T"),
    Objetivo  = c("lambda", "lambda", "exp(-lambda)", "exp(-lambda)"),
    Media     = c(mean(W1), mean(WT1), mean(W2), mean(WT2)),
    Varianza  = c(var(W1),  var(WT1), var(W2), var(WT2))
  )

  # ---- Valores teóricos ----
  # Caso 1:
  var_W1_theo  <- lambda           # Var(X1)
  var_WT1_theo <- lambda / n       # Var(T/n)

  # Caso 2:
  # Var(1{X1=0}) = e^{-lambda}(1 - e^{-lambda})
  var_W2_theo  <- exp(-lambda) * (1 - exp(-lambda))

  # Var(((n-1)/n)^T) = E[a^{2T}] - (E[a^T])^2 con T ~ Poisson(n*lambda), a=(n-1)/n
  a <- 1 - 1/n
  E_aT  <- exp(n * lambda * (a - 1))           # = exp(-lambda)
  E_a2T <- exp(n * lambda * (a^2 - 1))         # = exp(-2*lambda + lambda/n)
  var_WT2_theo <- E_a2T - E_aT^2               # = e^{-2lambda}(e^{lambda/n} - 1)

  tab_teo <- data.frame(
    Estimador = c("W = X1", "W_T = T/n", "W = 1{X1=0}", "W_T = ((n-1)/n)^T"),
    Var_teor  = c(var_W1_theo, var_WT1_theo, var_W2_theo, var_WT2_theo)
  )

  # ---- Data para gráficas ----
  df1 <- rbind(
    data.frame(valor = W1,  tipo = "W = X1",        objetivo = "lambda"),
    data.frame(valor = WT1, tipo = "W_T = T/n",     objetivo = "lambda")
  )

  df2 <- rbind(
    data.frame(valor = W2,  tipo = "W = 1{X1=0}",           objetivo = "exp(-lambda)"),
    data.frame(valor = WT2, tipo = "W_T = ((n-1)/n)^T",     objetivo = "exp(-lambda)")
  )

  # ---- Gráficas ----
  p_dens_1 <- ggplot(df1, aes(x = valor, fill = tipo)) +
    geom_density(alpha = 0.35, adjust = 1.1) +
    geom_vline(xintercept = lambda, linetype = 2) +
    labs(
      title = bquote("Distribuciones de " * W ~ "vs" ~ W[T] ~ " para " ~ tau(lambda) == lambda),
      subtitle = bquote(lambda == .(lambda) ~ ", " ~ n == .(n) ~ ", " ~ N[sim] == .(Nsim)),
      x = "valor", y = "densidad"
    ) +
    theme_minimal() +
    theme(legend.position = "top")

  p_box_1 <- ggplot(df1, aes(x = tipo, y = valor, fill = tipo)) +
    geom_boxplot(alpha = 0.4, width = 0.6, outlier.alpha = 0.25) +
    geom_hline(yintercept = lambda, linetype = 2) +
    labs(
      title   = "Boxplots: W vs W_T (objetivo: lambda)",
      x = NULL, y = "valor"
    ) +
    theme_minimal() +
    theme(legend.position = "none")

  # Para el caso indicador es discreto; usamos barras y boxplot (WT2 es continuo en [0,1])
  p_bar_2 <- ggplot(df2, aes(x = factor(round(valor, 3)), fill = tipo)) +
    geom_bar(position = "dodge") +
    labs(
      title = bquote("Distribuciones de " * W ~ "vs" ~ W[T] ~ " para " ~ tau(lambda) == e^{-lambda}),
      subtitle = bquote(lambda == .(lambda) ~ ", " ~ n == .(n) ~ ", " ~ N[sim] == .(Nsim)),
      x = "valor (redondeado a 3 decimales)", y = "frecuencia"
    ) +
    theme_minimal() +
    theme(legend.position = "top")

  p_box_2 <- ggplot(df2, aes(x = tipo, y = valor, fill = tipo)) +
    geom_boxplot(alpha = 0.4, width = 0.6, outlier.alpha = 0.25) +
    geom_hline(yintercept = exp(-lambda), linetype = 2) +
    labs(
      title = "Boxplots: W vs W_T (objetivo: exp(-lambda))",
      x = NULL, y = "valor"
    ) +
    theme_minimal() +
    theme(legend.position = "none")

  if (show_plots) {
    print(p_dens_1); print(p_box_1)
    print(p_bar_2);  print(p_box_2)
  }

  if (save_plots) {
    ggsave(sprintf("%s_dens_lambda.png", prefix), p_dens_1, width = 7, height = 4.5, dpi = 150)
    ggsave(sprintf("%s_box_lambda.png",  prefix), p_box_1,  width = 6, height = 4.5, dpi = 150)
    ggsave(sprintf("%s_bar_expl.png",    prefix), p_bar_2,  width = 7, height = 4.5, dpi = 150)
    ggsave(sprintf("%s_box_expl.png",    prefix), p_box_2,  width = 6, height = 4.5, dpi = 150)
  }

  # ---- Salida ----
  list(
    resumen_empirico = tab,
    resumen_teorico  = tab_teo,
    plots = list(dens_lambda = p_dens_1, box_lambda = p_box_1,
                 bar_expl = p_bar_2, box_expl = p_box_2)
  )
}

# ===== Ejemplo de uso =====
out <- explore_LS(lambda = 2, n = 10, Nsim = 20000, seed = 123, show_plots = TRUE)

# Mira los resúmenes:
out$resumen_empirico
          Estimador     Objetivo     Media    Varianza
1            W = X1       lambda 1.9870500 1.961480372
2         W_T = T/n       lambda 1.9991000 0.199480164
3       W = 1{X1=0} exp(-lambda) 0.1354500 0.117109153
4 W_T = ((n-1)/n)^T exp(-lambda) 0.1354626 0.004083226
out$resumen_teorico
          Estimador    Var_teor
1            W = X1 2.000000000
2         W_T = T/n 0.200000000
3       W = 1{X1=0} 0.117019644
4 W_T = ((n-1)/n)^T 0.004055133

Qué muestra

  • Densidad y boxplot (objetivo \(\lambda\)):
    Compara \(W = X_1\) vs \(W_T = T/n\).
    Ambos se centran en \(\lambda\), pero verás menor varianza para \(W_T\).

  • Barras/boxplot (objetivo \(e^{-\lambda}\)):
    Compara \(W = \mathbf{1}\{X_1=0\}\) vs \(W_T = \Big(\tfrac{n-1}{n}\Big)^T\).
    Ambos centrados en \(e^{-\lambda}\), con varianza menor para \(W_T\).


Sugerencia: juega con \(\lambda\) y \(n\) en explore_LS(lambda, n, ...) para ver cómo cambia la varianza.
Por Lehmann–Scheffé, \(W_T\) es el UMVUE y único (bajo varianza finita) cuando \(T\) es suficiente y completo.

3.12.2 Normalidad asintótica

En muchas ocasiones sólo es posible realizar estudios del comportamiento asintótico (cuando n tiende a infinito) de los estimadores. Ya hemos estudiado una propiedad asintótica: la consistencia. Veremos ahora que es posible medir la velocidad de convergencia de estimadores consistentes y así seleccionar los que convergen al verdadero valor del parámetro más rápidamente.

Example 3.11 Sea \(X_1,\ldots, X_n\) una m.a.s. de \(X\sim Pois(\lambda)\), con \(\lambda>0\). En este modelo, el estimador de momentos de \(\lambda\) coincide con el máximo verosimil: \(\hat{\lambda}_n=\bar{X}_n\). Para determinar la distribución del estimador \(\hat{\lambda}\) resulta mucho más útil aproximarla por una distribución más sencilla a la que se acerca asintóticamente.

La versión del teorema central del límite para variables aleatorias independientes e idénticamente distribuidas puede aplicarse porque \(V(X) = \lambda < \infty\).

\[\begin{align} \frac{\sqrt{n}(\hat{\lambda}_n-\lambda)}{\sqrt{\lambda}}\longrightarrow N(0,1) \mbox{ débilmente,} \end{align}\] es decir, para todo \(\lambda\in \Theta\) y para todo \(w\in\mathbb{R}\), \[\begin{align} P_\lambda(\hat{\lambda}\leq w)\approx \phi\left(\frac{\sqrt{n}(w-\lambda)}{\sqrt{\lambda}}\right) \end{align}\] donde \(\phi\) es la función de distribución de la normal estándar. La aproximación es tanto mejor cuanto mayores son \(n\) o \(\lambda\).

Obsérvese que \(\hat{\lambda}_n\) es consistente pues, por las leyes de los grandes números, \[\hat{\lambda}_n = \bar{X}_n \stackrel{P}{\longrightarrow} E(X) = \lambda.\] Así, \(\hat{\lambda}_n -\lambda \stackrel{P}{\longrightarrow} 0\) y también en distribución. Esta convergencia a la distribución degenerada en 0 no nos informa de la velocidad a la que \(\hat{\lambda}_n\) se acerca a \(\lambda\) ni de cómo lo hace.

El hecho de que \(V(\sqrt{n}(\hat{\lambda}_n -\lambda)) = \lambda\) para todo n indica que la velocidad a la que \(\hat{\lambda}_n\) se acerca a \(\lambda\) es la misma con la que \(\frac{1}{\sqrt{n}}\) se acerca a 0: multiplicar por \(\sqrt{n}\) es la forma de estabilizar las diferencias \((\hat{\lambda}_n-\lambda)\), es la estandarización adecuada.

El resultado derivado del teorema central del límite, la distribución asintótica de \(\sqrt{n}(\hat{\lambda}_n -\lambda)\) es \(N(0,\lambda)\), responde a la pregunta de cómo es la aproximación \(\hat{\lambda}_n\) a \(\lambda\) : los valores del estimador se distribuyen alrededor del verdadero valor del parámetro igual que los valores de una variable aleatoria \(N(0,\lambda)\) se distribuyen alrededor de 0.

Definition 3.5 En la práctica la gran mayoría de los estimadores usuales, convenientemente centrados y normalizados, tienen distribución asintótica normal. Se dice que presentan y se denota \[\begin{align} \hat{\theta}_n\sim AN(\theta,v_n) \end{align}\] cuando \[\begin{align} \frac{1}{\sqrt{v_n}}(\hat{\theta}_n-\theta)\stackrel{d}{\longrightarrow}N(0,1) \end{align}\]

A la cantidad \(v_n\) se la llama varianza asintótica de \(\hat{\theta}_n\). El teorema central del límite es el responsable de la normalidad asintótica de muchos estimadores.

3.12.3 Método delta

En muchos casos, sólo será de interés el comportamiento del estimador alrededor del verdadero valor del parámetro. Si además el estimador es una función suave de un estadístico cuyo comportamiento asintótico es conocido, esa función podrá linealizarse en un entorno del verdadero valor del parámetro, lo cuál facilitará enormemente el estudio asintótico del estimador.

Example 3.12 Queremos estimar \(\theta=P(X=0) = \exp\{-\lambda\}\). Por el principio de invariancia, el estimador máximo verosímil de \(\theta\) es \(\hat{\theta}_n = \exp\{-\bar{X}_n\}\) , dado que \(\bar{X}_n\) es el estimador máximo verosímil de \(\lambda\). Por otro lado, \(\bar{X}_n\) es consistente para \(\lambda\) y \(g(\lambda) = \exp\{-\lambda\}\) es una función continua, luego \(\hat{\theta}_n\) es consistente para \(\theta\). Estamos interesados ahora en encontrar la distribución asintótica de \[\begin{align} \sqrt{n}(\hat{\theta}_n-\theta)=\sqrt{n}(\exp\{-\bar{X}_n\}-\exp\{-\lambda\}) \end{align}\] La herramienta en la que nos basaremos para hallar esa distribución asintótica es el método delta.

Theorem 3.13 (Método Delta) Sea \(\{a_n\}_n\) una sucesión de número reales tales que \(a_n\rightarrow\infty\) cuando \(n\rightarrow\infty\) y con \(a_n\neq0\) para todo \(n\). Sea \(\hat{\theta}_n\) una sucesión de estimadores de \(\theta\) tales que \[\begin{align} a_n(\hat{\theta}_n-\theta)\stackrel{d}{\longrightarrow}N(0,\sigma_\theta^2) \end{align}\] y sea \(g(x)\) una función con primera derivada continua en un intervalo que contiene a \(\theta\). \[\begin{align} a_n(g(\hat{\theta}_n)-g(\theta))\stackrel{d}{\longrightarrow}N(0,(g^{'}(\theta))^2\sigma_\theta^2) \end{align}\]

Example 3.13 Estimamos \(\theta=P(X=0) = \exp\{-\lambda\}\) mediante \(\hat{\theta}_n = \exp\{-\bar{X}_n\}\). Por otra parte, \(\sqrt{n}(\hat{\lambda}_n-\lambda)\stackrel{d}{\longrightarrow}N(0,\lambda)\). Además \(g(\lambda)=\exp\{-\lambda\}\) es derivable con derivada continua \(g^{'}(\lambda)=-\exp\{-\lambda\}\)

Aplicando el método delta. \[\begin{align} \sqrt{n}(\hat{\theta}_n-\theta)=\sqrt{n}(\exp\{-\bar{X}_n\}-\exp\{-\lambda\})\stackrel{d}{\longrightarrow}N(0,\lambda\exp\{-2\lambda\}) \end{align}\]

3.12.4 Eficiencia relativa asintótica

Sea \(T_n(\underset{\sim}{X}) = T_n(X_1 , \ldots, X_n)\) una sucesión de estimadores de una función \(\tau(\theta)\) que verifica lo siguiente: \[\begin{align} \sqrt{n}(T_n(\underset{\sim}{X})-\tau(\theta))\stackrel{d}{\longrightarrow}N(b(\theta),\sigma^2(\theta)) \end{align}\] Si \(b(\theta) = 0\) diremos que \(T_n(\underset{\sim}{X})\) es asintóticamente insesgado. En caso contrario, diremos que \(T_n(\underset{\sim}{X})\) es asintóticamente sesgado.

Sean dos sucesiones \(T_n(X)\) y \(S_n(X)\) de estimadores de \(\tau(\theta)\) asintóticamente normales:

\[\begin{align} \sqrt{n}(T_n(\underset{\sim}{X})-\tau(\theta))\stackrel{d}{\longrightarrow}N(0,\sigma_T^2(\theta)) \end{align}\]

\[\begin{align} \sqrt{n}(S_n(\underset{\sim}{X})-\tau(\theta))\stackrel{d}{\longrightarrow}N(0,\sigma_S^2(\theta)) \end{align}\] Se define la eficiencia relativa asintótica de \(S_n(X)\) respecto a \(T_n(X)\) como \[\begin{align} ARE(\theta, S_n , T_n)=\frac{\frac{1}{\sigma_S^2(\theta)}}{\frac{1}{\sigma_T^2(\theta)}}=\frac{\sigma_T^2(\theta)}{\sigma_S^2(\theta)} \end{align}\]

El valor de la eficiencia relativa asintótica puede interpretarse como el cociente de los tamaños de muestra necesarios para obtener la misma precisión asintótica (o la misma varianza asintótica) mediante los dos estimadores en la estimación de \(\tau(\theta)\). En efecto, si elegimos tamaño muestral m para \(T\) y n para \(S\), las varianzas asintóticas son, respectivamente, \(\frac{\sigma_T^2(\theta)}{m}\) y \(\frac{\sigma_S^2(\theta)}{n}\). Si forzamos aque ambas sean iguales, se tiene que \[\begin{align} \frac{\sigma_T^2(\theta)}{m}=\frac{\sigma_S^2(\theta)}{n}\Longleftrightarrow\frac{m}{n}=\frac{\sigma_T^2(\theta)}{\sigma_S^2(\theta)}=ARE(\theta, S_n , T_n) \end{align}\]

Es decir, si \(ARE(\theta, S_n , T_m) = 0.5\) entonces \(S\) es menos eficiente que \(T\) asintóticamente: para tener la misma precisión con el estimador \(S\) hace falta una muestra el doble de grande que si utilizásemos \(T\) (\(ARE(\theta, S_n , T_m)= 0.5=\frac{m}{n}\Longrightarrow n = 2m\)).


3.13 Ejercicios

Exercise 3.1 Se cuenta con una muestra de variables aleatorias independientes e idénticamente distribuidas \(x=(x_1,\ldots,x_n)\) como cada uno de los modelos que se indican abajo.

  1. Poisson\((\lambda)\)
  2. Normal con media conocida \(\mu\) y varianza desconocida \(\sigma^2\)
  3. \(f(x;\theta)=\theta x^{\theta-1}I_{\left[0,1\right]}(x)\)

Encuentra un estimador de momentos y el estimador de máxima verosimilitud. Determine si son insesgados asintóticamente (ve la Sección 11.7 del Kalbfleisch donde se define este concepto) o no y suficientes. Encuentra el vector de estadísticas suficientes de menor dimensión. Demuestre si es suficiente minimal o no.

Exercise 3.2 Se toma una muestra de tamaño \(n=1\) de una m.a.s de \(X\sim exp(\lambda)\), es decir \(f(x|\lambda)=\lambda\exp(-\lambda x)I_{\left[0,\infty\right)}(x), \lambda>0\) y se observa \(x=3\). Determine la verosimilitud \(L(\lambda|3)\) y establezca el estimador máximo verosil.

Exercise 3.3 Sea \(X_1,\ldots, X_n\) una sucesión de variables aleatorias iid Bernoulli con parámetro \(\theta\), con \(0<\theta<1\). Muestre que \(T(\underset{\sim}{X})=X_1+\ldots+ X_n\) es un estadístico suficiente para \(\theta\).

Exercise 3.4 Sea \(X_1,\ldots, X_n\) una sucesión de v.a.s iid de \(X\sim N(\mu,\sigma^2)\), donde \(\sigma^2\) es conocida. Utilizando el teorema 1.3 muestre que \(T(\underset{\sim}{X})=T(X_1,\ldots, X_n)=\bar{X}=\frac{\sum_{i=1}^{n}X_i}{n}\) es una estadística suficiente para \(\mu\).

Exercise 3.5 Sea \(X_1,\ldots, X_n\) una sucesión de v.a.s iid de \(X\sim N(\mu,\sigma^2)\), donde \(\sigma^2\) es conocida. Utilizando el teorema 1.3 muestre que \(T(\underset{\sim}{X})=T(X_1,\ldots, X_n)=\bar{X}=\frac{\sum_{i=1}^{n}X_i}{n}\) es una estadística suficiente para \(\mu\).

Exercise 3.6 Sea \(X_1,\ldots, X_n\) una sucesión de v.a.s iid de $X$ definida de la siguiente manera. \[\begin{align} f(x|\theta)= \left\{ \begin{array}{lcc} \frac{1}{\theta} & si & x=1, 2, \ldots, \theta \\ \\ 0 & en & \mbox{cualquier otro caso} \end{array} \right. \end{align}\] Utilice el teorema de factorización y determine que \(T(\underset{\sim}{X})=\max_iX_i\) es un estadístico suficiente para \(\theta\).

Exercise 3.7 Muestre que \(T_2=X_1+X_2\) es un estadístico suficiente para \(\theta\), donde \(\theta\) es el parámetro de la distribución Bernoulli, la cual es la distribución de una muestra aleatoria simple de tamaño dos \(X_1, X_2\).

Exercise 3.8 Muestre que \(T_2^{´}=X_1X_2\) NO es un estadístico suficiente para \(\theta\), donde \(\theta\) es el parámetro de la distribución Bernoulli, la cual es la distribución de una muestra aleatoria simple de tamaño dos \(X_1, X_2\).

Exercise 3.9 Sea \(X_1,\ldots, X_n\) una mas de \(X\sim N(\mu, \sigma^2)\). Demuestre que \(\bar{X}\) y \(S^2\) son estimadores insesgados para \(\mu\) y \(\sigma^2\) respectivamente

Exercise 3.10 Sea \(X_1,\ldots, X_n\) una mas de \(X\sim N(\mu, \sigma^2)\). Muestre que el error cuadrático medio de \(\bar{X}\) y \(S^2\) coindien con su varianza respectivamente.

Exercise 3.11 Sea \(X_1,\ldots, X_n\) m.a.s. de \(X\sim N(\mu, 1)\), \(-\infty<\mu<\infty\). La media \(\bar{X}\) y la mediana muestral \(M_n\) son estimadores insesgados de \(\mu\). Para \(n=21\), muestre que la mediana muestral \(M_n\) es inanmisible para \(\mu\), frente a \(\bar{X}\).

Exercise 3.12 Suponga que \(\hat{\theta}\) es un estimador para un parámetro \(\theta\) y \(E(\hat{\theta}) = a\theta+b\) para algunas constantes diferentes de cero a y b.

  1. En términos de a, b, y \(\theta\), ¿cuál es sesgo \(B(\hat{\theta})\)?

  2. Encuentre una función de \(\hat{\theta}\), por ejemplo \(\hat{\theta}^{*}\) que es un estimador insesgado para \(\theta\).

Exercise 3.13 Usando la identidad \[\begin{align} (\hat{\theta}-\theta) = [\hat{\theta}-E( \hat{\theta})] + [E(\hat{\theta})-\theta] = [\hat{\theta}-E(\hat{\theta})] + B(\hat{\theta}), \end{align}\]

demuestre que \[\begin{align} MSE(\hat{\theta}) = E[(\hat{\theta}-\theta)^2]= V(\hat{\theta})+(B(\hat{\theta}))^2. \end{align}\]

Exercise 3.14 Consulte los Ejercicios en (Exercise 3.13) y considere el estimador insesgado \(\hat{\theta}^*\) que usted propuso en el Ejercicio (Exercise 3.12). Exprese \(MSE(\hat{\theta}^*)\) como función de \(V(\hat{\theta}^*)\).

Exercise 3.15 Suponga que \(E(\hat{\theta}_1) = E(\hat{\theta}_2) = \theta\), \(V(\hat{\theta}_1) = \sigma_1^2\) , y \(V(\hat{\theta}_2)=\sigma_2^2\). Considere el estimador \(\hat{\theta}_3=a\hat{\theta}_1+(1-a)\hat{\theta}_2\).

  1. Demuestre que \(\hat{\theta}_3\) es un estimador insesgado para \(\theta\).

  2. Si \(\hat{\theta}_2\) son independientes, ¿cómo debe escogerse la constante \(a\) para minimizar la varianza de \(\hat{\theta}_3\)?

Exercise 3.16 Suponga que \(E(\hat{\theta}_1) = E(\hat{\theta}_2) = \theta\), \(V(\hat{\theta}_1) = \sigma_1^2\) , y \(V(\hat{\theta}_2)=\sigma_2^2\). Considere el estimador \(\hat{\theta}_3=a\hat{\theta}_1+(1-a)\hat{\theta}_2\).

  1. ¿Cómo debe elegirse la constante a para minimizar la varianza de \(\hat{\theta}_3\) , si \(\hat{\theta}_1\) y \(\hat{\theta}_2\) no son independientes pero son tales que \(Cov(\hat{\theta}_1 , \hat{\theta}_2) = c \neq 0\)?

Exercise 3.17 Suponga que \(Y_1,\, Y_2,\, Y_3\) denotan una muestra aleatoria ind de una distribuci'on exponencial con funci'on de densidad \[f(y)=\left\{\begin{array}{ll}\left(\frac{1}{\theta}\right)e^{-y/\theta},& y>0,\\0,&\mbox{ en cualquier otro punto}\end{array}\right.\] Considere los siguientes estimadores de \(\theta\): \[\hat{\theta}_1= Y_1,\,\,\,\,\,\,\,\,\,\,\hat{\theta}_2=\frac{Y_1+Y_2}{2},\,\,\,\,\,\,\,\,\,\,\hat{\theta}_3=\frac{Y_1+2Y_2}{3},\,\,\,\,\,\,\,\,\,\,\hat{\theta}_4=\min(Y_1,\,Y_2,\,Y_3),\,\,\,\,\,\,\,\,\,\,\hat{\theta}_5=\overline{Y}\]

  1. ¿Cuáles son insesgados?
  2. Entre los estimadores insesgados, ¿cuál tiene la varianza más pequeña?

Exercise 3.18 Sean \(Y_1,\,Y_2,\cdots,Y_n\) una muestra con \(E(Y_i)=\mu\) y \(V(Y_i)=\sigma^2\). Demuestre que \[S'^2=\frac{1}{n}\sum_{i=1}^n(Y_i-\overline{Y})^2\] es un estimador sesgado para \(\sigma^2\) y que \[S^2=\frac{1}{n-1}\sum_{i=1}^n(Y_i-\overline{Y})^2\] es un estimador insesgado para \(\sigma^2\).

Exercise 3.19 Sea \(\{X_n\}\) para \(n\leq 1\) , sucesión de variables aleatorias independientes e idénticamente distribuidas definidas en el espacio de probabilidad \((\Omega, \mathcal{A}, P)\) con función de distribución común \(F\) . Se denota por \(F_n\) la función de distribución empírica obtenida de las \(n\) primeras variables aleatorias \(X_1 ,\ldots, X_n\) . Sea \(x\in\mathbb{R}\). Verifique lo siguiente:

  1. \(P(nF_n(x)=j)=P(F_n(x)=\frac{j}{n})=\binom{n}{j}(F(x))^j(1-F(x))^{n-j}\), \(j=1,\ldots,n\)\

  2. \(E(F_n(x))=F(x)\); \(Var(F_n(x))=\frac{1}{n}F(x)(1-F(x))\).

  3. \(\lim _{n\to \infty }F_{n}(x)=F(x)\)

  4. \(\lim _{n\to \infty }\frac{F_{n}(x)-F(x)}{\sqrt{\frac{F(x)(1-F(x))}{n}}}=Z\), donde \(Z\) es una variable aleatoria con distribución normal estándar y la convergencia es convergencia en distribución.

Exercise 3.20 Se toma una muestra de tamaño tres de una población con distribución de Poisson de parámetro \(\theta\) cuyos los resultados son, \(x_1 = 2\), \(x_2 = 0\), \(x_3 = 5\). Determinar la estimación máximoverosimil de \(\theta\).

3.14 Referencias

  • Gómez, Guadalupe, & Delicado, Pedro (2006). Curso de Inferencia y Decisión. Departament d’Estadística i Investigació Operativa, Universitat Politècnica de Catalunya.

  • Wackerly, D. D., Mendenhall, W., & Scheaffer, R. L. (2008). Estadística matemática con aplicaciones (7ª ed.). Cengage Learning.

  • Roussas, G. G. (1997). A Course in Mathematical Statistics (2nd ed.). Academic Press.

  • Kalbfleisch, J. G. Probability and Statistical Inference. Springer-Verlag, 1985.