Función de Verosimilitud

Recordemos que si \(x_1,x_2,\ldots,x_n\) es una m.a.s. extraída de una población \(X\) con distribución de probabilidad dada por \(f(x|\theta)\), siendo \(\theta\) un (o más) parárametro poblacional deconocido, entonces la función de verosimilitud de la muestra se define mediante la expresión:

\[\begin{equation} L(\theta)=f(x_1,x_2,\ldots,x_n|\theta) \end{equation}\]

Ahora bien, como cada \(x_i\) es una realización de la v.a. \(X_i\) y estas son independientes e idénticamente distribuidas (iid), la función de verosimilitud puede escribirse como sigue:

\[\begin{equation} L(\theta)=\prod\limits_{i=1}^n f(x_i | \theta) \end{equation}\]

Estimación por el Método de Máxima Verosimilitud

La teoría nos indica que para encontrar el estimador máximo verosímil (EMV) para el parámetro poblacional desconocido \(\theta\), debemos resolver el siguiente problema de maximización:

\[\begin{equation} \max_{\widehat{\theta}}\; L(\theta)\quad\text{o bien} \quad \max_{\widehat{\theta}}\; \ln L(\theta) \end{equation}\]

Sin embargo, computacionalmente hablando, el problema que realmente se resuelve es minimizar el negativo de la función de log-verosimilitud: \[\begin{align} -\ln L(\theta)&= -\ln \prod\limits_{i=1}^n f(x_i | \theta)\\ &=-\sum\limits_{i=1}^n \ln f(x_i|\theta) \end{align}\]

Lo anterior es posible gracias a que \(\max \,Z = \min \,(-Z)\), siendo \(Z\) una función objetivo cualquiera.

Implementación en R

Antes de comenzar, es necesario instalar y cargar algunos paquetes relevantes.

#install.packages("bbmle")
#install.packages("stats4")
library(stats4) # para la función mle
library(bbmle) # para la función mle2

Función mle del paquete stats4

Comenzaremos generando 50 datos aleatorios según una distribución normal con media 10 y desviación estándar 2. Guardamos estos datos en un vector llamado x. Nótese que previamente usamos la opción set.seed(123) para fijar los valores generados.

# removemos los warnings para no visualizarlos

set.seed(123)
x = rnorm(n = 50,mean = 10,sd = 2)

Ahora construimos el negativo de la función de log-verosimilitud. Esta función, llamada NegLogLik en este caso, depende de los parámetros poblacionales desconocidos media y desviación estándar, mu y sigma respectivamente:

NegLogLik = function(mu,sigma){-sum(dnorm(x,mu,sigma,log = TRUE))}

Luego usamos la función mle del paquete stats4 para realizar la estimación por máxima verosimilitud:

EMV1 = mle(NegLogLik, start = list(mu=10, sigma=5))

El argumento star=list() es para ingresar los valores iniciales con los cuales R construirá numéricamente una aproximación de la Matriz Hessiana asociada al problema de optimización (recuerde el curso de Cálculo III).

Finalmente, visualizamos los resultados obtenidos en la estimación con la opción summary():

summary(EMV1)
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(mu = 10, sigma = 5))
## 
## Coefficients:
##       Estimate Std. Error
## mu    10.06881  0.2592437
## sigma  1.83313  0.1833128
## 
## -2 log L: 202.4963

Los resultados obtenidos nos indican que la estimación puntual obtenida por el método de máxima verosimilitud para \(\mu\) y \(\sigma\) es

\[ \begin{array}{lcr} \widehat{\mu}&=&10.07\\ \widehat{\sigma}&=&1.83 \end{array} \]

También se puede observar que los errores estándares computados son bastante pequeños, lo cual siempre es un buen síntoma. Nótese además que los valores obtenidos son muy cercamos a los declarados al momento de generar aleatoriamente los datos utilizados para construir la estimación puntual.

Podemos además generar un intervalo de confianza del 95% para los verdaderos valores de los parámetros poblacionales desconocidos:

confint(EMV1)
## Profiling...
##          2.5 %    97.5 %
## mu    9.550743 10.586871
## sigma 1.525189  2.260819

Por otro lado, en el ejemplo anterior se utilizaron datos generados aleatoriamente con ayuda de las distribuciones de probabilidad integradas en R. Ahora, si solo contamos con datos sueltos como por ejemplo 9, 10, 5, 3, 6, 8 y 10, los cuales suponemos siguen una distribución normal, el proceso para realizar la estimación por máxima verosimilitud es similar al anterior.

Lo primero que hacemos es construir un vector y en el cual guardamos los datos muestrales:

y = c(9,10,5,3,6,8,10)

Construimos el negativo de la función de log-verosimilitud respectiva:

NegLogLik2 = function(mu,sigma){-sum(dnorm(y,mu,sigma,log = TRUE))}

Realizamos la estimación por máxima verosimilitud:

EMV2 = mle(NegLogLik2, start = list(mu=10, sigma=5))

Y finalmente visualizamos los resultados obtenidos en la estimación:

summary(EMV2)
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(mu = 10, sigma = 5))
## 
## Coefficients:
##       Estimate Std. Error
## mu    7.285717  0.9414340
## sigma 2.490800  0.6656942
## 
## -2 log L: 32.64159

Función mle2 del paquete bbmle

La función mle2 del paquete bbmle en su forma más básica es prácticamente igual a la función mle del paquete stats4. La principal diferencia radica en que debemos indicar en el argumento data de donde se tomarán los datos para realizar los cálculos:

EMV3 = mle2(NegLogLik,start = list(mu=10,sigma=5), data = list(x))

Visualizamos los datos y vemos que los resultados son prácticamente iguales a los obtenidos anteriormente. No obstante, la función mle2 entrega más información que la función mle, particularmente, el término Pr(z), el cual será discutido más adelante.

summary(EMV3)
## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(mu = 10, sigma = 5), 
##     data = list(x))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## mu    10.06881    0.25924  38.839 < 2.2e-16 ***
## sigma  1.83313    0.18331  10.000 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 202.4963

La principal ventaja de la función mle2 es que con ella se puede prescindir de construir la función de verosimilitud, ya que es posible pedirle a la función que la calcule internamente usando el símbolo ~.

A modo de ejemplo, generemos una secuencia de 30 números según una distribución normal con media igual a 20 y desviación típica igual a 5:

set.seed(12345)
x1=rnorm(30, mean=20, sd=5)
EMV4 = mle2(x1~dnorm(mu,sigma),start = list(mu=15,sigma=4), data=data.frame(x1))
summary(EMV4)
## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = x1 ~ dnorm(mu, sigma), start = list(mu = 15, 
##     sigma = 4), data = data.frame(x1))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## mu    20.39395    0.84206 24.2190 < 2.2e-16 ***
## sigma  4.61218    0.59541  7.7462 9.467e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 176.8597

Notemos que la notación x~dnorm(mu,sigma) significa que x es en realidad el primer argumento de la función dnorm(). De igual forma, es importante destacar que al usar esta metodología, los datos que se utilizarán para realizar los cómputos se deben ingresar como un data.frame().

Ejercicios

Use R para implementar y resolver los siguientes problemas: