3 Estimación puntual

Lurking behind a parametric model is usually some notion of causality. (…) Abandoning mathematical models comes close to abandoning the historic scientific goal of understanding nature.

Bradley Efron

3.1 Modelos estadísticos

Sea \(X_1,\ldots,X_n\) un conjunto de \(n\) variables aleatorias independientes e idénticamente distribuidas con distribución \(F\), y \(x_1,\ldots,x_n\) los \(n\) valores observados o realizaciones correspondientes. La distribución \(F\) no es conocida pero se supone que pertenece a una familia de posibles distribuciones. Formular un modelo estadístico es simplemente especificar cuál es la familia de posibles distribuciones. Un modelo es paramétrico si cada distribución de la familia es totalmente conocida salvo por el valor de un parámetro \(\theta\in\mathbb{R}^d\). Es decir, la familia de posibles distribuciones es \(\{F_\theta:\, \theta\in\Theta \subset \mathbb{R}^d\}\). El conjunto \(\Theta\) de posibles valores del parámetro se llama espacio paramétrico. En los modelos que consideramos en estas notas, habitualmente \(d=1\) (por ejemplo, la familia de distribuciones de Poisson) o \(d=2\) (por ejemplo, la familia de distribuciones gaussianas). Siempre vamos a suponer que se da la siguiente condición de identificabilidad: si \(\theta\neq \theta'\), entonces \(F_\theta\neq F_{\theta'}\). Si esta condición no se cumpliera sería imposible poder distinguir a partir de la muestra entre dos valores \(\theta\) y \(\theta'\) para los que la distribución de probabilidad es la misma.

Pueden considerarse también modelos no paramétricos. Por ejemplo, \(\{F:\, F\ \mbox{tiene función de densidad}\ f\}\). En este caso, el espacio paramétrico es el conjunto de todas las funciones de densidad, que no tiene dimensión finita. Una posible clase de estimadores en este contexto es la de estimadores de núcleo que definimos en la sección 1.5.4.

Cuando se trabaja con un modelo paramétrico, el objetivo general de la inferencia estadística es extraer información sobre el parámetro \(\theta\) a partir de las observaciones \(x_1,\ldots,x_n\). En este capítulo consideramos métodos generales para obtener estimadores (aproximaciones al verdadero valor de \(\theta\)) a partir de los datos muestrales. Estudiaremos tres métodos bastante diferentes:

El método de momentos
El método de máxima verosimilitud
Los estimadores bayesianos

3.2 El método de momentos

Conceptualmente este es el método más sencillo para construir estimadores, aunque en general las propiedades de los estimadores obtenidos no suelen ser óptimas. No obstante, en algunos casos particulares importantes, este método proporciona los mismos resultados que otros procedimientos más sofisticados.

Sean \(X_1,\ldots,X_n\) v.a.i.i.d. con distribución determinada por \(f(\cdot;\theta)\), donde \(\theta=(\theta_1,\ldots,\theta_d)\) es un parámetro \(d\)-dimensional. Consideramos los momentos de la distribución de \(X\), \(\alpha_k(\theta):=\mbox{E}_\theta(X^k)\), \(k=1,\ldots,d\), que son funciones de \(\theta\). Un procedimiento natural para obtener un estimador de \(\theta\) es aproximar los momentos anteriores por los promedios análogos de los datos muestrales. El problema entonces se reduce a resolver en \(\theta_1,\ldots,\theta_d\) el sistema de ecuaciones \[\hat m_1=\alpha_1(\theta) \; , \; \ldots \; , \; \hat m_d=\alpha_d(\theta),\] donde \(\hat m_k\) es el momento muestral de orden \(k\), es decir, \(\hat m_k=\frac{\sum_{i=1}^nX_i^k}{n}\).

La idea es estimar el parámetro como aquel valor de \(\theta\) que hace que los momentos poblacionales (tantos como componentes tenga \(\theta\)) coincidan con los correspondientes momentos muestrales. En general, si \(\theta\) es el verdadero valor del parámetro, NO sucederá que \(\hat m_k=\alpha_k(\theta)\) (de hecho, \(\hat m_k\) es aleatorio) pero es de esperar que \(\hat m_k\approx \alpha_k(\theta)\) para tamaños muestrales grandes.

La principal ventaja del método de los momentos es que proporciona estimadores sencillos que además suelen ser consistentes ya que, por las leyes de los grandes números, los momentos muestrales convergen a los poblacionales.

Ejemplos

Distribución uniforme en \((0,\theta)\). Supongamos que \(X_1,\ldots,X_n\) es una muestra de v.a.i.i.d. de una distribución uniforme en \((0,\theta)\). Entonces, \(\mu = \theta/2\). El estimador de momentos debe verificar \(\bar{x} = \hat{\theta}/2\), es decir, \(\hat{\theta}=2\bar{x}\).
Un caso particular de la distribución beta. Supongamos que \(X_1,\ldots,X_n\) es una muestra de v.a.i.i.d. de una distribución con densidad \(f(x;\theta) = (\theta + 1)x^\theta\), con \(x\in [0,1]\), \(\theta > -1\).
- Representa gráficamente la densidad para \(\theta=1, 2, 3\). A medida que \(\theta\) aumenta, ¿cómo es de esperar que vayan variando los datos?
- ¿Cuál es el estimador de momentos de \(\theta\)? Compara el resultado con las conclusiones del apartado anterior.

3.3 El método de máxima verosimilitud

Este método da estimadores con ciertas propiedades de optimalidad, por lo que es el método de elección dentro de la inferencia clásica.

3.3.1 Estimadores de máxima verosimilitud

Las realizaciones muestrales \(x_1,\ldots,x_n\) son más o menos probables en función de lo que valga el parámetro \(\theta\). Por ejemplo, si sabemos que los datos vienen de una distribución \(\mbox{N}(\theta, 1)\), entonces si fuese \(\theta = 4\) la probabilidad de que la mitad de los datos sea mayor que 4 es 1/2. Sin embargo, si fuese \(\theta = 0\), es prácticamente imposible que algún dato sea mayor que 4. Por lo tanto, si nos encontráramos con una muestra para la que varias observaciones son mayores que 4 y tenemos que apostar entre los valores \(\theta = 0\) y \(\theta = 4\), ¿cuál sería la mejor opción?

El método de máxima verosimilitud consiste en estimar el parámetro \(\theta\) mediante el valor que hace más verosímiles los datos que realmente hemos observado. Veamos un ejemplo adicional para insistir en la idea antes de dar una definición más formal.

Ejemplo

En una urna cerrada hay 4 bolas, \(\theta\) de ellas son blancas y \(4-\theta\) son negras (\(\theta\) desconocido). Se llevan a cabo dos extracciones de bolas con reemplazamiento. ¿Cuál será el estimador de máxima verosimilitud si una de las bolas extraídas es blanca y la otra es negra? Dicho más coloquialmente, si tuviéramos que apostar por el número de bolas blancas, ¿por qué valor apostaríamos dados los resultados obtenidos? La respuesta se deduce de la siguiente tabla, en la que aparecen las probabilidades de observar una bola blanca y otra negra en función de \(\theta\):

\(\theta\)	\(\mbox{P}\{x_1 = B,x_2=N\}\)
0	0
1	3/16
2	4/16
3	3/16
4	0

Está claro que nadie apostaría por \(\theta=0\) o \(\theta=4\) ya que para estos dos valores es imposible obtener una bola negra y otra blanca. Observando la tabla, vemos que si \(\theta=2\) la probabilidad de obtener una bola negra y otra blanca es máxima y por lo tanto \(\hat{\theta}=2\) es el estimador de máxima verosimilitud. Supongamos que se extrae una tercera bola y resulta ser negra. ¿Cuál es el estimador de máxima verosimilitud en este caso?

La función de verosimilitud

Sea \(x_1,\ldots,x_n\) una realización de una muestra \(X_1,\ldots,X_n\) con función de densidad o de probabilidad conjunta \(g(x_1,\ldots,x_n;\theta)\) en \((x_1,\ldots,x_n)\), donde \(\theta\in\Theta\subset \mathbb{R}^d\). Para esta muestra, la función de verosimilitud \(L:\Theta \to \mathbb{R}\) se define como \[L(\theta) = L(\theta; x_1,\ldots,x_n) = g(x_1,\ldots,x_n;\theta).\]

Si \(X_1,\ldots,X_n\) son v.a.i.i.d. con densidad o probabilidad \(f(x;\theta)\), entonces \(L(\theta) =\prod_{i=1}^n f(x_i;\theta)\)

Hay una función de verosimilitud diferente para cada muestra, pero lo que interesa es cómo varía la verosimilitud al variar \(\theta\) puesto que la muestra ya ha sido observada y, por lo tanto, está fija. A partir de \(L(\theta)\) ya podemos definir el estimador de máxima verosimilitud.

Sea \(x_1,\ldots,x_n\) una muestra y sea \(L(\theta)\) la correspondiente función de verosimilitud. Un estimador de máxima verosimilitud (EMV) de \(\theta\) es un valor \(\hat\theta\in\Theta\) tal que \[L(\hat\theta)= \sup_{\theta\in\Theta}L(\theta).\]

En el caso en el que \(\Theta\) es un conjunto abierto y \(L(\theta)\) es derivable y cóncava en \(\Theta\), para calcular el EMV basta derivar e igualar a cero. Sin embargo, como es más fácil derivar sumas que productos, en vez de maximizar directamente \(L(\theta)\) resulta más conveniente maximizar \[\ell(\theta; x)=\ell(\theta)=\log L(\theta)=\sum_{i=1}^n\log f(x_i,\theta).\] Como el logaritmo es creciente, tanto \(L(\theta)\) como \(\ell(\theta)\) alcanzan el máximo en el mismo punto.

3.3.2 Ejemplos y cuestiones

Distribución de Poisson

El número de llamadas por hora que se reciben en una centralita telefónica por hora sigue una distribución de Poisson de parámetro \(\lambda\). En cuatro horas diferentes se han observado de forma independiente los siguientes números de llamadas: \(x_1=2\), \(x_2=3\), \(x_3=5\), \(x_4=5\). Calcula el EMV de \(\lambda\). Da una expresión general del EMV si se dispone de una muestra \(x_1,\ldots,x_n\).

Calculamos la función de verosimilitud. Usando la independencia de las observaciones: \[L(\lambda) = e^{-\lambda}\frac{\lambda^2}{2!}e^{-\lambda}\frac{\lambda^3}{3!}e^{-\lambda}\frac{\lambda^5}{5!}e^{-\lambda}\frac{\lambda^5}{5!} = e^{-4\lambda}\frac{\lambda^{15}}{2!3!5!5!}.\] Tomamos logaritmos: \[\ell(\lambda) = -4\lambda + 15\log\lambda - c,\] donde \(c=\log(2!3!5!5!)\). Representamos gráficamente esta función (salvo la constante, que es irrelevante para el punto en el que se alcanza el valor máximo) en la figura 3.1.

Figura 3.1: Logaritmo de la verosimilitud (salvo constantes) para una muestra de una distribución de Poisson y valor en el que se alcanza el máximo.

Derivando e igualando a cero tenemos que \(-4+15/\hat{\lambda}=0\), de donde \(\hat{\lambda} = 15/4\). Este es el valor en el que se ha situado la línea vertical en la figura 3.1.

Es fácil demostrar siguiendo los mismos pasos (ejercicio) que \(\hat{\lambda} = \bar{x}\) es el EMV de \(\lambda\) en una distribución de Poisson para una muestra genérica \(x_1,\ldots,x_n\). Obsérvese que coincide con el estimador de momentos (¿por qué?).

Distribución exponencial

Sea \(X_1,\ldots,X_n\) una muestra de v.a.i.i.d. de una distribución exponencial de parámetro \(\theta\) (densidad: \(f(x)=\theta e^{-\theta x}\), si \(x\geq 0\)). Comprueba que el EMV de \(\theta\) es \(\hat\theta = 1/\bar{x}\). ¿Coincide con el estimador obtenido con el método de momentos?

Distribución normal

Sea \(X_1,\ldots,X_n\) una muestra de v.a.i.i.d. de una distribución normal de media \(\mu\) y varianza \(\sigma^2\). No es difícil comprobar que el EMV de \(\mu\) es \(\hat\mu = \bar{X}\) y el de \(\sigma^2\) es \(\hat\sigma^2 = n^{-1}\sum_{i=1}^n (X_i-\bar{X})^2\) (ejercicio). En este caso el EMV de la varianza no es insesgado (ya que, como sabemos, el estimador insesgado requiere dividir por \(n-1\) en lugar de \(n\)).

Observaciones truncadas

La idea de maximizar la verosimilitud se puede extender a situaciones más generales. Se puede aplicar por ejemplo a situaciones mixtas en las que algunos de los factores que componen la verosimilitud son funciones de densidad (como en el ejemplo de la distribución exponencial) y otros son probabilidades (como en el caso de la distribución de Poisson). Depende del tipo de datos y de la información disponible. En la siguiente situación hay que considerar de forma natural factores de ambos tipos.

El tiempo de vida de los ratones con cierta enfermedad sometidos a un tratamiento es una variable aleatoria con distribución exponencial de parámetro \(\theta\). Se lleva a cabo un experimento con \(n\) ratones, se observa el tiempo de vida de \(m\) de ellos (\(x_1,\ldots,x_m\)) pero se interrumpe el experimento una vez transcurrido un tiempo \(T\), de manera que de los \(n-m\) restantes solo se sabe que su tiempo de vida es superior a \(T\). Se dice que estas observaciones están truncadas ya que solo tenemos una cota inferior del tiempo de vida de los ratones cuya supervivencia es mayor que \(T\). Se supone que todos los tiempos son independientes. En este caso la verosimilitud se descompone en dos partes: la correspondientes a los datos no truncados y la correspondiente a los datos truncados. Si \(f(x;\theta)=\theta e^{-\theta x}\), \(x>0\), es la función de densidad exponencial y \(F(x;\theta)=1-e^{-\theta x}\) es la función de distribución, la función de verosimilitud es \[L(\theta) = (1-F(T;\theta))^{n-m}\prod_{i=1}^m f(x_i;\theta)=e^{-(n-m)\theta T} \theta^{m} e^{-\theta (x_1+\cdots + x_m)}.\] Hemos usado que la probabilidad de que un ratón viva un tiempo mayor que \(T\) es \(1-F(T;\theta)=e^{-\theta T}\). Tomando logaritmos, \[\ell(\theta)=\log L(\theta) = -(n-m)\theta T+ m\log\theta - \theta (x_1+\cdots+x_m).\] Derivando e igualando a cero obtenemos el EMV: \[\hat\theta = \frac{m}{x_1+\cdots+x_m + (n-m)T}.\] Si comparamos con el caso en que no hay truncamiento, se han sustituido en el denominador los \(n-m\) tiempos desconocidos de los datos truncados por su cota inferior, \(T\). La infraestimación de la suma de tiempos de vida se compensa usando \(m\) en el numerador en lugar de \(n\). Así pues, el resultado es razonable, y se ha obtenido de forma automática al aplicar el principio general de maximizar la verosimilitud.

Distribución uniforme en \((0,\theta)\)

Supongamos que \(X_1,\ldots,X_n\) es una muestra de v.a.i.i.d. de una distribución uniforme en \((0,\theta)\). Veamos cuál es el EMV de \(\theta\). En este caso, la función de verosimilitud es \[L(\theta) = \left\{\begin{array}{cc} \frac{1}{\theta^n}, & \mbox{si} \ \theta\geq X_{(n)}, \\ 0, & \mbox{si} \ \theta< X_{(n)}. \end{array} \right.\] Esta función de verosimilitud no es derivable (ni siquiera es continua). Si la representamos gráficamente es claro que el EMV es \(\hat{\theta} = X_{(n)}\), el máximo de las observaciones. En este ejemplo, el método de momentos y el de máxima verosimilitud dan estimadores muy diferentes. ¿Cuál crees que es mejor?

Salvo en modelos muy sencillos el EMV no se puede calcular por el método estándar de derivar e igualar a cero. Lo habitual es tener que usar algún algoritmo de optimización. Cuidado especial requieren los casos (como el último ejemplo) en los que el soporte de la distribución es función del parámetro.

3.3.3 El algoritmo EM

Los estimadores de máxima verosimilitud suelen tener excelentes propiedades, pero a veces involucran un problema de optimización difícil de resolver. El algoritmo EM (Expectation-Maximization) es un método para calcular el estimador de máxima verosimilitud particularmente fácil de implementar en casos en los que falta cierta información tal que el problema sería mucho más sencillo si se dispusiera de ella. Esto ocurre, por ejemplo, cuando hay datos faltantes, los datos están censurados o corresponden a una mezcla de distribuciones, sin que se sepa a qué distribución corresponde cada observación. El algoritmo fue propuesto por Dempster, Laird y Rubin (1977) en un famoso artículo que acumula decenas de miles de citas.

Planteamiento y algoritmo

Disponemos de una muestra \(x=(x_1,\ldots,x_n)\) procedente de una distribución que depende de un parámetro \(\theta\) y queremos calcular el estimador de máxima verosimilitud de \(\theta\).

Llamamos \(y\) a otros datos que proporcionarían una información adicional que convertiría el problema en otro mucho más fácil. Normalmente estos datos son en realidad inaccesibles, pero eso no importa en la aplicación del método.

En esta situación podemos plantearnos dos problemas:

El problema que realmente tenemos que resolver, que requiere maximizar el logaritmo de la verosimilitud para \(x\), \(\ell (\theta; x)\).
El problema simplificado, que requiere maximizar el logaritmo de la verosimilitud completa para \((x,y)\) que denotamos \(\ell_c(\theta; x,y)\).

Aunque \(y\) sea en realidad inobservable, si maximizar \(\ell_c\) es más fácil que maximizar \(\ell\) puede ser útil aproximar \(\ell_c(\theta;X,Y)\) calculando su esperanza condicionada dado \(X=x\), cuando el valor del parámetro es el más reciente obtenido en la aplicación del algoritmo.

La aplicación de la idea anterior lleva al algoritmo EM:

Algoritmo EM

Iniciamos el método con algún valor sensato \(\theta_0\). Sea \(\theta_k\) el valor del parámetro tras \(k\) iteraciones. En cada cada iteración hay que dar los dos pasos siguientes:

PASO E (Esperanza) Se sustituye la log-verosimilitud completa por su esperanza condicionada a \(x\), suponiendo que \(\theta_k\) fuese el verdadero valor del parámetro:

\[Q(\theta;\theta_k,x):= \mbox{E}_{Y|x;\theta_k}[\ell_c(\theta;x,Y)].\]

PASO M (Maximización) Se maximiza en \(\theta\) la función \(Q(\theta;\theta_k,x)\), lo que da lugar al siguiente valor del parámetro, \(\theta_{k+1}\).

Se reiteran los pasos E y M hasta que se cumple algún criterio de convergencia.

Bajo condiciones de regularidad (hay muchos artículos que analizan el comportamiento de este método) puede demostrarse que en cada paso se mejora el objetivo del problema que realmente queremos resolver, es decir, \(\ell(\theta_{k+1}; x)\geq \ell(\theta_k; x)\) y que la sucesión de valores generada converge a \(\hat\theta\), el estimador de máxima verosimilitud correspondiente a la muestra \(x\) que realmente tenemos.

Ejemplo

Supongamos que las observaciones \(x_1,\ldots,x_n\) son independientes y vienen con probabilidad \(1-p\) de una distribución con densidad \(f_0\) y con probabilidad \(p\) de otra distribución \(f_1\). Tanto \(f_0\) como \(f_1\) son conocidas y el objetivo es estimar la proporción de observaciones \(p\) que vienen de \(f_1\). La función de densidad de los datos es \(f(x) = (1-p)f_0(x) + pf_1(x)\) (se puede comprobar con la fórmula de la probabilidad total). Por lo tanto, la log-verosimilitud que tenemos que maximizar es \[\ell(p; x) = \sum_{i=1}^n \log [(1-p)f_0(x_i) + p f_1(x_i)].\] No es posible encontrar explícitamente el valor de \(p\) que maximiza esta función por lo que hay que recurrir a algún método numérico.

Si alguien nos dijera a qué distribución pertenece cada observación el problema sería inmediato de resolver. Escribimos \(y_i=1\) si la observación \(x_i\) viene de \(f_1\) e \(y_i=0\) si viene de \(f_0\). La verosimilitud completa correspondiente a \((x,y)\) es \[ L_c(p;x,y) := \prod_{i=1}^n [(1-p)f_0(x_i)]^{1-y_i}[p f_1(x_i)]^{y_i}, \] y tomando logaritmos, \[ \ell_c(p; x,y) = \sum_{i=1}^n [(1-y_i)\log(1-p)+y_i\log p ] + C, \] donde \(C\) no depende del parámetro \(p\), por lo que es irrelevante para nuestros propósitos. Derivando respecto a \(p\) es fácil ver que el óptimo se alcanza en el valor \(n^{-1}\sum_{i=1}^n y_i\). Como no podía ser de otra manera, si conocemos de dónde viene cada observación estimamos \(p\) como la proporción de valores muestrales que vienen de \(f_1\).

Nosotros no sabemos de qué distribución viene cada \(x_i\) pero si aplicamos la fórmula de Bayes y usamos el valor \(p_k\) del parámetro,

\[ y_{i,k} :=\mbox{E}(Y_i|x_i;p_k) = \frac{p_kf_1(x_i)}{p_kf_1(x_i) + (1-p_k)f_0(x_i)}. \tag{3.1}\]

Como consecuencia, los dos pasos del algoritmo se reducen a:

PASO E:

Salvo sumandos que no dependen de \(p\),

\[Q(p;p_k) = \sum_{i=1}^n [(1-y_{i,k})\log(1-p) + y_{i,k}\log p],\] donde \(y_{i,k}\) se calcula mediante la ecuación (3.1).

Paso M:

Maximizamos \(Q(p;p_k)\). Para ello derivamos e igualamos a cero y obtenemos \[ p_{k+1} = \frac{1}{n}\sum_{i=1}^n y_{i,k} \] Un ejemplo numérico con datos simulados se muestra a continuación. Primero simulamos \(n=100\) datos que pueden venir de una normal estándar (con probabilidad \(p=0.8\)), o de una normal de media 4 y varianza 1 (con probabilidad \(1-p=0.2\)):

# Generamos datos de una N(0,1) o de una N(4,1)
set.seed(156)
n <- 100
p <- 0.8  # objetivo a estimar
media <- 4

y <- rbinom(n, 1, p)  
x <- rnorm(n)
n0 <- sum(y==0)
x[y==0] <- rnorm(n0, mean = media)

# Representa histograma de los datos

library(ggplot2)
ggplot(data.frame(x)) +
  geom_histogram(aes(x = x, y = after_stat(density)), col = 'black', fill = 'olivedrab4', bins = 20)

En el histograma se observa claramente la bimodalidad correspondiente a las dos posibles distribuciones de las que pueden venir los datos. Ahora fingimos que no conocemos que \(p=0.8\) y aplicamos el algoritmo EM a la muestra x (no usamos para nada y). Usamos el valor inicial \(p_0=1/2\).

# Iniciamos el algoritmo
p_k <- 0.5
tol <- 0.1

# Pasos E-M
while(tol > 0.001){
  p_0 <- p_k
  y_k <- p_0*dnorm(x)/(p_k*dnorm(x) + (1-p_0)*dnorm(x, mean = media))  # Paso E
  p_k <- mean(y_k) # Paso M
  tol <- abs(p_k-p_0)
}
p_k
#> [1] 0.8179144

Se obtiene el valor 0.818, cercano al valor verdadero. Hemos supuesto las dos distribuciones conocidas. Una situación de más interés y más complicada sería no suponer las medias conocidas y tratar de estimar por máxima verosimilitud tanto \(p\) como las dos medias.

3.3.4 Información de Fisher y cota de Fréchet-Cramer-Rao (FCR)

En esta sección y la siguiente vamos a profundizar en las propiedades estadísticas de los EMV. La muestra está formada por v.a.i.i.d. \(X_1,\ldots,X_n\) con distribución en la familia \(\{F_\theta:\, \theta\in\Theta \subset \mathbb{R}\}\). Denotamos por \(f(\cdot;\theta)\) a la función de densidad o de probabilidad de \(F_\theta\). Para simplificar utilizaremos -salvo que se diga lo contrario y si no hay riesgo de ambigüedad- la notación \(X\equiv (X_1,\ldots,X_n)\) y \(x\equiv (x_1,\ldots,x_n)\). Recordamos también la notación \(g(x;\theta)=g(x_1,\ldots,x_n;\theta)\), \(\ell(\theta;X)=\log L(\theta;X) = \log g(X; \theta)\). En el caso i.i.d., \(\ell(\theta;X) = \sum_{i=1}^n \log f(X_i;\theta)\).

Los resultados que vamos a estudiar requieren ciertas condiciones de regularidad para cumplirse. Aunque estas condiciones parecen complicadas, en realidad las cumplen prácticamente todas las distribuciones de probabilidad habituales, con la excepción de aquellas cuyo soporte depende del parámetro.

Condición R1: \(\Theta\) es abierto.
Condición R2: El soporte \(S = \{t:\, f(t;\theta)>0\}\) no depende de \(\theta\).
Condición R3: Para todo \(t\in S\), la función \(f(t;\theta)\) es derivable una, dos o tres veces respecto a \(\theta\) con derivada continua (según las derivadas que aparezcan en cada resultado).

Conviene insistir en que todas las derivadas que aparecen son siempre respecto al parámetro (con \(X\) fijo) y las esperanzas y varianzas se calculan respecto a \(X\) para cada valor del parámetro.

Muchas propiedades de los EMV dependen de la llamada cantidad de información de Fisher contenida en la muestra \(X\). Esta importante función se define de la siguiente forma: \[I_n(\theta) = \mbox{Var}[\ell'(\theta;X)].\] Dado que las observaciones son i.i.d., \(\ell'(\theta;X)=\sum_{i=1}^n \ell'(\theta;X_i)\) y por lo tanto \(I_n(\theta) = nI(\theta)\), donde \(I(\theta)\), sin subíndice, denota la información de Fisher de cada observación.

Condición R4: Para todo \(\theta\in \Theta\), \(0<I(\theta)<\infty\).
Condición R5: Para todo \(\theta\in \Theta\), \(\mbox{E}[\ell'(\theta;X)] = 0\).
Condición R6: Para todo \(\theta\in \Theta\), \(I_n(\theta) = \mbox{E}[-\ell''(\theta;X)]\).

Observación

Para que se cumplan las condiciones R5 y R6 lo único que hace falta es poder intercambiar las derivadas con las integrales en el siguiente razonamiento (lo hacemos solo en el caso continuo, el discreto es totalmente análogo).

En primer lugar veamos que si se puede hacer el intercambio de la derivada con la integral, se cumple R5. Observemos que \[\frac{\partial}{\partial\theta}\log g(x;\theta) = \frac{1}{g(x;\theta)}\frac{\partial}{\partial\theta} g(x; \theta),\] y por lo tanto \[\frac{\partial}{\partial\theta} g(x; \theta) = \left[\frac{\partial}{\partial\theta}\log g(x;\theta)\right] g(x;\theta). \tag{3.2}\] Como consecuencia, \[\mbox{E}(\ell'(\theta;X)) = \int\left[\frac{\partial}{\partial\theta}\log g(x;\theta)\right]g(x;\theta)dx = \int \frac{\partial}{\partial\theta} g(x; \theta) dx = \frac{\partial}{\partial\theta}\int g(x; \theta) dx = 0.\] ¿A qué se debe la última igualdad?

Veamos ahora que R6 también se sigue de un cambio de orden entre las integrales y las derivadas. Si derivamos de nuevo respecto a \(\theta\) en la ecuación \[0 = \mbox{E}(\ell'(\theta;X)) = \int\left[\frac{\partial}{\partial\theta}\log g(x;\theta)\right]g(x;\theta)dx\] tenemos que \[0 = \int\left[\frac{\partial^2}{\partial\theta^2}\log g(x;\theta)\right]g(x;\theta)dx +\int\left[\frac{\partial}{\partial\theta}\log g(x;\theta)\right]\left[\frac{\partial}{\partial\theta} g(x; \theta)\right]dx,\] y sustituyendo el valor del último corchete por la expresión en (3.2) se sigue \[0 = \int\left[\frac{\partial^2}{\partial\theta^2}\log g(x;\theta)\right]g(x;\theta)dx+\int\left[\frac{\partial}{\partial\theta}\log g(x;\theta)\right]^2 g(x; \theta)dx.\] Cambiando la notación y aplicando R5, la expresión anterior es equivalente a: \[0=\mbox{E}[\ell''(\theta, X)] + \mbox{E}[\ell'(\theta; X)^2] = \mbox{E}[\ell''(\theta, X)] + I_n(\theta),\] lo que demuestra R6.

Bajo estas condiciones, la información de Fisher se puede calcular como la varianza de \(\ell'(\theta; X)\) o como la esperanza de \(-\ell''(\theta; X)\). Este segundo método suele ser más sencillo.

Ejemplo

Veamos un ejemplo suponiendo que tenemos una muestra de \(n\) observaciones i.i.d. de una distribución de Poisson de parámetro \(\lambda\). En este caso tenemos, para cada observación, \(\ell(\lambda, X) = -\lambda + X\log\lambda - \log X!\). Como consecuencia, \(\ell'(\lambda;X)=-1+X/\lambda\) y \(-\ell''(\lambda;X) = X/\lambda^2\). Recordando que \(\mbox{E}(X)=\lambda\) para una v.a. con distribución de Poisson, tenemos que \[I_n(\lambda) = nI(\lambda) = n\mbox{E}[-\ell''(\lambda;X)]= \frac{n}{\lambda}.\]

Una primera aplicación de la información de Fisher es que interviene en una cota inferior de la varianza de cualquier estimador insesgado. De esta forma, si se encuentra un estimador insesgado tal que su varianza alcanza la cota, en términos de sesgo-varianza es el estimador óptimo.

Teorema 3.1 (Cota de Fréchet-Cramer-Rao) Sea \(T(X)=T(X_1,\ldots,X_n)\) un estimador insesgado de \(h(\theta)\), donde \(h\) es una función derivable. Entonces, bajo las condiciones de regularidad R1-R5 y si \(h'(\theta)=\mbox{E}[T(X)\ell'(\theta;X)]\), se cumple \[\mbox{Var}[T(X)] \geq \frac{h'(\theta)^2}{I_n(\theta)}.\]

Prueba. La hipótesis sobre \(h'(\theta)\) lo único que requiere es de nuevo poder intercambiar la derivada con la integral. Como \(T\) es insesgado, \[h(\theta) = \mbox{E}[T(X)] = \int T(x) g(x; \theta) dx.\] Derivando respecto al parámetro y aplicando (3.2), \[h'(\theta) = \int T(x) \left[\frac{\partial}{\partial\theta}\log g(X;\theta)\right] g(X;\theta)dx = \mbox{E}[T(X)\ell'(\theta;X)] = \mbox{Cov}[T(X)\ell'(\theta;X)].\] Elevando al cuadrado y usando que el coeficiente de correlación en valor absoluto siempre es menor o igual a uno, \[h'(\theta)^2 = \mbox{Cov}[T(X)\ell'(\theta;X)]^2 \leq \mbox{Var}[T(X)]\mbox{Var}[\ell'(\theta;X)] = \mbox{Var}[T(X)]I_n(\theta).\]

Cuanto mayor es la información de Fisher, la cota deja más margen de mejora para encontrar un estimador insesgado con menor varianza. Si un estimador insesgado alcanza la cota FCR se dice que es eficiente.

Si las observaciones son i.i.d. y \(h(\theta) =\theta\), entonces \[\mbox{Var}[T(X)] \geq \frac{1}{nI(\theta)},\] donde \(I(\theta)\) es la información de Fisher de cada una de las observaciones.

Ejemplo: distribución de Poisson (continuación)

Cualquier estimador insesgado \(T(X)\) de \(\lambda\) a partir de v.a.i.i.d. de una distribución de Poisson debe verificar: \[\mbox{Var}[T(X)] \geq \frac{\lambda}{n},\] ya que hemos visto que \(I_n(\lambda)= n/\lambda\).

Sabemos que la media \(\bar{X}\) es un estimador insesgado para \(\lambda\) y, además, \(\mbox{Var}(\bar{X}) = \sigma^2/n = \lambda/n\), por lo que alcanza la cota. Esto significa que la media muestral es óptimo en términos de sesgo y varianza para estimar el parámetro de una distribución de Poisson con muestras de observaciones i.i.d..

Cuestión: distribución de Bernoulli

Dada una muestra de v.a.i.i.d. de una distribución de Bernoulli de parámetro \(p\), calcula la cota de FCR que satisface la varianza de cualquier estimador insesgado obtenido a partir de dicha muestra.

3.3.5 Comportamiento asintótico del EMV

En este apartado vamos a demostrar que bajo las condiciones de regularidad, los estimadores de máxima verosimilitud son asintóticamente normales. La información de Fisher tiene un papel importante en la distribución asintótica de los EMV.

Teorema 3.2 (Normalidad asintótica del EMV en el caso i.i.d.) Sea \(\hat\theta\) el EMV de \(\theta\) basado en una muestra de v.a.i.i.d. \(X_1,\ldots,X_n\). Suponemos que se cumplen las condiciones de regularidad R1-R6 y además las dos siguientes:

Condición R7: Para todo \(\theta\) y \(\delta>0\), \(|\ell'''(t,x)|\leq M(x)\) para \(|\theta-t|\leq \delta\), donde \(M\) es una función tal que \(\mbox{E}[M(X_i)]<\infty\).
Condición R8: (Consistencia) \(\hat{\theta} \overset{\mbox{p}}{\to} \theta\).

Entonces, \[\sqrt{n}(\hat\theta-\theta) \overset{\mbox{d}}{\to} \mbox{N}\left(0,\frac{1}{I(\theta)}\right),\] donde \(I(\theta)\) es la información de Fisher correspondiente a cada observación.

Prueba. Se basa en el siguiente desarrollo de Taylor: \[0 = \sum_{i=1}^n \ell'(\hat \theta, x_i) = \sum_{i=1}^n \ell'(\theta, x_i) + (\hat\theta-\theta)\sum_{i=1}^n \ell''(\theta, x_i) + \frac{1}{2}(\hat\theta-\theta)^2\sum_{i=1}^n \ell'''(\tilde\theta, x_i),\] donde \(\tilde\theta\) es un valor situado entre \(\theta\) y \(\hat{\theta}\).

Despejando, tenemos \[\sqrt{n}(\hat\theta-\theta) = - \frac{\sqrt{n}\frac{1}{n}\sum_{i=1}^n \ell'(\theta, x_i)}{\frac{1}{n}\sum_{i=1}^n \ell''(\theta, x_i) + \frac{1}{2n}(\hat\theta-\theta)\sum_{i=1}^n \ell'''(\tilde\theta, x_i)} = \frac{ \sqrt{n}\bar{U}_n}{\bar{H}_n + R_n},\] donde \(U_i= \ell'(\theta;X_i)\), \(H_i=-\ell''(\theta;X_i)\) y \(R_n = -(2n)^{-1}(\hat\theta-\theta)\sum_{i=1}^n \ell'''(\tilde\theta, x_i)\).

Supongamos que podemos probar que \(R_n\to_P 0\) (¿por qué es razonable esperarlo?). Vamos a estudiar los otros dos términos.

Como \(\mbox{E}(U_i) = 0\) (R5) y \(0<\mbox{Var}(U_i) = I(\theta)<\infty\) (R4), por el TCL tenemos \(\sqrt{n}\bar{U}\overset{\mbox{d}}{\to} \mbox{N}(0, I(\theta))\). Por otra parte, \(\mbox{E}(H_i) = I(\theta)\) (R6) y aplicando la LDGN, \(\bar{H} \overset{\mbox{p}}{\to} I(\theta)\). Finalmente, usando el lema de Slutsky (proposición 2.2), \[\frac{ \sqrt{n}\bar{U}_n}{\bar{H}_n + R_n} \overset{\mbox{d}}{\to} \mbox{N}\left(0,\frac{1}{I(\theta)}\right).\]

Para terminar la demostración basta probar que \(R_n\overset{\mbox{p}}{\to} 0\). Para ello usaremos las condiciones R7 y R8. Sea \(\epsilon>0\). Para cualquier \(\delta >0\), \[\mbox{P}(|R_n|>\epsilon) = \mbox{P}(|R_n|>\epsilon,|\hat\theta-\theta|>\delta) + \mbox{P}(|R_n|>\epsilon,|\hat\theta-\theta|\leq\delta).\] El primer sumando converge a cero, \[\mbox{P}(|R_n|>\epsilon,|\hat\theta-\theta|>\delta)\leq \mbox{P}(|\hat\theta-\theta|>\delta) \to 0,\] dado que el EMV es débilmente consistente. Para el segundo sumando,la condición R7 implica que si \(|\hat\theta-\theta|\leq\delta\), entonces \[|R_n| \leq \frac{\delta}{2n} \sum_{i=1}^n M(X_i) \Leftrightarrow \frac{1}{n}\sum_{i=1}^n M(X_i)\geq \frac{2|R_n|}{\delta}.\] Por lo tanto, \[\mbox{P}(|R_n|>\epsilon,|\hat\theta-\theta|\leq\delta) \leq \mbox{P}\left(\frac{1}{n}\sum_{i=1}^n M(X_i) > \frac{2\epsilon}{\delta} \right).\] Pero \[\mbox{P}\left(\frac{1}{n}\sum_{i=1}^n M(X_i) > \frac{2\epsilon}{\delta} \right)=\mbox{P}\left(\frac{1}{n}\sum_{i=1}^n M(X_i) - \mbox{E}[M(X_1)]> \frac{2\epsilon}{\delta} - \mbox{E}[M(X_1)]\right)\to 0,\] para \(\delta>0\) suficientemente pequeño, puesto que \(n^{-1} \sum_{i=1}^n M(X_i) \overset{\mbox{p}}{\to} \mbox{E}[M(X_1)]<\infty\).

Observaciones

Para \(n\) suficientemente grande el resultado anterior muestra que \[\hat\theta \cong \mbox{N}\left(\theta,\frac{1}{nI(\theta)}\right).\] Esto significa que para tamaños muestrales grandes podemos esperar que el EMV sea aproximadamente insesgado y que -también aproximadamente- alcance la cota de FCR. Por ello, a veces se dice que el EMV es asintóticamente óptimo o asintóticamente eficiente.
La variabilidad asintótica del EMV es inversamente proporcional a la información de Fisher. Dado que \(I(\theta) = \mbox{E}[-\ell''(\theta,X)]\), la información de Fisher es una medida de la curvatura esperada de \(\ell(\theta;x)\), el logaritmo de la verosimilitud. Esto es bastante intuitivo puesto que si (en términos esperados) la log-verosimilitud tiene poca curvatura y es bastante plana, valores muy diferentes del parámetro dan lugar a modelos de probabilidad similares y, por lo tanto, difíciles de distinguir entre sí. Numéricamente, un pequeño cambio en una función con poca curvatura (el que ocurre si en lugar de una muestra usamos otra) puede hacer que el punto en el que se alcanza su máximo varíe bastante. Sin embargo, si hay mucha curvatura el punto donde se alcanza el máximo es más estable.
Los EMV son consistentes bajo condiciones de regularidad. La consistencia resulta bastante plausible: hemos visto que \(\mbox{E}[\ell'(\theta, X)] = 0\) mientras que, por otra parte, el EMV resuelve la ecuación \(n^{-1}\sum_{i=1}^n \ell'(\theta, X_i) = 0\), que es justo la versión muestral de la ecuación resultante de sustituir la esperanza por el promedio. Es de esperar que para valores grandes de \(n\), ambas ecuaciones se parezcan y por lo tanto sus soluciones también, con lo que \(\hat\theta\approx\theta\). Sin embargo, la convergencia puntual de una sucesión de funciones no implica la convergencia de sus ceros a los de la función límite. Hace falta convergencia uniforme. Además, esta por sí misma no es suficiente ya que también hay que garantizar que la sucesión de soluciones esté eventualmente incluida en un conjunto cerrado y acotado. Por todo ello, demostrar en general la consistencia de los EMV es un resultado bastante técnico y que queda fuera del alcance de estas notas. A menudo es más fácil probar directamente la consistencia para cada caso particular en el que estemos interesados.

3.4 Metodología bayesiana

3.4.1 Estimadores bayesianos

Hay dos interpretaciones clásicas de la probabilidad de un suceso: la primera, que podemos llamar frecuentista, consiste en interpretar la probabilidad de un suceso como el límite de la frecuencia relativa de veces que ocurre este suceso cuando un experimento aleatorio se va repitiendo más y más veces. La segunda interpretación, la bayesiana, consiste en interpretar la probabilidad de un suceso como el grado de creencia subjetiva en que tal suceso ocurra.

Hasta ahora hemos interpretado el parámetro \(\theta\) como una cantidad fija pero desconocida que tratamos de estimar a partir de los datos \(x_1,\ldots,x_n\). Sin embargo, si adoptamos un enfoque bayesiano, tiene sentido describir la incertidumbre sobre el parámetro mediante una distribución de probabilidad definida en el espacio paramétrico \(\Theta\), es decir, tratar a \(\theta\) como si fuera una variable aleatoria en lugar de una constante desconocida. Las probabilidades que asigna esta distribución a los diferentes subconjuntos de \(\Theta\) reflejan la creencia que tenemos de que \(\theta\) pertenezca a cada uno de estos subconjuntos.

El método bayesiano opera entonces de la forma siguiente: en un principio, como hemos dicho, tenemos que establecer una distribución sobre \(\Theta\) que no dependa de los datos del experimento, que sea previa a la observación de la muestra y que refleje, por ejemplo, la opinión de un experto sobre los posibles valores del parámetro. Esta distribución podría estar basada en otros experimentos o muestras recogidas anteriormente. Llamaremos a esta distribución la distribución a priori y denotaremos por \(\pi(\theta)\) su función de densidad o probabilidad. La información sobre \(\theta\) contenida en \(\pi(\theta)\) se puede combinar con el modelo estadístico para los datos \(g(x_1,\ldots,x_n;\theta)\) (que en el contexto bayesiano se suele denotar como \(g(x_1,\ldots,x_n|\theta)\), ya que se interpreta como una distribución condicionada) mediante el teorema de Bayes para calcular la llamada distribución a posteriori \(\pi(\theta|x_1,\ldots,x_n)\), que representa la incertidumbre sobre el parámetro después de observar los datos: \[\pi(\theta \vert x_1,\ldots,x_n)=\frac{g(x_1,\ldots,x_n \vert \theta)\pi(\theta)}{\int_\Theta g(x_1,\ldots,x_n \vert \theta)\pi(\theta)d\theta}.\] Finalmente, se toma como estimador de \(\theta\) alguna medida numérica que resuma la distribución a posteriori, por ejemplo, la esperanza, es decir \[\hat\theta = \mbox{E}(\theta|x_1,\ldots,x_n)=\int_{\Theta}\theta\, \pi(\theta|x_1,\ldots,x_n)\,d\theta.\] Otra posibilidad podría ser usar la moda de la distribución a posteriori: \[\hat\theta =\mbox{arg} \max_{\theta\in\Theta}\pi(\theta|x_1,\ldots,x_n).\]

3.4.2 Ejemplos

Estimación bayesiana del parámetro de una distribución de Poisson con espacio paramétrico finito

El número de llegadas por hora a las emergencias de un hospital sigue una distribución de Poisson de parámetro \(\lambda\). A priori se supone que el valor de \(\lambda\) toma alguno de los siguientes valores con las siguientes probabilidades:

\[\begin{array}{l|rrrrr} \mbox{Valores} & 3 & 3.5 & 4 & 4.5 & 5 \\ \hline \mbox{Probabilidades} & 0.1 & 0.2 & 0.4 & 0.2 & 0.1 \end{array}\]

Se han registrado las visitas durante \(n=10\) periodos de una hora independientes y el número total de visitas ha sido de 31. Veamos cuál es el estimador bayesiano de \(\lambda\).

La función de verosimilitud es \[L(\lambda) = e^{-10\lambda}\lambda^{-31}/c,\] donde \(c\) es una constante que no depende de \(\lambda\). Por tanto, la distribución a posteriori verifica \[\pi(\lambda|x_1,\ldots,x_n) \propto \pi(\lambda) \times e^{-10\lambda}\lambda^{-31}.\] La observación anterior permite obtener la siguiente tabla:

\(\lambda\)	\(\pi(\lambda)\)	Producto	\(\pi(\lambda\|x_1,\ldots,x_n)\)
3	0.10	5.78	0.24
3.50	0.20	9.27	0.39
4.00	0.40	7.84	0.33
4.50	0.20	1.02	0.04
5.00	0.10	0.09	0.00

Si calculamos la media de la distribución a posteriori resulta un valor del parámetro que no coincide con ninguno de los cinco que estamos considerando como posibles a priori. Por ello, en este caso es mejor usar la moda de la distribución a posteriori, esto es, el valor más probable a posteriori (que como se puede comprobar no coincide con el valor más probable a priori):

\[\hat\lambda = \mbox{arg} \max_\lambda \pi(\lambda|x_1,\ldots,x_n) = 3.5.\] No es necesario calcular las constantes de normalización (independientes del parámetro) para resolver el problema.

Estimación de una proporción con distribución a priori beta

Consideramos ahora un problema en el que el espacio paramétrico es un intervalo. Sea \(X_1,\ldots,X_n\) una muestra de v.a.i.i.d. de una distribución \(\mbox{B}(1,\theta)\). Queremos estimar la probabilidad de éxito \(\theta\). Se supone que la distribución a priori de \(\theta\) es una distribución beta de parámetros \(\alpha\) y \(\beta\) adecuados. A los parámetros de la distribución a priori se les suele llamar hiperparámetros. La elección adecuada de los hiperparámetros permite reflejar un conjunto muy amplio de posibles creencias a priori. Por ejemplo, si \(\alpha=\beta > 1\), entonces pensamos que a priori los valores del parámetro alrededor de 1/2 son más probables que cerca de 0 o 1.

Veamos ahora cómo se calcula el estimador bayesiano definido como la media a posteriori \(\hat\theta = \mbox{E}(\theta|x_1,\ldots,x_n)\). Despreciando las constantes que no dependen de \(\theta\), \[\pi(\theta|x_1,\ldots,x_n) \propto \pi(\theta) g(x_1,\ldots,x_n|\theta) = \theta^{\alpha-1}(1-\theta)^{\beta-1} \theta^{\sum_{i=1}^n x_i}(1-\theta)^{n - \sum_{i=1}^n x_i}.\] Por lo tanto (¿Por qué?), \(\theta|x_1,\ldots,x_n \equiv \mbox{Beta}(\alpha + n\bar{x},n+\beta - n\bar{x})\). Como consecuencia, \[\hat\theta = \mbox{E}(\theta|x_1,\ldots,x_n) = \frac{\alpha + n\bar{x}}{n+\beta - n\bar{x}} = \frac{\alpha + \beta}{\alpha +\beta + n}\frac{\alpha}{\alpha + \beta} + \frac{n}{\alpha +\beta +n} \bar{x}= w_n \frac{\alpha}{\alpha + \beta} + (1-w_n)\bar{x}.\]

El estimador final es una mezcla entre lo que pensábamos a priori y la información proporcionada por los datos.

Media de una normal con distribución a priori normal (varianza conocida)

Sea \(X_1,\ldots,X_n\) una muestra de v.a.i.i.d. de una distribución \(\mbox{N}(\mu,\sigma^2)\). Se supone que la distribución a priori de \(\mu\) es también una distribución normal, \(\mbox{N}(\mu_0,\sigma^2_0)\). Vamos a comprobar que la distribución a posteriori es de nuevo normal, \(\mbox{N}(\mu_1,\sigma^2_1)\), donde \[\frac{1}{\sigma^2_1} = \frac{1}{\sigma^2_0} + \frac{n}{\sigma^2}\] y \[\mu_1 = \mu_0\frac{1/\sigma_0^2}{1/\sigma_0^2 + n/\sigma^2} + \bar{x}\frac{n/\sigma^2}{1/\sigma_0^2 + n/\sigma^2}.\] En la situación de este ejemplo tenemos \[\pi(\mu|x_1,\ldots,x_n) \propto \exp\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n (x_i-\mu)^2\right\}\exp\left\{-\frac{1}{2\sigma_0^2}(\mu-\mu_0)^2\right\}.\] El doble del exponente de la expresión anterior cambiado de signo es igual a \[\frac{1}{\sigma^2}(\sum_{i=1}^n x_i^2 + n\mu^2 -2\mu\sum_{i=1}^n x_i) + \frac{1}{\sigma_0^2}(\mu^2+\mu_0^2-2\mu\mu_0) = A\mu^2-2B\mu + C,\] donde \(A:=n/\sigma^2 + 1/\sigma_0^2\), \(B:=n\bar{x}/\sigma^2 +\mu_0/\sigma_0^2\) y \(C\) no depende de \(\mu\). Como consecuencia, si completamos el cuadrado del exponente, \[\pi(\mu|x_1,\ldots,x_n) \propto \exp\left\{-\frac{1}{2}(A\mu^2-2B\mu)\right\}\propto \exp\left\{-\frac{A}{2}(\mu^2-2\frac{B}{A}\mu+\frac{B^2}{A^2})\right\}.\] Finalmente, \[\pi(\mu|x_1,\ldots,x_n) \propto \exp\left\{-\frac{1}{2/A}(\mu - B/A)^2\right\}.\]

De la última expresión se deduce que \(\mu|x_1,\ldots,x_n\) tiene distribución normal de media \(B/A\) y varianza \(1/A\). Sustituyendo los valores de \(A\) y \(B\) tenemos el resultado.

3.4.3 El uso de métodos de simulación para aproximar los estimadores bayesianos

Desde el punto de vista teórico, el enfoque bayesiano es muy simple. Sin embargo, en función de cuál sea la distribución a priori especificada y el modelo estadístico utilizado, los problemas computacionales que presenta el cálculo de la distribución a posteriori y su esperanza pueden ser muy difíciles, especialmente si \(\theta\) es un vector de alta dimensión.

Tradicionalmente, con el fin de simplificar los cálculos, se elegía una distribución a priori de tal forma que la distribución a posteriori se pudiera identificar fácilmente. Las familias conjugadas para un modelo dado son familias paramétricas de distribuciones a priori tales que la distribución a posteriori pertenece a la misma familia paramétrica, cuando los datos siguen ese modelo. En los ejemplos anteriores hemos comprobado que la distribución beta es conjugada para la binomial en la estimación de una proporción, y que la distribución normal es conjugada para la normal para estimar la media con varianza conocida.

Cuando solo hay un parámetro que toma valores en un intervalo, aunque la distribución a posteriori no sea conocida siempre podemos considerar una aproximación basada en un conjunto de posibles valores representativos de la totalidad del espacio paramétrico \(\{\theta_1,\ldots,\theta_N\}\subset \Theta\). Una vez nos restringimos a este conjunto finito de valores podemos resolver un problema muy similar al primero de los ejemplos que vimos en la sección anterior.

Otra posibilidad es usar métodos de simulación. Veamos un ejemplo en el que un método de simulación permite aproximar un estimador bayesiano muy sencillamente, mientras que el uso de los metodos numéricos tradicionales sería engorroso:

Estimación bayesiana de una diferencia de proporciones

Supongamos que observamos \(X\equiv \mbox{B}(n,p_1)\) e \(Y\equiv\mbox{B}(m,p_2)\). Estamos interesados en estimar la diferencia de proporciones \(\theta = p_1-p_2\). Por ejemplo, \(X\) e \(Y\) podrían representar el número de individuos que se recuperan de cierta enfermedad tras recibir un tratamiento y un placebo, respectivamente. En este caso, \(\theta\) es la diferencia entre las probabilidades de recuperación en ambos grupos. Suponemos a priori que \(p_1\) y \(p_2\) son independientes y uniformes en \([0,1]\), es decir, \(\pi(p_1,p_2)=1\). Por lo tanto, la distribución a posteriori verifica \[\pi(p_1,p_2|X,Y)\propto p_1^X (1-p_1)^{n-X}p_2^Y(1-p_2)^{m-Y},\] lo que implica que \(p_1\) y \(p_2\) son independientes a posteriori, con \(p_1|X,Y\equiv \mbox{Beta}(X+1,n-X+1)\) y \(p_2|X,Y\equiv \mbox{Beta}(Y+1,m-Y+1)\).

Para aproximar la distribución a posteriori de \(\theta\) podemos generar un número grande de pares \((p_1^{*(j)},p_2^{*(j)})\), \(j=1,\ldots,R\), con las distribuciones beta marginales que hemos calculado. Para cada par, se obtiene la diferencia \(\theta^{*(j)} = p_1^{*(j)} - p_2^{*(j)}\). Por la ley de los grandes números, una aproximación de la media a posteriori es el promedio: \[\hat{\theta} = \frac{1}{R}\sum_{j=1}^R \theta^{*(j)}.\] En el código siguiente se implementa este procedimiento si \(X=8\), \(Y=6\), \(n=m=10\).

set.seed(100)

# Datos

X <- 8
Y <- 6
n <- 10
m <- 10

# Generación de números aleatorios

R <- 1000 # número de valores que se generan
p1 <- rbeta(R, X+1, n-X+1)
p2 <- rbeta(R, Y+1, m-Y+1)
theta <- p1 - p2

# Distribución a posteriori aproximada

ggplot(data.frame(theta)) +
  geom_histogram(aes(x = theta), fill = 'olivedrab4', col = 'black', bins = 15)


# Media a posteriori

mean(theta)
#> [1] 0.1595184

El estimador bayesiano es aproximadamente 0.16, algo menor que la diferencia de las proporciones muestrales 0.2. Resulta un valor razonable si tenemos en cuenta que el valor esperado a priori de \(\theta\) es 0. Como ejercicio se puede pensar en cómo se podría resolver el problema analíticamente, sin usar simulación.

Si la dimensión de \(\Theta\) es grande los métodos de simulación tradicionales se encuentran con diversas dificultades. Para estos casos se han desarrollado métodos numéricos basados en simulación de cadenas de Markov (Gibbs sampling y, más en general, métodos MCMC (Markov chain Monte Carlo)) que permiten extender la aplicación de los métodos bayesianos a modelos muy complejos con un número grande de parámetros.

3.5 Ejercicios

Ejercicio 3.1 El siguiente código genera una muestra de 100 datos de una distribución de Cauchy con parámetro de posición \(\theta\):

set.seed(123)
theta <- 10
n <- 100
muestra <- rt(n, 1) + theta

Esta distribución tiene la siguiente función de densidad \[ f(x;\theta) = \frac{1}{\pi}\frac{1}{1+(x-\theta)^2}, \]

que coincide de hecho con la de una distribución \(t\) de Student con 1 grado de libertad trasladada una magnitud \(\theta\).

Calcula el estimador de máxima verosimilitud de \(\theta\). ¿Se parece al valor verdadero?
Lleva a cabo algún experimento de simulación para aproximar la varianza del estimador de máxima verosimilitud.
Compara el resultado de la simulación con la aproximación de la varianza que se obtendría a partir de la teoría asintótica para los estimadores de máxima verosimilitud.

(Indicación: Busca información de funciones de R que sirvan para optimizar funciones o resolver ecuaciones no lineales: optimize, nlm, uniroot, etc.)

Ejercicio 3.2 Se observa una muestra de una distribución de Poisson de parámetro \(\theta\) en la que el parámetro \(\theta\) solo puede tomar alguno de los valores del espacio paramétrico \(\Theta = \{1,\, 2,\, 3\}\). Antes de observar la muestra, la opinión de un experto sobre los distintos valores que puede tomar \(\theta\) queda reflejada en la siguiente distribución de probabilidad a priori sobre \(\Theta\):

\(\theta\)	1	2	3
\(\pi(\theta)\)	1/2	1/4	1/4

Si definimos el estimador bayesiano de \(\theta\) como el valor del parámetro más probable a posteriori tras observar la muestra, y ésta consiste en dos observaciones independientes \(x_1,x_2\) tales que \(x_1+x_2=2\), ¿cuánto vale el estimador bayesiano de \(\theta\)?
Para una muestra que cumple la condición del apartado anterior, ¿cuál es el estimador de máxima verosimilitud de \(\theta\)?

Ejercicio 3.3 Para estimar el número desconocido \(N\) de ejemplares de cierta especie animal en un área determinada se suelen emplear los procedimientos llamados de captura y recaptura. Lo que sigue es un ejemplo del esquema más sencillo posible. Se capturan \(n\) ejemplares de la población, se identifican con una etiqueta, y se vuelven a soltar. Posteriormente, se seleccionan aleatoriamente con reemplazamiento \(m\) ejemplares de la población, y se observa que hay \(r\) de ellos que están etiquetados. Con esta información, ¿cuál es el estimador de máxima verosimilitud de \(N\)?

Ejercicio 3.4 Sean \(X_1,\ldots,X_n\) v.a.i.i.d. de una distribución \(\mbox{Beta}(\theta,1)\), donde \(\theta>0\).

Demuestra que \(-\log X_i\) tiene distribución exponencial de parámetro \(\theta\).
Calcula el estimador de máxima verosimilitud de \(\theta\). ¿Es insesgado?
Calcula la varianza del estimador de máxima verosimilitud. ¿Alcanza la cota de Fréchet-Cramer-Rao?
Considera el siguiente estimador de \(\theta\): \[ \tilde{\theta} = -\frac{n-1}{\sum_{i=1}^n\log X_i}. \] ¿Es \(\tilde{\theta}\) insesgado? ¿Alcanza la cota de Fréchet-Cramer-Rao?

3.6 Referencias

Knight (1999) trata la teoría sobre estimación de máxima verosimilitud de manera parecida a como se ha hecho aquí, y en general es un libro recomendable de nivel intermedio sobre inferencia estadística. Para ampliar información sobre el enfoque bayesiano se puede consultar Clyde et al (2021). El ejemplo acerca del uso de simulación en la inferencia bayesiana para la diferencia de proporciones lo he tomado de Wasserman (2013). Como sugiere su título, en el libro se traza una visión global pero concisa de muchos métodos de inferencia, paramétricos y no paramétricos.

Clyde, M., Çetinkaya-Rundel, M., Rundel, C., Bankd, D., Chai, C., Huang, L. (2021). An Introduction to Bayesian Thinking.
Knight, K. (1999). Mathematical statistics. CRC Press.
Wasserman, L. (2013). All of statistics: a concise course in statistical inference. Springer.