Distribución Beta

Sea \(X_1, \ldots, X_n\) es una muestra aleatoria de una distribución \(Be(\alpha, \beta)\), con \(\alpha > 0\) y \(\beta > 0\) desconocidos. La función de verosimilitud de \(\alpha\) y \(\beta\) es

\[ L(\alpha, \beta\,\vert\,\mathbf{x}) = B(\alpha, \beta)^{-n}\prod_{j=1}^nx_j^{\alpha-1} (1-x_j)^{\beta-1}. \]

La log-verosimilitud es

\[ \ell(\alpha, \beta\,\vert\,\mathbf{x}) = -n\log B(\alpha, \beta) + (\alpha-1)\sum_{j=1}^n\log{x_j} + (\beta-1)\sum_{j=1}^n\log{(1-x_j)} \]

Al intentar maximizar analíticamente nos encontramos con el problema de derivar \(\log B(\alpha, \beta)\).

Los estimadores de máxima verosimilitud no tienen una expresión análitica, pero pueden ser calculados numéricamente para cada caso particular. Para estimar numéricamente en R, utilizamos las funciones nlm y constrOptim.

nlm

La función nlm se utiliza para minimizar. Requiere como mínimo dos parámetros para operar:

  • La funición objetivo.
  • Un vector de valores iniciales.

Consultando la ayuda ??nlm se pueden ver más detalles acerca de esta función.

constrOptim

La funicón constrOptim también se utiliza para minimizar. La principal diferencia es que permite minimzar con restricciones lineales de desigualdad. Esto es de particular interés puesto que en la mayoría de los casos tenemos restricciones en los valores posibles de los parámetros. Por ejemplo, en el caso de la distribución \(Beta\), \(\alpha > 0\) y \(\beta > 0\). Esta función requiere más argumentos para operar:

  • Un vector de valores iniciales.
  • La función objetivo.
  • El gradiente de la función objetivo (aunque se puede delcarar como NULL).
  • Una matriz ui y un vector ci que determinan las restricciones de la forma

\[ \mathbf{U}_i \boldsymbol{\theta} \geq \mathbf{c}_i \]

donde las desigualdades aplican componente a componente.

En nuestro caso, la función objetivo es el negativo de la log-verosimilitud (negativo porque nosotros queremos maximizar). Como valores iniciales podemos tomar los estimadores de momentos. Las restricciones dependerán de cada caso particular. Por ejemplo, si \(\boldsymbol{\theta}^T = (\alpha, \beta)\) y las resticciones son \(\alpha >0\) y \(\beta >0\), entonces \[ \mathbf{U}_i = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) \qquad \text{y} \qquad \mathbf{c}_i = \left( \begin{array}{c} 0 \\ 0 \end{array} \right). \]

Ejemplo

Generamos una muestra de tamaño \(n = 147\) de una distribución \(Be(3, 5)\). El objetivo es estimar por máxima verosimilitud.

set.seed(1010)
muestra <- rbeta(147, 3, 5)

Como valores iniciales podemos tomar los estimadores de momentos, que son:

\[ \hat{\alpha} = \bar{X}_n\left(\frac{\bar{X}_n(1-\bar{X}_n)}{S^2_n}-1\right) \qquad \text{y} \qquad \hat{\beta} = (1-\bar{X}_n) \left(\frac{\bar{X}_n(1-\bar{X}_n)}{S^2_n}-1\right) \]

Primero calculamos los estadísticos \(\bar{X}_n\) y \(S^2_n\)

xb <- mean(muestra)
s2 <- mean((muestra-mean(muestra))^2)
xb; s2
## [1] 0.3574178
## [1] 0.02426896

Ahora calculamos los estimadores de momentos

a_mom <- xb * (xb*(1-xb)/s2 - 1)
b_mom <- (1-xb) * (xb*(1-xb)/s2 - 1)

Definimos la log-verosimilitud

log.vero <- function(theta) - sum(log(dbeta(muestra, theta[1], theta[2])))

Maximizamos con la función nlm

Lmax.1 <- nlm(log.vero, c(a_mom, b_mom))
Lmax.1
## $minimum
## [1] -69.27177
## 
## $estimate
## [1] 3.067443 5.492221
## 
## $gradient
## [1]  4.169522e-08 -5.174903e-08
## 
## $code
## [1] 1
## 
## $iterations
## [1] 5

Maximizamos con la función constrOptim

Lmax.2 <- constrOptim(theta = c(a_mom, b_mom), f = log.vero, grad = NULL,
                      ui = diag(1, 2), ci = c(0, 0))
Lmax.2
## $par
## [1] 3.067274 5.492083
## 
## $value
## [1] -69.27177
## 
## $counts
## function gradient 
##       43       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $outer.iterations
## [1] 2
## 
## $barrier.value
## [1] -0.0004234404

Comparamos resultados

alpha beta
Momentos 3.025021 5.438522
MV nlm 3.067443 5.492221
MV constrOptim 3.067274 5.492083

Distribución Gamma

Sea \(X_1, \ldots, X_n\) es una muestra aleatoria de una distribución \(Ga(\alpha, \lambda)\), con \(\alpha > 0\) y \(\beta > 0\) desconocidos. La función de verosimilitud de \(\alpha\) y \(\lambda\) es

\[ L(\alpha, \beta\,\vert\,\mathbf{x}) = \Gamma(\alpha)^{-n}\lambda^{n\alpha} \left(\prod_{j=1}^nx_j\right)^{\alpha-1} e^{-\lambda \sum_{j=1}^nx_j} . \]

La log-verosimilitud es

\[ \ell(\alpha, \beta\,\vert\,\mathbf{x}) = -n\log \Gamma(\alpha) + n\alpha\log\lambda + (\alpha-1)\sum_{j=1}^n\log{x_j} - \lambda\sum_{j=1}^nx_j \]

Al intentar maximizar analíticamente nos encontramos con el problema de derivar \(\log \Gamma(\alpha)\). Podemos maximizar numéricamente.

Ejemplo

Generamos una muestra de tamaño \(n = 213\) de una distribución \(Ga(3, 5)\). El objetivo es estimar por máxima verosimilitud.

xb <- mean(muestra2) s2 <- mean((muestra-mean(muestra2))^2)

set.seed(1010)
muestra2 <- rgamma(213, 3, 5)

Como valores iniciales podemos tomar los estimadores de momentos, que son:

\[ \hat{\alpha} = \frac{\bar{X}^2_n}{S^2_n} \qquad \text{y} \qquad \hat{\lambda} = \frac{\bar{X}_n}{S^2_n}. \]

Primero calculamos los estadísticos \(\bar{X}_n\) y \(S^2_n\)

xb <- mean(muestra2)
s2 <- mean((muestra2-mean(muestra2))^2)
xb; s2
## [1] 0.6308564
## [1] 0.1565441

Ahora calculamos los estimadores de momentos

al_mom <- xb^2 / s2
la_mom <- xb / s2

Definimos la log-verosimilitud

log.vero2 <- function(theta) - sum(log(dgamma(muestra2, theta[1], theta[2])))

Maximizamos con la función nlm

Lmax2.1 <- nlm(log.vero2, c(al_mom, la_mom))
Lmax2.1
## $minimum
## [1] 65.20799
## 
## $estimate
## [1] 2.857121 4.528954
## 
## $gradient
## [1] -7.883532e-06  7.599699e-06
## 
## $code
## [1] 1
## 
## $iterations
## [1] 7

Maximizamos con la función constrOptim

Lmax2.2 <- constrOptim(theta = c(al_mom, la_mom), f = log.vero2, grad = NULL,
                      ui = diag(1, 2), ci = c(0, 0))
Lmax2.2
## $par
## [1] 2.857211 4.528986
## 
## $value
## [1] 65.20799
## 
## $counts
## function gradient 
##       59       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $outer.iterations
## [1] 2
## 
## $barrier.value
## [1] -0.0002453668

Comparamos resultados

alpha lambda
Momentos 2.542286 4.029897
MV nlm 2.857121 4.528954
MV constrOptim 2.857211 4.528986