En este capítulo se mostrará como usar R para obtener la función de log-verosimilitud y estimadores por el método de máxima verosimilitud.

El concepto de verosimilitud fue propuesto por @Fisher22 en el contexto de estimación de parámetros. En la Figura \@ref(fig:fisher) se muestra una fotografía de Ronald Aylmer Fisher.

```{r fisher, echo=FALSE, fig.cap='Fotografía de Ronald Aylmer Fisher (1890-1962).', dpi=400, fig.align='center', out.width = '40%'}

La función de verosimilitud para un vector de parámetros $\boldsymbol{\theta}$ dada una muestra aleatoria $\boldsymbol{x}=(x_1, \ldots,x_n)^\top$ con una distribución asumida se define usualmente como:

\begin{equation}

L(\boldsymbol{\theta} | \boldsymbol{x}) = \prod_{i=1}^{n} f(x_i | \boldsymbol{\theta}),

(\#eq:lik)

\end{equation}

donde $x_i$ representa uno de los elementos de la muestra aleatoria y $f(\cdot)$ es la función de masa/densidad de la distribución de la cual se obtuvo $\boldsymbol{x}$.

La función de log-verosimilitud $l(\boldsymbol{\theta} | \boldsymbol{x})$ se define como el logaritmo de la función de verosimilitud $L(\boldsymbol{\theta} | \boldsymbol{x})$, es decir

\begin{equation}

l(\boldsymbol{\theta} | \boldsymbol{x}) = \log L(\boldsymbol{\theta} | \boldsymbol{x}) = \sum_{i=1}^{n} \log f(x_i | \boldsymbol{\theta})

(\#eq:loglik)

\end{equation}

El método de máxima verosimilitud se usa para estimar los parámetros de una distribución. El objetivo de este método es encontrar los valores de $\boldsymbol{\theta}$ que maximizan a $L(\boldsymbol{\theta} | \boldsymbol{x})$ o a $l(\boldsymbol{\theta} | \boldsymbol{x})$, los valores encontrados se representan usualmente por $\hat{\boldsymbol{\theta}}$.

```{block2, type='rmdnote'}

En este ejemplo vamos a considerar la distribución binomial cuya función de masa de probabilidad está dada por:

La distribución binomial anterior tiene sólo un parámetro $p$, por lo tanto en este caso se $\theta=p$.

Suponga que se tiene el vector `rta` que corresponde a una muestra aleatoria de una distribución binomial con parámetro $n=5$ conocido.

1) Calcular el valor de log-verosimilitud $l(\theta)$ si asumiendo que $p=0.30$ en la distribución binomial.

Para obtener el valor de $l(\theta)$ en el punto $p=0.30$ se aplica la definición dada en la expresión \@ref(eq:loglik). Como el problema trata de una binomial se usa entonces la función de masa `dbinom` evaluada en la muestra `rta`, el parámetro `size` como es conocido se reemplaza por el valor de cinco y en el parámetro `prob` se cambia por 0.3. Como interesa la función de log-verosimilitud se debe incluir `log=TRUE`. A continuación el código necesario.

Por lo tanto $l(\theta)= -24.55$

2) Construir una función llamada `ll` a la cual le ingrese valores del parámetro $p$ de la binomial y que la función entregue el valor de log-verosimilitud.

La función solicitada tiene un cuerpo igual al usado en el numeral anterior, a continuación el código necesario para crearla.

Vamos a probar la función en dos valores arbitrarios $p=0.15$ y $p=0.80$ que pertenezcan al dominio del parámetro $p$ de la distribución binomial.

El valor de log-verosimilitud para $p=0.15$ fue de -25.54 mientras que para $p=0.80$ fue de -98.46.

Vamos ahora a chequear si la función `ll` está vectorizada y para esto usamos el código mostrado a continuación y deberíamos obtener un vector con los valores `c(-25.54, -98.56)`.

No obtuvimos el resultado esperado, eso significa que nuestra función no está vectorizada. Ese problema lo podemos solucionar así:

Vemos que ahora que cuando se ingresa un vector a la función `ll` se obtiene un vector.

```{block2, type='rmdnote'}

3) Dibujar la curva log-verosimilitud $l(\theta)$, en el eje X debe estar el parámetro $p$ del cual depende la función de log-verosimilitud.

En la Figura \@ref(fig:loglik1) se presenta la curva solicitada.

```{r loglik1, fig.align='center', fig.cap='Función de log-verosimilitud para el ejemplo sobre binomial.', fig.asp=0.9, fig.width=5, echo=TRUE}

4) Observando la Figura \@ref(fig:loglik1), ¿cuál esl el valor de $p$ que maximiza la función de log-verosimilitud?

Al observar la Figura \@ref(fig:loglik1) se nota que el valor de $p$ que maximiza la función log-verosimilitud está muy cerca de 0.2.

5) ¿Cuál es el valor exacto de $p$ que maximiza la función log-verosimilitud?

En R existe la función `optimize` que sirve para encontrar el valor que __minimiza__ una función uniparamétrica en un intervalo dado, sin embargo, aquí interesa es maximimizar la función de log-verosimilitud, por esa razón se construye la función `minusll` que es el negativo de la función `ll` para así poder usar `optimize`. A continuación el código usado. \index{optimize}

Del resultado anterior se observa que cuando $p=0.23$ el valor máximo de log-verosimilitud es -23.32 (negativo de `minusll`).

Suponga que la estatura de una población se puede asumir como una normal $N(170, 25)$. Suponga también que se genera una muestra aleatoria de 50 observaciones de la población con el objetivo de recuperar los valores de la media y varianza poblacionales a partir de la muestra aleatoria.

La muestra se va a generar con la función `rnorm` pero antes se fijará una semilla con la intención de que el lector pueda replicar el ejemplo y obtener la misma muestra aleatoria aquí generada, el código para hacerlo es el siguiente.

1) Construya la función de log-verosimilitud para los parámetros de la normal dada la muestra aleatoria `y`.

Abajo se muestra la forma de construir la función de log-verosimilitud.

```{block2, type='rmdwarning'}

2) Dibujar la función de log-verosimilitud.

En la Figura \@ref(fig:loglik2) se muestra el gráfico de niveles para la superficie de log-verosimilitud. De esta figura se nota claramente que los valores que maximizan la superficie están alrededor de $\mu=170$ y $\sigma=5$.

```{r loglik2, fig.cap='Gráfico de niveles para la función de log-verosimilitud para el ejemplo sobre normal.', fig.asp=0.6, fig.width=9}

3) Obtenga los valores de $\mu$ y $\sigma$ que maximizan la función de log-verosimilitud.

Para obtener los valores solicitados vamos a usar la función `nlminb` que es un optimizador. A la función `nlminb` se le debe indicar por medio del parámetro `objective` la función que queremos optimizar (minimizar); el parámetro `start` es un vector con los valores iniciales para comenzar la búsqueda de $\mu$ y $\sigma$; los parámetros `lower` y `upper` sirven para delimitar el espacio de búsqueda. A continuación se muestra el código usado para obtener los valores que minimizan a `minusll`, es decir, los valores que maximizan la función de log-verosimilitud. \index{nlminb}

De la salida anterior podemos observar que los valores óptimos de $\mu$ y $\sigma$ son 170.338 y 5.424 respectivamente, resultado que coincide con lo observado en la Figura \@ref(fig:loglik2) y con los valores reales de simulación de la muestra. Esto indica que el procedimiento de estimación de parámetros por máxima verosimilitud entrega valores insesgados de los parámetros a estimar.

Un resultado interesante de la salida anterior es que se reporta el valor mínimo que alcanza la función `minusll`, este valor fue de 155.5, por lo tanto, se puede afirmar que el valor máximo de log-verosimilitud es -155.5.

Otros resultados importantes de la salida anterior son el valor de `convergence=0` que indica que la búsqueda fue exitosa; `iterations=13` indica que se realizaron 13 pasos desde el punto inicial `start` hasta las coordenadas de optimización.

```{block2, type='rmdnote'}

4) ¿Hay alguna función para obtener directamente el valor que maximiza la función log-verosimilitud?

La respuesta es si. Si la distribución estudiada es una de las distribuciones básicas se puede usar la función `fitdistr` del paquete básico **MASS**. Esta función requiere de los datos que se ingresan por medio del parámetro `x`, y de la distribución de los datos que se ingresa por medio del parámetro `densfun`. La función `fitdistr` admite 15 distribuciones diferentes para hacer la búsqueda de los parámetros que caracterizan una distribución, se sugiere consultar la ayuda de la función `fitdistr` escribiendo en la consola `help(fitdistr)`. A continuación el código usado. \index{fitdistr}

```{r, message=FALSE}

El objeto `res` contiene los resultados de usar `fitdistr`. En la primer línea están los valores de los parámetros que maximizan la función de log-verosimilitud, en la parte de abajo, dentro de paréntesis, están los errores estándar o desviaciones de éstos estimadores.

Al objeto `res` es de la clase _fitdistr_ y por lo tanto se le puede aplicar la función genérica `logLik` para obtener el valor de la log-verosimilitud. Se sugiere consultar la ayuda de la función `logLik` escribiendo en la consola `help(logLik)`. A continuación el código para usar `logLik` sobre el objeto `res`.

De esta última salida se observa que el valor coincide con el obtenido cuando se usó `nlminb`.

Generar $n=100$ observaciones de una gamma con parámetro de forma igual a 2, parámetro de tasa igual a 0.5 y luego responder las preguntas.

1) ¿Cómo se puede generar la muestra aleatoria solicitada?

Para generar la muestra aleatoria (`ma`) solicitada se fijó la semilla con el objetivo de que el lector pueda obtener los mismos resultados de este ejemplo.

2) Asumiendo que la muestra aleatoria proviene de una normal (lo cual es incorrecto), estime los parámetros de la distribución normal.

3) Asumiendo que la muestra aleatoria proviene de una gamma estime los parámetros de la distribución gamma.

```{r, warning=FALSE}

En la salida anterior están los valores estimados de los parámetros de la distribución por el método de máxima verosimilitud, observe la cercanía de éstos con los verdaderos valores de 2 y 0.5 para forma y tasa respectivamente.

4) Dibuje dos qqplot, uno asumiendo distribución normal y el otro distribución gamma. ¿Cuál distribución se ajusta mejor a los datos simulados?

Para dibujar el qqplot se usa la función genérica `qqplot`, recomendamos consultar @hernandez_correa para los detalles de cómo usar esta función. Al usar `qqplot` para obtener el qqplot normal y gamma es necesario indicar los valores $\hat{\boldsymbol{\theta}}$ obtenidos en el numeral anterior, por eso es que en el código mostrado a continuación aparece `mean=4.3083, sd=2.8085` en el qqplot normal y `shape=2.23978, rate=0.51988` en el qqplot gamma.

```{r normgamma, fig.cap='Gráfico cuantil cuantil normal y gamma para la muestra simulada.', fig.asp=0.6, fig.width=9}

En la Figura \@ref(fig:normgamma) se muestran los qqplot solicitados. Se observa claramente que al asumir normalidad (lo cual es incorrecto), los puntos del qqplot no están alineados, mientras que al asumir distribución gamma (lo cual es correcto), los puntos si están alineados. De esta figura se concluye que la muestra `ma` puede provenir de una $Gamma(2.23978, 0.51988)$.

```{block2, type='rmdtip'}

```{block2, type='rmdwarning'}

5) Para comparar modelos se puede utilizar el _Akaike information criterion_ ($AIC$) propuesto por @Akaike74 que sirve para medir la calidad relativa de los modelos estadísticos, la expresión para calcular el indicador es $AIC=-2 \, \hat{l}+2 \, df$, donde $\hat{l}$ corresponde al valor de $\log$-verosimilitud y $df$ corresponde al número de parámetros estimados del modelo. Siempre el modelo elegido es aquel modelo con el menor valor de $AIC$. Calcular el $AIC$ para los modelos asumidos normal y gamma.

De los resultados anteriores se concluye que entre los dos modelos, el mejor es el gamma porque su $AIC=466$ es el menor de toos los $AIC$.

```{block2, type='rmdnote'}

En esta sección se explican los conceptos y utilidad de la función Score y la Información de Fisher.

La función Score denotada por $S(\theta)$ se define como la primera derivada de la función de log-verosimilitud así:

y el estimador de máxima verosimilitud $\hat{\theta}$ se encuentra solucionando la igualdad

En el valor máximo $\hat{\theta}$ la curva $l(\theta)$ es cóncava hacia abajo y por lo tanto la segunda derivada es negativa, así la curvatura $I(\theta)$ se define como

Una curvatura grande $I(\hat{\theta})$ está asociada con un gran pico en la función de log-verosimilitud y eso significa una menor incertidumbre sobre el parámetro $\theta$ [@pawitan_2013]. En particular la varianza del estimador de máxima verosimilitud está dada por

Suponga que se desea estudiar una variable que tiene distribución Poisson con parámetro $\lambda$ desconocido. Suponga además que se tienen dos situciones:

- Un solo valor 5 para estimar $\lambda$.

- Cuatro valores 5, 10, 6 y 15 para estimar $\lambda$.

Dibujar la función $l(\lambda)$ para ambos casos e identificar la curvatura.

A continuación el código para evaluar la función $l(\lambda)$ para cada caso.

En la Figura \@ref(fig:loglikpoisson) se muestran las dos curvas $l(\lambda)$ para cada uno de los casos. De la figura se observa claramente que cuando se tienen 4 observaciones la curva es más puntiaguda y por lo tanto menor incertibumbre sobre el parámetro $\lambda$ a estimar.

```{r loglikpoisson, fig.cap='Curvas de log-verosimilitud para los dos casos.', fig.height=4, fig.width=5, echo=FALSE}

En esta sección se mostrará como estimar los parámetros de un modelo de regresión general.

Considere el modelo de regresión mostrado abajo. Simule 1000 observaciones del modelo y use la función `optim` para estimar los parámetros del modelo.

\begin{align*}

y_i &\sim N(\mu_i, \sigma^2), \\

\mu_i &= -2 + 3 x_1, \\

\sigma &= 5, \\

x_1 &\sim U(-5, 6).

\end{align*}

El código mostrado a continuación permite simular un conjunto de valores con la estructura anteior.

El vector de parámetros del modelo anterior es $\boldsymbol{\theta}=(\beta_0, \beta_1, \sigma)^\top=(-2, 3, 5)^\top$, el primer elemento corresponde al intercepto, el segundo a la pendiente y el último a la desviación.

Ahora vamos a usar la función `optim` para encontrar los valores que maximizan la función de log-verosimilitud, el código para hacer eso se muestra a continuación. En el parámetro `par` se coloca un vector de posibles valores de $\boldsymbol{\theta}$ para iniciar la búsqueda, en `fn` se coloca la función de interés, en `lower` y `upper` se colocan vectores que indican los límites de búsqueda de cada parámetro, los $\beta_k$ pueden variar entre $-\infty$ y $\infty$ mientras que el parámetro $\sigma$ toma valores en el intervalo $(0, \infty)$. Como la función `minusll` tiene argumentos adicionales `y` e `x1`, estos pasan a la función `optim` al final como se muestra en el código.

En el objeto `res1` está el resultado de la optimización, para explorar los resultados usamos

De la salida anterior se observa que el vector de parámetros estimado es $\hat{\beta}_0 = `r res1$par[1]`$, $\hat{\beta}_1 = `r res1$par[2]`$ y $\hat{\sigma} = `r res1$par[3]`$, se observa también que el valor de la máxima log-verosimilitud fue de `r -res1$value`. Vemos entonces que el vector estimado está muy cerca del verdadero $\boldsymbol{\theta}=(\beta_0=-2, \beta_1=3, \sigma=5)^\top$.

```{block2, type='rmdnote'}

En algunas ocasiones es mejor hacer la búsqueda de los parámetros en el intervalo $(-\infty, \infty)$ que en una región limitada como por ejemplo $(0, \infty)$ o $(-1, 1)$, ya que las funciones de búsqueda podrían tener problemas en los bordes de esos intervalos. Una estrategia usual en este tipo de casos es aplicar una transformación apropiada al parámetro que tiene el dominio limitado. En el presente ejemplo $\sigma$ sólo puede tomar valores mayores que cero y una transformación de tipo $\log$ podría ser muy útil ya que $\log$ relaciona los reales positivos con todos los reales. La transformación para este problema sería $\log(\sigma)=\beta_3$ o escrita de forma inversa $\sigma=\exp(\beta_3)$. El nuevo parámetro $\beta_3$ puede variar en $(-\infty, \infty)$ pero al ser transformado por la función exponencial este se volvería un valor apropiado para $\sigma$. Para implementar esta variación lo único que se debe hacer es modificar la línea 3 de la función `minusll` como se muestra a continuación:

Para hacer la búsqueda se procede de forma similar, abajo el código necesario.

De la salida anterior se observa que el vector de parámetros estimado es $\hat{\beta}_0 = `r res2$par[1]`$, $\hat{\beta}_1 = `r res2$par[2]`$ y $\hat{\sigma} = \exp(`r res2$par[3]`)=`r exp(res2$par[3])`$, se observa también que el valor de la máxima log-verosimilitud fue de `r -res2$value`. Vemos entonces que el vector estimado está muy cerca del verdadero $\boldsymbol{\theta}=(\beta_0=-2, \beta_1=3, \sigma=5)^\top$.

El paquete **maxLik** tiene un conjunto de funciones útiles para estación por máxima verosimilitud y varias viñetas que pueden ser consultadas [en este enlace](https://cran.r-project.org/web/packages/maxLik/index.html).

En este ejemplo vamos a repetir un ejemplo anterior con la distribución normal.

Ahora vamos a crear la función de log-verosimilitud de la siguiente manera:

Ahora vamos a estimar los parametros con el paquete **maxLik**.

```{r message=FALSE}

1) Al inicio del Capítulo \@ref(central) se presentó la base de datos sobre medidas del cuerpo, consulte la explicación sobre la base de datos y responda lo siguiente.

+ Si se asume que la `edad` tiene distribución normal, ¿cuáles son los estimadores de máxima verosimilitud para $\mu$ y $\sigma$?

+ Como el histograma para la edad muestra un sesgo a la derecha se podría pensar que la distribución gamma sería una buena candidata para explicar las edades observadas. Asumiendo una distribución gamma, ¿cuáles son los estimadores de máxima verosimilitud para los parámetros?

+ ¿Cuál de los dos modelos es más apropiado para explicar la variable de interés? Calcule el $AIC$ para decidir.

2) En el capítulo \@ref(discretas) se presentó un ejemplo donde se usó la base de datos sobre cangrejos hembra. Consulte la explicación sobre la base de datos y responda lo siguiente.

+ Suponga que el número de satélites sobre cada hembra es una variable que se distribuye Poisson. Construya en R la función de log-verosimilitud $l$, dibuje la función $l$ y encuentre el estimador de máxima verosimilitud de $\lambda$.

+ Repita el ejercicio anterior asumiendo que el número de satélites se distribuye binomial negativo.

+ ¿Cuál de los dos modelos es más apropiado para explicar la variable de interés? Calcule el $AIC$ para decidir.

3) Al inicio del Capítulo \@ref(varia) se presentó la base de datos sobre apartamentos usados en Medellín, consulte la explicación sobre la base de datos y responda lo siguiente.

+ Dibuje una densidad para la variable área del apartamento.

+ Describa lo encontrado en esa densidad.

+ ¿Qué distribuciones de 2 parámetros podrían explicar el comportamiento del área de los apartamentos? Mencione al menos 3.

+ Para cada una de las distribuciones anteriores dibuje un gráfico de contornos o calor para la función de log-verosimilitud y estime los parámetros de la distribución elegida.

+ ¿Cuál de los dos modelos es más apropiado para explicar la variable de interés? Calcule el $AIC$ para decidir.

4) Considere el siguiente modelo de regresión.

\begin{align*}

y_i &\sim Gamma(shape_i, scale_i), \\

\log(shape_i) &= 3 - 7 x_1, \\

\log(scale_i) &= 3 - 1 x_2, \\

x_1 &\sim U(0, 1), \\

x_2 &\sim Poisson(\lambda=3)

\end{align*}

+ Simule 100 observaciones del modelo anterior.

+ Escriba el vector de parámetros del problema.

+ Construya la función `minusll` para el problema.

+ Use la función `optim` para estimar los parámetros del problema.