15 Verosimilitud

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.

15.1 Función de verosimilitud

El concepto de verosimilitud fue propuesto por Fisher (1922) en el contexto de estimación de parámetros. En la Figura 15.1 se muestra una fotografía de Ronald Aylmer Fisher.

Fotografía de Ronald Aylmer Fisher (1890-1962).

Figure 15.1: Fotografía de Ronald Aylmer Fisher (1890-1962).

Asumiendo un modelo estadístico parametrizado por una cantidad fija y desconocida \(\theta\), la verosimilitud \(L(\theta)\) es la probabilidad de los datos observados \(x\) como una función de \(\theta\) (Pawitan 2013). Si la variable de interés es discreta se usa la probabilidad y si es continua se usa la densidad para obtener la verosimilitud.

La función de verosimilitud para un vector de parámetros \(\boldsymbol{\Theta}\) dada una muestra aleatoria \(\boldsymbol{x}\) 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}), \tag{15.1} \end{equation}\]

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

15.2 Función de log-verosimilitud

La función de log-verosimilitud \(l\) se define como el logaritmo de la función de verosimilitud \(L\), 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}) \tag{15.2} \end{equation}\]

15.3 Método de máxima verosimilitud para estimar parámetros

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 \(L\) o \(l\) y valores encontrados se representan por \(\hat{\boldsymbol{\Theta}}\).

Ejemplo

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

\[f(x)=P(X=x)=\binom{n}{x} p^x (1-p)^{n-x}, \quad 0<p<1, \quad n \leq 1, 2, \ldots, \quad 0 \leq x \leq n\]

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

rta <- c(2, 2, 1, 1, 1, 1, 0, 2, 1, 2,
         1, 0, 1, 2, 1, 0, 0, 2, 2, 1)
  1. Calcular el valor de log-verosimilitud \(l\) si asumiendo que \(p=0.30\) en la distribución binomial.

Para obtener el valor de \(l\) en el punto \(p=0.30\) se aplica la definición dada en la expresión (15.2). 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.

sum(dbinom(x=rta, size=5, prob=0.3, log=TRUE))
## [1] -24.55

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

  1. 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.

ll <- function(prob) sum(dbinom(x=rta, size=5, prob=prob, log=T))

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.

ll(prob=0.15)  # Individual para p=0.15
## [1] -25.54
ll(prob=0.80)  # Individual para p=0.80
## [1] -98.46

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).

ll(prob=c(0.15, 0.80))
## [1] -57.32

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

ll <- Vectorize(ll)
ll(prob=c(0.15, 0.80))
## [1] -25.54 -98.46

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

Necesitamos que la función ll esté vectorizada para poder dibujarla y para poder optimizarla.
  1. Dibujar la curva log-verosimilitud \(l\), en el eje X debe estar el parámetro \(p\) del cual depende la función de log-verosimilitud.

En la Figura 15.2 se presenta la curva solicitada.

Función de log-verosimilitud para el ejemplo sobre binomial.

Figure 15.2: Función de log-verosimilitud para el ejemplo sobre binomial.

  1. Observando la Figura 15.2, ¿cuál esl el valor de \(p\) que maximiza la función de log-verosimilitud?

Al observar la Figura 15.2 se nota que el valor de \(p\) que maximiza la función log-verosimilitud está muy cerca de 0.2.

  1. ¿Cuál es el valor exacto de \(p\) que maximiza la función log-verosimilitud?

En 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.

minusll <- function(x) -ll(x)

optimize(f=minusll, interval=c(0, 1))
## $minimum
## [1] 0.23
## 
## $objective
## [1] 23.32

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

Ejemplo

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.

set.seed(1235)  # La semilla es 1235
y <- rnorm(n=50, mean=170, sd=5)
y[1:7]  # Para ver los primeros siete valores generados
## [1] 166.5 163.6 174.9 170.6 170.6 178.5 170.2
  1. Construya la función de log-verosimilitud para los parámetros de la normal dada la muestra aleatoria y.
ll <- function(param) {
  media <- param[1]  # param es el vector de parámetros
  desvi <- param[2] 
  sum(dnorm(x=y, mean=media, sd=desvi, log=TRUE))
}
Siempre que el interés sea encontrar los valores que maximizan una función de log-verosimilitud, los parámetros de la distribución deben ingresar a la función ll como un vector. Esto se debe hacer para poder usar las funciones de búsqueda optim y nlminb.
  1. Dibujar la función de log-verosimilitud.

En la Figura 15.3 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\). Para ver el código usado en la creación de la Figura 15.3 se recomienda consultar la sección sobre contour en Hernández (2018).

Gráfico de niveles para la función de log-verosimilitud para el ejemplo sobre normal.

Figure 15.3: Gráfico de niveles para la función de log-verosimilitud para el ejemplo sobre normal.

  1. 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.

minusll <- function(x) -ll(x)
nlminb(objective=minusll, start=c(163, 3.4),
       lower=c(160, 3), upper=c(180, 7))
## $par
## [1] 170.338   5.424
## 
## $objective
## [1] 155.5
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 13
## 
## $evaluations
## function gradient 
##       16       35 
## 
## $message
## [1] "relative convergence (4)"

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 15.3 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.

En se tienen dos funciones básicas para optimizar funciones, es decir, para encontrar los valores que minimizan una función dada. Esas dos funciones son nliminb y optim. Para optimizar en una sola dimensión se usa la función optimize.
  1. ¿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. La función fitdistr 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.

require(MASS) # El paquete ya está instalado, solo se debe cargar
res <- fitdistr(x=y, densfun='normal')
res
##      mean        sd   
##   170.3384     5.4235 
##  (  0.7670) (  0.5424)

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.

logLik(res)
## 'log Lik.' -155.5 (df=2)

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

Ejemplo

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

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.

n <- 100
set.seed(12345)
ma <- rgamma(n=n, shape=2, rate=0.5)
  1. Asumiendo que la muestra aleatoria proviene de una normal (lo cual es incorrecto) estime los parámetros de la distribución normal.
fit1 <- fitdistr(x=ma, densfun='normal')
fit1
##     mean      sd  
##   4.3083   2.8085 
##  (0.2808) (0.1986)
  1. Asumiendo que la muestra aleatoria proviene de una gamma estime los parámetros de la distribución.
fit2 <- fitdistr(x=ma, densfun='gamma')
fit2
##     shape     rate  
##   2.23978   0.51988 
##  (0.29620) (0.07703)

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.

  1. 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 Hernández (2018) 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.

par(mfrow=c(1, 2))

qqplot(y=ma, pch=19,
       x=qnorm(ppoints(n), mean=4.3083, sd=2.8085),
       main='Normal Q-Q Plot',
       xlab='Theoretical Quantiles',
       ylab='Sample Quantiles')

qqplot(y=ma, pch=19,
       x=qgamma(ppoints(n), shape=2.23978, rate=0.51988),
       main='Gamma Q-Q Plot',
       xlab='Theoretical Quantiles',
       ylab='Sample Quantiles')
Gráfico cuantil cuantil normal y gamma para la muestra simulada.

Figure 15.4: Gráfico cuantil cuantil normal y gamma para la muestra simulada.

En la Figura 15.4 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)\).

Para obtener el gráfico cuantil cuantil bajo normalidad se puede usar directamente la función qqnorm, consultar Hernández (2018) para mayores detalles.
En este ejemplo se eligió la mejor distribución entre dos candidatas usando una herramienta gráfica, lo que se recomienda usar algún método menos subjetivo (cuantitativo) para tomar decisiones.
  1. Para comparar modelos se puede utilizar el Akaike information criterion (\(AIC\)) propuesto por Akaike (1974) 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.
-2 * logLik(fit1) + 2 * 2  # AIC para modelo normal
## 'log Lik.' 494.3 (df=2)
-2 * logLik(fit2) + 2 * 2  # AIC para modelo gamma
## 'log Lik.' 466 (df=2)

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\).

Modelos anidados pueden ser comparados por medio del global deviance (\(GD\)) dado por \(GD=-2 \, \hat{l}\) y modelos no anidados por medio del Generalized Akaike information criterion (\(GAIC\)) propuesto por Akaike (1983) y dado por \(GAIC=-2 \, \hat{l} + \sharp \, df\) siendo \(\sharp\) el valor de penalidad por cada parámetro adicional en el modelo; cuando \(\sharp = 2\), el \(GAIC\) coincide con el \(AIC\) y el Schwarz Bayesian criterion (\(SBC\)) propuesto por Schwarz (1978) se dá cuando el valor de penalidad es \(\sharp = \log(n)\) donde \(n\) es el número de observaciones del modelo; siempre el modelo elegido es aquel modelo con el menor valor de cualquiera de los criterios de información anteriores.

Ejemplo

¿Será posible explorar entre muchas distribuciones estadísticas aquella que mejor explica una variable de interés?

La respuesta es si, se puede usar la función fitDist del paquete gamlss . Esta función tiene la siguiente estructura:

fitDist(y, k = 2, 
      type = c("realAll", "realline", "realplus",
               "real0to1", "counts", "binom"))

El parámetro y sirve para ingresar el vector con la información; k=2 es la penalización por cada parámetro estimado para calcular el \(GAIC\), por defecto es 2; y el parámetro type sirve para indicar el tipo de distribución, los posibles valores son:

  • realAll: para hacer la búsqueda en todas las distribuciones disponibles en gamlss.
  • realline: para variables en \(\Re\).
  • realplus: para variables en \(\Re^+\).
  • real0to1: para variables en el intervalo \((0, 1)\).
  • counts: para variables de conteo.
  • binom: para variables de tipo binomial.
  1. Usar la función fitDist para elegir la distribución que mejor se ajusta a los datos simulados en el ejemplo anterior, usar una penalizacion de \(k=2\) para calcular el \(AIC\).

Se usa la función fitDist con type='realplus' porque se sabe que la muestra aleatoria tiene valores sólo en \(\Re^+\), los resultados de almacenan en el objeto modelos. Para obtener la lista de los mejores modelos con su respectivo \(BAIC\) se escribe en la consola modelos$fits. Abajo el código usado.

require(gamlss)
modelos <- fitDist(y=ma, type='realplus', k=2)
## Error in solve.default(oout$hessian) : 
##   Lapack routine dgesv: system is exactly singular: U[4,4] = 0
modelos$fits
##       GA     WEI3     WEI2      WEI       GG      GIG 
##    466.0    466.8    466.8    466.8    468.0    468.0 
##    BCCGo     BCCG     BCPE    BCPEo      GB2     BCTo 
##    468.4    468.4    469.4    469.4    470.0    470.4 
##      BCT   exGAUS   LOGNO2    LOGNO       IG      EXP 
##    470.4    471.2    474.6    474.6    483.9    494.1 
## PARETO2o       GP  PARETO2   IGAMMA 
##    496.1    496.1    496.1    504.2

De la lista anterior se observa que la función gamma está en el primer lugar con un \(BAIC=466\), valor que coincide con el \(AIC\) obtenido en el ejercicio anterior. Esto significa que como la gamma tiene el menor \(BAIC\), es la distribución que mejor se ajusta a los datos de la muestra, esto coincide con la realidad, ya que la muestra fue generada de una distribución gamma.

¿Por queé \(AIC\) y \(BAIC\) coinciden en este ejemplo?

  1. ¿Cuáles son los valores estimados de los parámetros para la distribución gamma identificada en el paso anterior?

En gamlss los parámetros de las distribuciones se nombran de una forma especial, el parámetro de posición se denomina \(\mu\), el de dispersión se denomina \(\sigma\) y los de simetría y curtosis con \(\nu\) y \(\tau\). La distribución gamma sólo tiene 2 parámetros y para obtener sus valores estimados se usa el siguiente código.

modelos$mu
## [1] 4.308
modelos$sigma
## [1] 0.6682

Si comparamos los anteriores resultados, \(\hat{\mu}= 4.308\) y \(\hat{\sigma}=0.6682\), con los encontrados al usar la función fitdist, \(\widehat{shape}=2.23978\) y \(\widehat{rate}=0.51988\), vemos que no coinciden, esto se debe a que la parametrización de la gamma usada por gamlss es diferente a la parametrización usada por la base de . Lo anterior no es motivo de preocupación, lo que se recomienda es que el usuario elija una de las parametrizaciones y trabaje con ella, no debe mezclarlas.

  1. Dibujar el histograma para la muestra aleatoria y agregar la densidad de la distribución gamma identificada como la distribución que mejor explica el comportamiento de la variable.

Para hacer lo solicitado se usa la función histDist del paquete gamlss, sólo es necesario ingresar los datos y la familia o distribución de la densidad requerida. Abajo el código usado.

histDist(y=ma, family=GA, main='', xlab='x', ylab='Densidad',
         line.col='blue', line.wd=4)
Histograma para la muestra simulada con la densidad de una Gamma(mu=4.308, sigma=0.6682).

Figure 15.5: Histograma para la muestra simulada con la densidad de una Gamma(mu=4.308, sigma=0.6682).

## 
## Family:  c("GA", "Gamma") 
## Fitting method: "nlminb" 
## 
## Call:  
## gamlssML(formula = ma, family = "GA", data = sys.parent()) 
## 
## 
## Mu Coefficients:
## [1]  1.46
## Sigma Coefficients:
## [1]  -0.403
## 
##  Degrees of Freedom for the fit: 2 Residual Deg. of Freedom   98 
## Global Deviance:     462 
##             AIC:     466 
##             SBC:     471.3

En la Figura 15.5 se presenta el histograma para muestra aleatoria y la densidad de la gamma que mejor explica estos datos. Se observa claramente que la curva de densidad azul acompaña el patrón de los datos.

15.4 Score e Información de Fisher

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

15.4.1 Score e Información de Fisher en el caso univariado.

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

\[ S(\theta) \equiv \frac{\partial}{\partial \theta} l(\theta) \]

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

\[ S(\theta) = 0 \]

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

\[ I(\theta) \equiv - \frac{\partial^2}{\partial \theta^2} l(\theta) \]

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

\[Var(\hat{\theta})=I^{-1}(\hat{\theta})\]

Ejemplo

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.

# Caso 1
w <- c(5)
ll1 <- function(lambda) sum(dpois(x=w, lambda=lambda, log=T))
ll1 <- Vectorize(ll1)

# Caso 2
y <- c(5, 10, 6, 15)
ll2 <- function(lambda) sum(dpois(x=y, lambda=lambda, log=T))
ll2 <- Vectorize(ll2)

En la Figura 15.6 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.

Curvas de log-verosimilitud para los dos casos.

Figure 15.6: Curvas de log-verosimilitud para los dos casos.

15.5 Método de máxima verosimilitud para estimar parámetros en modelos de regresión

En la sección anterior se estudió cómo estimar los parámetros de una distribución, en este capítulo se aprenderá a estimar los parámetros de un modelo de regresión general.

Ejemplo

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.

n <- 1000
x1 <- runif(n=n, -5, 6)
y <- rnorm(n=n, mean=-2 + 3 * x1, sd=5)

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.

A continuación se define la función de menos log-verosimilitud para el modelo anterior. A pesar de que nos interesa maximizar la función de log-verosimilitud hemos creado su negativo, esto porque la mayoría de las funciones de optimización minimizan y no maximizan; maximizar \(f(x)\) es equivalente a minimizar \(-f(x)\).

minusll <- function(theta, y, x1) {
  media <- theta[1] + theta[2] * x1  # Se define la media
  desvi <- theta[3]                  # Se define la desviación.
  - sum(dnorm(x=y, mean=media, sd=desvi, log=TRUE))
}

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.

res1 <- optim(par=c(0, 0, 1), fn=minusll,
              method='L-BFGS-B',
              lower=c(-Inf, -Inf, 0),
              upper=c(Inf, Inf, Inf), y=y, x1=x1)

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

res1
## $par
## [1] -1.905  3.080  5.014
## 
## $value
## [1] 3031
## 
## $counts
## function gradient 
##       19       19 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

De la salida anterior se observa que el vector de parámetros estimado es -1.9046, 3.0796, 5.0142 y que el valor de la máxima log-verosimilitud fue de 3031.2092. Vemos entonces que el vector estimado está muy cerca del verdadero \(\boldsymbol{\Theta}=(\beta_0, \beta_1, \sigma)^\top=(-2, 3, 5)^\top\).

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 límites numéricos. 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:

minusll <- function(theta, y, x1) {
  media <- theta[1] + theta[2] * x1  
  desvi <- exp(theta[3])  # <<<<<---- El cambio fue aquí
  - sum(dnorm(x=y, mean=media, sd=desvi, log=TRUE))
}

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

res2 <- optim(par=c(0, 0, 0), fn=minusll, y=y, x1=x1)

Para obtener los valores ajustados y almacenados en el objeto res2 usamos el siguiente código.

c(res2$par[1:2], exp(res2$par[3]))
## [1] -1.905  3.079  5.015

De la salida anterior se observa que el vector estimado está muy cerca del verdadero \(\boldsymbol{\Theta}=(\beta_0, \beta_1, \sigma)^\top=(-2, 3, 5)^\top\).

EJERCICIOS

  1. Al inicio del Capítulo 9 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.
  1. En la sección ?? se presentó 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 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.
  1. Al inicio del Capítulo 10 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.
  1. 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.

References

Fisher, R. A. 1922. “On the Mathematical Foundations of Theoretical Statistics.” Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 222. The Royal Society: 309–68. http://www.jstor.org/stable/91208.

Pawitan, Yudi. 2013. In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford University Press.

Hernández, J., F. & Correa. 2018. Gráficos Con R. Primera. Medellín, Colombia: Universidad Nacional de Colombia. http://yihui.name/knitr/.

Akaike, H. 1974. “A new look at the statistical model identification.” IEEE Transactions on Automatic Control 19 (6): 716–23.

Akaike, H. 1983. “Information measures and model selection.” Bulletin of the International Statistical Institute 50: 277–90.

Schwarz, G. 1978. “Estimating the dimension of a model.” Annals of Statistics 6 (2): 461–64.