1. Introducción

2. Función de verosimilitud

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} \mathcal{L}(\boldsymbol{\theta} | \boldsymbol{x}) = \prod_{i=1}^{n} f(x_i| \boldsymbol{\theta}) \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}\).

3. Función de log-verosimilitud

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

\[\begin{equation} \ell(\boldsymbol{\theta} | \boldsymbol{x}) = \log \mathcal{L}(\boldsymbol{\theta} | \boldsymbol{x}) = \sum_{i=1}^{n} \log f(x_i | \boldsymbol{\theta}) \end{equation}\]

4. 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 a \(\mathcal{L}(\boldsymbol{\theta} | \boldsymbol{x})\) o a \(\ell(\boldsymbol{\theta} | \boldsymbol{x})\), los valores encontrados se representan usualmente por \(\hat{\boldsymbol{\theta}}\).

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

A continuación se muestra el procedimiento para obtener el Estimador de Máxima Verosimilitud (EMV) de \(\theta_j\) data una muestra particular

1. Escribir la verosimilitud: \(\mathcal{L}(\theta)=\mathcal{L}\left(x_{1}, \ldots, x_{n} ; \theta\right)\)

2. Escribir el logaritmo de la verosimilitud: \(\ell(\theta)=\ln\mathcal{L}(\theta)\), \(\ell(\theta)\) se llama función soporte o log-verosimilitud.

3. Obtener el \(\theta_{j}\) tal que \[ \frac{\partial}{\partial \theta_{j}} \ell(\theta)=0 \] y denotarlo por \(\hat{\theta}_{j}\).

4. Comprobar que realmente es un máximo, es decir, que \[ \left.\frac{\partial^{2}}{\partial \theta_{j}^{2}} \ell(\theta)\right|_{\theta_{j}=\hat{\theta}_{j}}<0 \]

Ejemplo 1

Determinar la función de verosimilitud y la función log-verosimilitud de una población \(Bernoulli(\theta)\). Representar gráficamente estas funciones cuando \(n=20\) y \(\sum x_i=12\).

Solución: Supongamos \(X_{1}, \ldots X_{n} \sim\) Bernoulli \((\theta)\). La función de densidad correspondiente es \(p(x ; \theta)=\theta^{x}(1-\theta)^{1-x},\) por lo que:

\[ \mathcal{L}(p)=\prod_{i=1}^{n} p\left(x_{i} ; \theta\right)=\prod_{i=1}^{n} \theta^{x_{i}}(1-\theta)^{1-x_{i}}=\theta^{\sum x_{i}}(1-\theta)^{n-\sum x_{i}} \]

denotemos \(S=\sum x_{i},\) entonces

\[ l(\theta)=S \log \theta+(n-S) \log (1-\theta) \]

y obtenemos \(\hat{\theta}=S / n\).

Defiendo la función de verosimilitud

L_bernoulli <- function(n, S){
  function(theta){
    theta ^ S * (1 - theta) ^ (n - S)
  }  
}
#Gráfico de la función de verosimilitud
library(ggplot2)
xy <- data.frame(x = 0:1, y = 0:1)
ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = L_bernoulli(n = 20, S = 12)) +
  xlab(expression(theta)) +
  ylab(expression(L(theta))) +
  ggtitle("Verosimilitud (n=20, S = 12)")

Defiendo la función de log-verosimilitud

# log-verosimilitud
l_bernoulli <- function(n, S){
  function(theta){
    S * log(theta) + (n - S) * log(1 - theta)
  }  
}
#Gráfico de la función de log-verosimilitud
library(ggplot2)
ggplot(xy, aes(x = x, y = y)) +
  stat_function(fun = l_bernoulli(n = 20, S = 12))+
  xlab(expression(theta)) +
  ylab(expression(l(theta))) +
  ggtitle("log-verosimilitud (n=20, S = 12)")

Ejemplo 2

Ejercicio: El tiempo de realización en minutos de una determinada tarea dentro de un proceso industrial es una variable aleatoria con función de densidad \[ f(x)=\frac{x}{\theta^{2}} e^{-x / \theta} \quad \text { si } x>0 \] donde \(\theta>0\).

a) Calcular el estimador máximo-verosímil de \(\theta\) para una muestra aleatoria simple de tamaño \(n .\)

  1. Escribir la verosimilitud: \[ \begin{aligned} \mathcal{L}(\theta) &=\quad \mathcal{L}\left(x_{1}, \ldots, x_{n} ; \theta\right) \stackrel{i n d}{=} \prod_{i=1}^{n} f_{X_{i}}\left(x_{i}\right) \stackrel{i d .}{=} \prod_{i=1}^{n} f_{X}\left(x_{i}\right) \\ &=\prod_{i=1}^{n}\left(\frac{x_{i}}{\theta^{2}} e^{-x_{i} / \theta}\right)=\left(\frac{1}{\theta^{2}}\right)^{n} \prod_{i=1}^{n} e^{-x_{i} / \theta} \prod_{i=1}^{n} x_{i} \\ &=\left(\frac{1}{\theta^{2}}\right)^{n} \exp \left\{-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}\right\} \prod_{i=1}^{n} x_{i} ; \quad x_{1}, \ldots, x_{n}>0 \end{aligned} \]

  2. Escribir el logaritmo de la verosimilitud: \[ \begin{aligned} \ell(\theta) & &=\quad & \ln \mathcal{L}(\theta)=\ln \left(\left(\frac{1}{\theta^{2}}\right)^{n} \exp \left\{-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}\right\} \prod_{i=1}^{n} x_{i}\right) \\ &= & & \ln \left(\left(\frac{1}{\theta^{2}}\right)^{n}\right)+\ln \left(\exp \left\{-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}\right\}\right)+\ln \left(\prod_{i=1}^{n} x_{i}\right) \\ &=& & \ln \left(\theta^{-2 n}\right)-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}+\ln \left(\prod_{i=1}^{n} x_{i}\right) \\ &=& &-2 n \ln \theta-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}+\ln \left(\prod_{i=1}^{n} x_{i}\right) ; \quad x_{1}, \ldots, x_{n}>0 \end{aligned} \]

  3. Obtener el \(\theta\) tal que \(\frac{\partial}{\partial \theta_{j}} \ell(\theta)=0\). \[ \begin{aligned} \frac{\partial}{\partial \theta} \ell(\theta) &=\frac{\partial}{\partial \theta}\left(-2 n \ln \theta-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}+\ln \left(\prod_{i=1}^{n} x_{i}\right)\right) \\ &=\frac{\partial}{\partial \theta}(-2 n \ln \theta)-\frac{\partial}{\partial \theta}\left(\frac{1}{\theta} \sum_{i=1}^{n} x_{i}\right)+\frac{\partial}{\partial \theta}\left(\ln \left(\prod_{i=1}^{n} x_{i}\right)\right) \\ &=\frac{-2 n}{\theta}-\left(-\frac{\sum_{i=1}^{n} x_{i}}{\theta^{2}}\right)+0 \\ &=\frac{-2 n}{\theta}+\frac{\sum_{i=1}^{n} x_{i}}{\theta^{2}} \end{aligned} \]

Por lo tanto: \[ \begin{aligned} \frac{\partial}{\partial \theta} \ell(\theta)=0 & \Leftrightarrow \frac{-2 n}{\theta}+\frac{\sum_{i=1}^{n} x_{i}}{\theta^{2}}=0 \Leftrightarrow \frac{\sum_{i=1}^{n} x_{i}}{\theta^{2}}=\frac{2 n}{\theta} \\ & \Leftrightarrow \frac{\sum_{i=1}^{n} x_{i}}{\theta}=2 n \Leftrightarrow \theta=\frac{\sum_{i=1}^{n} x_{i}}{2 n} \end{aligned} \] El candidato a EMV es \(\hat{\theta}_{M V}=\frac{\sum_{i=1}^{n} x_{i}}{2 n}\)

  1. Comprobar que realmente es un máximo, es decir, que \(\left.\frac{\partial^{2}}{\partial \theta^{2}} \ell(\theta)\right|_{\theta=\hat{\theta}_{M V}}<0\).

\[ \begin{aligned} \frac{\partial^{2}}{\partial \theta^{2}} \ell(\theta) &=\frac{\partial}{\partial \theta}\left(\frac{-2 n}{\theta}+\frac{\sum_{i=1}^{n} x_{i}}{\theta^{2}}\right) \\ &=(-2 n) \frac{-1}{\theta^{2}}+\sum_{i=1}^{n} x_{i} \frac{-2}{\theta^{3}} \\ &=\frac{2 n}{\theta^{2}}-\frac{2}{\theta^{3}} \sum_{i=1}^{n} x_{i}=\frac{2}{\theta^{2}}\left(n-\frac{1}{\theta} \sum_{i=1}^{n} x_{i}\right) \end{aligned} \]

Si ahora lo evaluamos en el candidato \(\hat{\theta}_{M V}=\frac{\sum_{i=1}^{n} x_{i}}{2 n},\) tenemos:

\[ \begin{aligned} \left.\frac{\partial^{2}}{\partial \theta^{2}} \ell(\theta)\right|_{\theta=\bar{\theta}_{M V}} &=\frac{2}{\hat{\theta}_{M V}^{2}}\left(n-\frac{1}{\hat{\theta}_{M V}} \sum_{i=1}^{n} x_{i}\right)=\frac{2}{\hat{\theta}_{M V}^{2}}\left(n-\frac{1}{\frac{\sum_{i=1}^{n} x_{i}}{2 n}} \sum_{i=1}^{n} x_{i}\right) \\ &=\frac{2}{\hat{\theta}_{M V}^{2}}(n-2 n)=\frac{2}{\hat{\theta}_{M V}^{2}}(-n)<0 \end{aligned} \] obtenemos que la segunda derivada de \(\ell(\theta)\) es negativa en \(\hat{\theta}_{M V}\) y por tanto es un máximo.

El estimador máximo verosímil de \(\theta\) es \(\hat{\theta}_{M V}=\frac{\bar{X}}{2}\). Para una muestra particular, la estimación máximo verosímil será \(\hat{\theta}_{M V}=\frac{\bar{x}}{2}\).

Ejemplo 3

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

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(\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 anterior. 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.55231

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.

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)  
## [1] -25.54468
ll(prob=0.80)  
## [1] -98.45598

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

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.54468 -98.45598

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.

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 siguiente figura se presenta la curva solicitada.

curve(ll, lwd=4, col='dodgerblue3',
      xlab='Probabilidad de éxito (p)', las=1,
      ylab=expression(paste("Probabilidad de éxito (p=", theta, ")"))
      )
grid()
Figura1: Función de log-verosimilitud para el ejemplo sobre binomial.

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

4) Observando la Figura anterior, ¿cuál es el valor de \(p\) que maximiza la función de log-verosimilitud?

Al observar la Figura anterior 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 maximizar 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.229993
## 
## $objective
## [1] 23.3246

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

Ejemplo 4

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)  
y <- rnorm(n=50, mean=170, sd=5)
y[1:7]  
## [1] 166.5101 163.5757 174.9498 170.5589 170.5710 178.4910 170.2392

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.

ll <- function(param) {
  media <- param[1]  
  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.

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

En la siguiente figura 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\).

ll1 <- function(a, b) sum(dnorm(x=y, mean=a, sd=b, log=TRUE))
ll1 <- Vectorize(ll1)
xx <- seq(from=160, to=180, by=0.5)
yy <- seq(from=3, to=7, by=0.5)
zz <- outer(X=xx, Y=yy, ll1)
filled.contour(x=xx, y=yy, z=zz, nlevels=20,
               xlab=expression(mu), ylab=expression(sigma),
               color = topo.colors)

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

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. 

minusll <- function(x) -ll(x)
nlminb(objective=minusll, start=c(163, 3.4),
       lower=c(160, 3), upper=c(180, 7))
## $par
## [1] 170.338374   5.423529
## 
## $objective
## [1] 155.4842
## 
## $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 anterior 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 R 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.

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. 

require(MASS) 
res <- fitdistr(x=y, densfun='normal')
res
##       mean           sd     
##   170.3383794     5.4235271 
##  (  0.7670026) (  0.5423527)

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.4842 (df=2)

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

Ejemplo 5

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.3082767   2.8084910 
##  (0.2808491) (0.1985903)

2) 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.23978235   0.51987909 
##  (0.29620136) (0.07702892)

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.

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

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

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

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.

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

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

Ejemplo 6

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

n <- 1000
x1 <- runif(n=n, min=-5, max=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.

minusll <- function(theta, y, x1) {
  media <- theta[1] + theta[2] * x1 
  desvi <- theta[3]                  
  - 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.904603  3.079600  5.014184
## 
## $value
## [1] 3031.209
## 
## $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 \(\hat{\beta}_0 = -1.9046025\), \(\hat{\beta}_1 = 3.0796\) y \(\hat{\sigma} = 5.0141842\), se observa también que el valor de la máxima log-verosimilitud fue de -3031.2092104. Vemos entonces que el vector estimado está muy cerca del verdadero \(\boldsymbol{\theta}=(\beta_0=-2, \beta_1=3, \sigma=5)^\top\).

Cuando se usa optim es necesario decirle que inicie la búsqueda de \(\boldsymbol{\theta}\) a partir de un lugar. Por esa razón se usó par=c(0, 0, 1), esto significa que la búsqueda inicia en el tripleta \(\beta_0=0\), \(\beta_1=0\) y \(\sigma=1\).

6. Referencias