A estimação por máxima verossimilhança possui

Há vários pacotes que ajudam a implementar a estimação por máxima verossimilhança. Neste post vou me ater apenas a dois pacotes: o optimx e o maxLik. O primeiro deles agrega funções de otimização de diversos outros pacotes numa sintaxe unificada centrada em algumas poucas funções. O último é feito especificamente para estimação de máxima verossimilhança então traz algumas comodidades como a estimação automática de erros-padrão.

library(maxLik)
library(optimx)
# Para reproduzir os resultados
set.seed(33)

Exemplo com distribuição Poisson

A ideia da estimação por máxima verossimilhança é de encontrar os parâmetros que maximizam a probabilidade de que um certo conjunto de dados sejam gerados por algum processo específico. Agora, em português: imagine que temos uma amostra aleatória; não sabemos qual distribuição os dados seguem, mas vamos supor que eles seguem alguma distribuição específica, por exemplo, uma Poisson. O formato da Poisson depende do parâmetro \(\lambda\) e a ideia então é de encontrar o valor de \(\lambda\) que melhor aproxima a distribuição empírica/observada dos dados. Isto é, dado uma amostra aleatória \(x_{1}, x_{2}, \dots , x_{n}\) buscamos o valor de \(\lambda\) que maximiza a probabilidade de que os dados tenham sido gerados por uma distribuição de Poisson.

A Poisson é uma distribuição discreta parametrizada por um valor \(\lambda > 0\). O formato da distribuição varia com o valor de \(\lambda\). Lembrando que a função de distribuição da Poisson é dada pela expressão:

\[ f(k, \lambda) = \text{Pr}(X = k) = \frac{\lambda^{k}e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \dots \]

O código abaixo faz o gráfico da função acima para diferentes valores de \(\lambda\).

Para montar a função de log-verossimilhança começamos com a função de verossimilhança de uma observação \(i\) qualquer e depois montamos a distribuição conjunta dos dados. Lembrando que a função de verossimilhança da Poisson para uma observação \(i\) é dada por:

\[ f(x_{i}, \lambda) = \frac{\lambda^{x_{i}}e^{-\lambda}}{x_{i}!} \]

Supondo que a amostra é de fato aleatória, a distribuição conjunta dos dados é simplesmente o produtório das distribuições individuais, isto é:

\[\begin{align} \prod_{k = 1}^{n} f(x_{k}, \lambda) & = \prod_{k = 1}^{n} \left( \frac{\lambda^{x_{k}}e^{-\lambda}}{x_{k}!} \right) \\ & = \frac{\lambda^{\sum x_{k}} e^{-n\lambda}}{\prod_{k = 1}^{n}(x_{k}!)} \end{align}\]

Aplicando log chegamos na função de log-verossimilhança:

\[\begin{align} \mathcal{L}(\lambda ; x_{1}, \dots , x_{n}) & = \text{ln}(\prod_{k = 1}^{n} f(x_{k}, \lambda)) \\ & = \text{ln}(\lambda^{\sum x_{k}}) + \text{ln}(e^{-n\lambda}) - \text{ln}(\prod_{k = 1}^{n}(x_{k}!)) \\ & = \text{ln}(\lambda)\sum_{k = 1}^{n}x_{k} - n\lambda - \sum_{k = 1}^{n} \text{ln}(x_{k}!) \end{align}\]

Dada uma amostra aleatória \(x_{1}, x_{2}, \dots , x_{n}\) buscamos o valor de \(\lambda\) que maximiza a função acima. Isto é, temos o seguinte problema de otimização:

\[ \underset{\lambda > 0}{\text{Max}} \quad \text{ln}(\lambda)\sum_{k = 1}^{n}x_{k} - n\lambda - \sum_{k = 1}^{n} \text{ln}(x_{k}!) \] onde os valores de \(x_{k}\) são os nossos dados observados. Para implementar este procedimento no R primeiro montamos a função log-verossimilhança no código abaixo. Em seguida vamos usar a função base optim que resolve problemas de minimização. Como a função optim serve para minimizar funções precisamos implementar o negativo da função de log-verossimilhança (lembrando que maximizar \(f(x)\) é equivalente à minimizar \(-f(x)\)).

ll_pois <- function(x, theta) {
    n <- length(x)
    ll <- log(theta) * sum(x) - n * theta - sum(log(factorial(x)))
    return(-ll)
}

Vamos simular 1000 observações \(x_{i} \sim \text{Poisson}(\lambda = 5)\).

amostra <- rpois(1000, lambda = 5)

Podemos tornar mais claro o procedimento de estimação olhando para o gráfico da função. O código abaixo plota o gráfico da função de log-verossimilhança para valores de \(\lambda\) no intervalo \([0, 20]\). Note que a função parece atingir um máximo em torno de \(5\), o valor verdadeiro do parâmetro.

eixo_lambda <- seq(0, 20, .001)
plot(eixo_lambda, -ll_pois(amostra, eixo_lambda), type = "l",
     xlab = bquote(lambda),
     ylab = bquote("L("~lambda~";"~x[1]~", ... ,"~x[n]~")"),
     main = "Função log-verossimilhança da Poisson")
abline(v = 5, col = "indianred")

Para estimar \(\lambda\) usamos a função optim. É preciso definir algum valor inicial para que o otimizador numérico comece a procurar pelo máximo global. Em casos mais simples, de funções globalmente côncavas, esta escolha não apresenta grandes consequências. Em casos mais complicados, o resultado da estimação pode ser bastante sensível à escolha de valor inicial porque o otimizador pode cair em máximos/mínimos locais. No final do post discuto brevemente algumas soluções. Neste exemplo escolho arbitrariamente começar no valor \(1\). Neste caso também poderíamos usar o fato de que a esperança da Poisson é igual a \(\lambda\), logo, pela Lei dos Grandes Números, a média amostral já deve estar próxima do verdadeiro valor de \(\lambda\). Assim, poderíamos também ter escolhido mean(amostra) como valor inicial.

fit <- optim(par = 1, fn = ll_pois, x = amostra,
             method = "BFGS", hessian = TRUE)
fit
## $par
## [1] 4.938
## 
## $value
## [1] 2202.837
## 
## $counts
## function gradient 
##       34        9 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
##          [,1]
## [1,] 202.5112

A saída da função traz várias informações técnicas sobre a otimização numérica. A estimativa pontual computada foi de \(4.938\) - bastante próxima do verdadeiro valor do parâmetro. Usando a função str podemos observar a estrutura do objeto fit criado acima. As principais informações que podemos extrair são as estimativas dos parâmetros fit$par e a hessiana calculada nestes parâmetros fit$hessian.

str(fit)
## List of 6
##  $ par        : num 4.94
##  $ value      : num 2203
##  $ counts     : Named int [1:2] 34 9
##   ..- attr(*, "names")= chr [1:2] "function" "gradient"
##  $ convergence: int 0
##  $ message    : NULL
##  $ hessian    : num [1, 1] 203

Usando a estimativa da hessiana podemos computar o erro-padrão da estimativa. Lembre que a variância assintótica do estimador de máxima verossimilhança é o inverso da matriz de informação de Fisher que pode ser expressa como o negativo do valor esperado da hessiana. Isto é, podemos encontrar a variância assintótica calculando o inverso da hessiana (como fizemos a minimização do negativo da função log-verossimilhança não precisamos calcular o negativo da hessiana).

(ep <- sqrt(1/fit$hess))
##           [,1]
## [1,] 0.0702709

Com o erro-padrão podemos calcular intervalos de confiança e estatística t.

(ic <- c(fit$par - 1.96 * ep, fit$par + 1.96 * ep))
## [1] 4.800269 5.075731
(est_t <- (fit$par - 5) / ep)
##            [,1]
## [1,] -0.8823045

Usando o pacote optimx

A função optim já é bastante antiga e um novo pacote, chamado optimx, foi criado. A ideia do pacote é de agregar várias funções de otimização que estavam espalhadas em diversos pacotes diferentes. As principais funções do pacote são optimx e optimr. Mais informações sobre o pacote podem ser encontradas aqui. A sintaxe das funções é muito similar à sintaxe original do optim. O código abaixo faz o mesmo procedimento de estimação que o acima.

summary(fit <- optimx(par = 1, fn = ll_pois, x = amostra))
##                 p1    value fevals gevals niter convcode kkt1 kkt2 xtime
## Nelder-Mead 4.9375 2202.837     32     NA    NA        0   NA   NA  0.02
## BFGS        4.9380 2202.837     34      9    NA        0   NA   NA  0.00

Uma das principais vantagens do optimx é a possibilidade de usar vários métodos de otimização numérica numa mesma função.

(fit <- optimx(par = 1, fn = ll_pois, x = amostra,
               method = c("nlm", "BFGS", "Rcgmin", "nlminb")))
##              p1    value fevals gevals niter convcode kkt1 kkt2 xtime
## nlm    4.937998 2202.837     NA     NA     8        0   NA   NA  0.00
## BFGS   4.938000 2202.837     34      9    NA        0   NA   NA  0.01
## Rcgmin 4.938000 2202.837     15     10    NA        0   NA   NA  0.00
## nlminb 4.938000 2202.837     10     12     9        0   NA   NA  0.00

Como este exemplo é bastante simples os diferentes métodos parecem convergir para valores muito parecidos.

Usando o pacote maxLik

A função maxLik (do pacote homônimo) traz algumas comodidades: primeiro ela maximiza as funções de log-verossimilhança, ou seja, não é preciso montar a função com sinal de menos como fizemos acima; segundo, ela já calcula erros-padrão e estatísticas-t dos coeficientes estimados. Além disso, ela também facilita a implementação de gradientes e hessianas analíticos e conta com métodos de otimização bastante populares como o BHHH. Mais detalhes sobre a função e o pacote podem ser encontradas aqui.

Para usar a função precisamos primeiro reescrever a função log-verossimilhança, pois agora não precisamos mais buscar o negativo da função. Como o R já vem com as funções de densidade de várias distribuições podemos tornar o código mais enxuto usando o dpois que implementa a função densidade da Poisson. O argumento log = TRUE retorna as probabilidades \(p\) como \(log(p)\).

ll_pois <- function(x, theta) {
    ll <- dpois(x, theta, log = TRUE)
    return(sum(ll))
}
summary(fit <- maxLik(ll_pois, start = 1, x = amostra))
## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 7 iterations
## Return code 1: gradient close to zero
## Log-Likelihood: -2202.837 
## 1  free parameters
## Estimates:
##      Estimate Std. error t value Pr(> t)    
## [1,]   4.9380     0.0703   70.25  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Podemos implementar manualmente o gradiente e a hessiana da função. Neste caso, a estimativa do parâmetro continua a mesma mas o erro-padrão diminui um pouco. Note que também podemos fornecer estas informações para a função optimx. Derivando a função de log-verossimilhança:

\[\begin{align} \frac{\partial \mathcal{L}(\lambda ; x)}{\partial\lambda} & = \frac{1}{\lambda}\sum_{k = 1}^{n}x_{k} - n \\ \frac{\partial^2 \mathcal{L}(\lambda ; x)}{\partial\lambda^2} & = -\frac{1}{\lambda^2}\sum_{k = 1}^{n}x_{k} \end{align}\]

O código abaixo implementa o gradiente e a hessiana e faz a estimação.

grad_pois <- function(x, theta) {
    (1 / theta) * sum(x) - length(x)
}
hess_pois <- function(x, theta) {
    -(1 / theta^2) * sum(x)
}
summary(fit2 <- maxLik(ll_pois, grad = grad_pois, hess = hess_pois,
                       start = 1, x = amostra))
## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 7 iterations
## Return code 1: gradient close to zero
## Log-Likelihood: -2202.837 
## 1  free parameters
## Estimates:
##      Estimate Std. error t value Pr(> t)    
## [1,]  4.93800    0.07027   70.27  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Exemplo com a distribuição normal

A distribuição normal tem dois parâmetros: \(\mu\) e \(\sigma\). Lembrando que o primeiro indica a média e o segundo a desvio-padrão. A função de densidade de probabilidade é dada por:

\[ f(x, \theta) = \frac{1}{\sqrt{2\pi\sigma^{2}}}\text{exp}\left(\frac{-(x - \mu)^2}{2\sigma^{2}}\right) \]

onde \(\theta = (\mu, \sigma)\). O gráfico abaixo mostra como o formato da função varia conforme o valor destes parâmetros.

Lembrando que a função de distribuição de probabilidade da normal para uma observação \(i\)

\[ f(x_{i}, \theta) = \frac{1}{\sqrt{2\pi\sigma^{2}}}e^{\frac{-(x_{i} - \mu)^2}{2\sigma^{2}}} \]

Fazendo o produtório da expressão acima para cada \(i\)

\[ \prod_{i = 1}^{N}f(x_{i}, \theta) = (2\pi\sigma^{2})^{-\frac{N}{2}}\text{exp}\left( -\frac{1}{2\sigma^{2}}\prod_{i = 1}^{N}(x_{i} - \mu)^2\right) \]

e agora passando log, temos:

\[ L(x_{i}, \theta) = -\frac{N}{2}\text{ln}(2\pi) -N\text{ln}(\sigma) -\frac{1}{2\sigma^{2}}\sum_{i = 1}^{N}(x_{i} - \mu)^{2} \]

Montamos a função log-verossimilhança usando a função dnorm

ll_norm <- function(theta) {
    ll <- dnorm(x, theta[1], theta[2], log = TRUE)
    -sum(ll)
}

Vamos simular uma amostra aleatória \(x_{1}, x_{2}, \dots, x_{1000}\) onde cada \(x_{i}\) segue uma distribuição normal com média \(2\) e desvio-padrão \(3\) (i.e., variância \(9\)).

Primeiro vamos usar a função optimx

# warnings por causa dos valores negativos no log
x <- rnorm(n = 1000, mean = 2, sd = 3)
(fit <- optimx(par = c(1, 1), fn = ll_norm, method = "BFGS"))
##            p1       p2    value fevals gevals niter convcode kkt1 kkt2
## BFGS 2.138082 3.058112 2536.736     48     20    NA        0 TRUE TRUE
##      xtime
## BFGS  0.02

Note que a função retorna mensagens de erro indicando que a função retornou NaNs. Isto acontece porque o otimizador experimenta valores não-positivos para \(\sigma\) e isto não é admissível pois \(\text{ln}(\sigma)\) não é definido para \(\sigma < 0\). Além disso, \(\sigma\) não pode ser igual a zero pois ele aparece no denominador do último termo à direita da expressão da log-verossimilhança. Intuitivamente isto é bastante óbvio: \(\sigma\) representa o desvio-padrão da distribuição e \(\sigma^{2}\): não podemos ter valores negativos ou nulos para esta expressão.

Podemos restringir os valores dos parâmetros para evitar as mensagens de erro, mas, na prática, isto não costuma fazer diferença no resultado final da estimação.

fit <- optimx(par = c(1, 1), fn = ll_norm, method = "L-BFGS-B",
                 upper = c(Inf, Inf), lower = c(-Inf, 0))

Propriedades dos estimadores de MV

A teoria dos estimadores de máxima verossimilhança nos garante que eles são consistentes (i.e. que eles aproximam o valor verdadeiro dos parâmetros) e normalmente assintóticos (a distribuição assintótica segue uma distribuição normal) desde que algumas condições de regularidade sejam atentdidas. Podemos visualizar isto através de um experimento simples: simulamos 5000 amostras aleatórias de tamanho 1000 seguindo uma distribuição \(N(2, 3)\); computamos as estimativas para \(\mu\) e \(\sigma\) e suas respectivas variâncias assintóticas e depois analisamos suas propriedades. O código abaixo implementa exatamente este experimento. Note que a matriz de informação de Fisher é aproximada pela hessiana.

r <- 5000; n <- 1000
estimativas <- matrix(ncol = 4, nrow = r)
for(i in 1:r) {
    x <- rnorm(n = n, mean = 2, sd = 3)
    fit <- optimr(par = c(1, 1), fn = ll_norm, method = "BFGS", hessian = TRUE)
    # Guarda o valor estimado do parâmetro
    estimativas[i, 1:2] <- fit$par
    estimativas[i, 3:4] <- diag(n * solve(fit$hess))
}

A consistência dos estimadores \(\hat{\theta}_{MV}\) significa que eles aproximam os valores verdadeiros do parâmetros \(\theta_{0}\) à medida que aumenta o tamanho da amostra. Isto é, se tivermos uma amostra grande então podemos ter confiança de que nossos estimadores estão muito próximos dos valores verdadeiros dos parâmetros.

\[ \hat{\theta}_{MV} \underset{N \to \infty}{\rightarrow} \theta_{0} \] O código abaixo calcula a média das estimativas para cada parâmetro - lembrando que \(\mu_{0} = 2\) e que \(\sigma_{0} = 3\). Além disso, o histograma das estimativas mostra como as estimativas concentram-se em torno do valor verdadeiro do parâmetro (indicado pela linha vertical).

# Consistência dos estimadores de MV
mu <- estimativas[, 1]; sigma <- estimativas[, 2]
mean(mu)
## [1] 2.00088
mean(sigma)
## [1] 2.997335
hist(mu, main = bquote("Estimativas para "~~mu), freq = FALSE, xlim = c(1.5, 2.5))
abline(v = 2, col = "indianred")

hist(sigma, main = bquote("Estimativas para "~~sigma), freq = FALSE, xlim = c(2.7, 3.3))
abline(v = 3, col = "indianred")

Dizemos que os estimadores de máxima verossimilhança são normalmente assintóticos porque a sua distribuição assintótica segue uma normal padrão. Especificamente, temos que:

\[ z_{\theta} = \sqrt{N}\frac{\hat{\theta}_{MV} - \theta}{\sqrt{\text{V}_{ast}}} \to \mathbb{N}(0, 1) \]

onde \(\text{V}_{ast}\) é a variância assintótica do estimador. No loop acima usamos o fato que a matriz de informação de Fisher pode ser estimada pela hessiana. O código abaixo calcula a expressão acima para os dois parâmetros e apresenta o histograma dos dados com uma normal padrão superimposta.

# Normalidade assintótica

# Centra a estimativa
mu_centrado <- estimativas[, 1] - 2 
sigma_centrado <- estimativas[, 2] - 3
# Computa z_mu z_sigma
mu_normalizado <- sqrt(n) * mu_centrado / sqrt(estimativas[, 3])
sigma_normalizado <- sqrt(n) * sigma_centrado / sqrt(estimativas[, 4])
# Eixo x
grid_x <- seq(-3, 3, 0.01)

hist(mu_normalizado, main = bquote("Histograma de 5000 simulações para z"[mu]),
     freq = FALSE, xlim = c(-3, 3), breaks = 30,
     xlab = bquote("z"[mu]), ylab = "Densidade")
lines(grid_x, dnorm(grid_x, mean = 0, sd = 1), col = "dodgerblue")
legend("topright", lty = 1, col = "dodgerblue", legend = "Normal (0, 1)")

hist(sigma_normalizado, main = bquote("Histograma de 5000 simulações para z"[sigma]),
     freq = FALSE, xlim =c(-3, 3), breaks = 30,
     xlab = bquote("z"[sigma]), ylab = "Densidade")
lines(grid_x, dnorm(grid_x, mean = 0, sd = 1), col = "dodgerblue")
legend("topright", lty = 1, col = "dodgerblue", legend = "Normal (0, 1)")

Escolha de valores inicias

O exemplo abaixo mostra o caso em que a escolha de valores iniciais imprópria leva a estimativas muito ruins.

# sensível a escolha de valores inicias
x <- rnorm(n = 1000, mean = 15, sd = 4)
(fit <- optim(par = c(1, 1), fn = ll_norm, method = "BFGS", hessian = TRUE))
## $par
## [1] 617.9780 962.9146
## 
## $value
## [1] 7985.067
## 
## $counts
## function gradient 
##      107      100 
## 
## $convergence
## [1] 1
## 
## $message
## NULL
## 
## $hessian
##              [,1]          [,2]
## [1,]  0.001078433 -0.0013512818
## [2,] -0.001351282  0.0001909939

Note que as estimativas estão muito distantes dos valores corretos \(\mu = 15\) e \(\sigma = 4\). Uma das soluções, já mencionada acima, é de usar os momentos da distribuição como valores iniciais.

O código abaixo usa os momentos empíricos como valores inicias para \(\mu\) e \(\sigma\):

\[\begin{align} \mu_{inicial} & = \frac{1}{n}\sum_{i = 1}^{n}x_{i} \\ \sigma_{inicial} & = \sqrt{\frac{1}{n} \sum_{i = 1}^{n} (x_{i} - \mu_{inicial})^2} \end{align}\]

(chute_inicial <- c(mean(x), sqrt(var(x))))
## [1] 14.859702  3.930849
(est <- optimx(par = chute_inicial, fn = ll_norm))
##                   p1       p2    value fevals gevals niter convcode kkt1
## Nelder-Mead 14.85997 3.929097 2787.294     47     NA    NA        0 TRUE
## BFGS        14.85970 3.928884 2787.294     19      2    NA        0 TRUE
##             kkt2 xtime
## Nelder-Mead TRUE  0.00
## BFGS        TRUE  0.02

Agora as estimativas estão muito melhores. Outra opção é experimentar com otimizadores diferentes. Aqui a função optimx se prova bastante conveniente pois admite uma grande variedade de métodos de otimizãção. Mais detalhes sobre os métodos podem ser encontrados na página de ajuda da função ?optimx.

# Usando outros métodos numéricos
optimx(par = c(1, 1), fn = ll_norm, hessian = TRUE,
       method = c("nlminb", "bobyqa", "Nelder-Mead", "BFGS"))
##                    p1         p2    value fevals gevals niter convcode
## nlminb       14.85970   3.928883 2787.294     23     47    19        0
## bobyqa       15.20015   8.993214 3211.554    108     NA    NA        0
## Nelder-Mead  14.85571   3.929621 2787.294     83     NA    NA        0
## BFGS        617.97795 962.914587 7985.067    107    100    NA        1
##              kkt1  kkt2 xtime
## nlminb       TRUE  TRUE  0.00
## bobyqa      FALSE FALSE  0.03
## Nelder-Mead  TRUE  TRUE  0.02
## BFGS         TRUE FALSE  0.03

Aplicações

  • Estimando um modelo de mínimos quadrados por máxima verossimilhança
  • Estimando um logit por máxima verossimilhança
  • Aproximando distribuições por máxima verossimilhança