Objetivo:

Construir função na linguagem R para calcular a log-verossimilhança de uma amostra e obter numericamente as estimativas dos parâmetros que maximizam a verossimilhança da amostra.

1. Estimação de Máxima Verossimilhança:

A estimação de máxima verossimilhança é um método estatístico para estimar os parâmetros de um modelo. Na estimativa de máxima verossimilhança, os parâmetros são escolhidos de maneira a maximizar a probabilidade de que o modelo assumido resulte nos dados observados.

Isto implica que, para implementar obter a estimativa de máxima verossimilhança, devemos:

  1. Supor um modelo paramétrico subjacente aos dados observados.

  2. Considerando o modelo assumido, ser capaz encontrar o máximo da função de verossimilhança dos dados observados, (A busca do máximo pode ser efetuada matematica ou numericamente).

Na verdade, a função de verossimilhança será simplificada tomando seu logarítmo, obtendo assim a função de log-verossimilhança. O objetivo é maximizar a função de verossimilhança. Como o logaritmo cresce monotonicamente, se maximizarmos o log da verossimilhança, automaticamente maximizamos a verossimilhança (uma vez que \(f(b) > f(a)\) somente quando \(b > a\)). Há muito material acessível à respeito desse tema, que será profundamente estudado nas disciplinas de Inferência Estatística. Como uma leitura rápida sugerida uma leitura (em diagonal) dos posts: Beginner’s Guide To Maximum Likelihood Estimation e Maximum Likelihood For Dummies.

2. Função de Log-verossimilhança de uma Amostra Exponencial:

Aqui tratamos de um problema bem simples, assumindo que o modelo subjacente tenha um único parâmetro a ser estimado por meio da maximização da função de log-verossimilhança.

Seja uma amostra aleatória \(X_1, X_2, \dots, X_n\) de uma população exponencial com parâmetro desconhecido \(\lambda\). O objetivo do procedimento é estimar o valor desse parâmetro desconhecido.

A expressão função de verossimilhança dessa amostra observada \(x_1, x_2, \dots, x_n\), suposta oriunda de uma população exponecial, com parâmetro \(\lambda\) desconhecido, é:

\[\begin{align} L(\lambda| x_1, x_2, \dots, x_n) & = \prod_{i=1}^n f(\lambda; x_i) \\ & = \prod_{i=1}^n \lambda \operatorname{e}^{-\lambda x_i} \\ & = \lambda^n \left(\operatorname{e}^{-\lambda x_1} + \operatorname{e}^{-\lambda x_2} + \dots + \operatorname{e}^{-\lambda x_n} \right) \\ & = \lambda^n \operatorname{e}^{-\lambda \sum_{i=1}^n x_i} \end{align}\]

A função de log-verossimilhança (\(l(\lambda| x_1, x_2, \dots, x_n)\)) é obtida aplicando-se o \(\log\) na expressão acima, ou seja:

\[\begin{align} l(\lambda| x_1, x_2, \dots, x_n) & = \log(\lambda^n \operatorname{e}^{-\lambda \sum_{i=1}^n x_i}) \\ & = \log(\lambda^n) + \log(\operatorname{e}^{-\lambda \sum_{i=1}^n x_i}) \\ & = n \log(\lambda) - \lambda \sum_{i=1}^n x_i \end{align}\]

Assim, dada a amostra observada, a estimativa de máxima verosssimilhança do parâmetro da exponencial (\(\hat{\lambda}\)) será o valor de \(\lambda\) que maximiza \(l(\lambda| x_1, x_2, \dots, x_n)\), ou seja \(\hat{\lambda}=\arg \max l(\lambda| x_1, x_2, \dots, x_n)\). Podemos determinar matematicamente uma expressão do valor de \(\lambda\) que zera a derivada de \(l(\lambda| x_1, x_2, \dots, x_n)\). Pode-se provar que a expressão do estimador de máxima verossimilhaça do parâmetro \(\lambda\) é:

\[ \hat{\lambda} = \frac{n}{\sum_{i=1} ^n x_i} = \frac{1}{\bar{x}} \]

Entretanto é possível também encontrar a estimativa de máxima verossimilhança de uma amostra observada otimizando numericamente a função \(l(\lambda| x_1, x_2, \dots, x_n)\). Para tanto, é necessário construir um código em R para determinar o valor da log-verossimilhança em função de lambda. Consideremos que os dados observados estejam armazenados no vetor amst. Assim, a codificação abaixo pode ser um primeiro rascunho do código desejado.

logVeros.exp <- function(lambda){ 
    ## Função para cálculo da logverossimilhança da exponencial para 
    ## amst (amostra dada).

      # saída será impressa com 3 dígitos
    op <- options(digits = 3)  
    
      # valor da logverossimilhança da amostra dada (amst) para um valor de lambda              
    LV <- sum(log(dexp(rate = lambda, x = amst)))
      # Saída será imporessa com 3 dígitos
    on.exit(options(op), add = TRUE)

      # valor da saída
    print(LV)
## As configurações globais serão preservadas
}

O código está inadequado, pois não usa os dados observados como argumento de entrada da função, mas sim um objeto global amst. O resultado será impresso no console options(digits = 3), mas a execução da função não afetará a configuração global do usuário, pois esse comando está sendo executado com o comando on.exit(). Essa função é executada quando a função é encerrada, independentemente de ter sido gerado algum tipo de erro. Lembre-se também que ao executar o comando return() a função é automaticamente encerrada, por isso está sendo utilizado o comando print(). Verifique os posts How and when should I use on.exit? e On on.exit(), para conhecer alguns detalhes sobre o uso desse comando. O uso de configurações temporárias é muito comum ao elaborarmos funções, pricipalmente nas funções gráficas, já que muitas configurações são ajustadas globalmente (R using temporary options settings inside a function). Recomendo também a leitura dos posts Return value from R function e Difference of print and return in R? para entender a mecânica de funcionamento dos comandos print() e return(), que são muito utilizado nas saídas das funções.

Para verificar a função construída (logVeros.exp()) , vamos gerar uma amostra aleatória de uma exponencial padrão (\(\lambda = 1\)), armazenando-a no vetor amostra. Utilizei a semente usual (666) para fins de conferência de valores.

  # geração da amostra amst
set.seed(666)
amostra <- rexp(15)

Note que a média dessa amostra é \(1.215\). Vamos transferir essa amostra para o objeto amst (esse objeto global é usado inadequadamente pela função que criamos!),

      # criação de objeto global utilizado pela função.
amst <- amostra

O valor da log-verossimilhança para \(\lambda = 1.5\) é:

  # verossimilhança da taxa lambda = 1.5 para a amostra amst
logVeros.exp(1.5)
## [1] -21.3

Perceba que a saída no console foi com 3 dígitos, mas com a preservação da opção global da quantidade de dígitos:

      # a opção global de quantidade de dígitos foi preservada
options()$digits
## [1] 7

O ideal é codificar ligeiramente diferente a função, de maneira a não nos preocuparmos com a maneira que serão apresentados os valores no console. Afinal, queremos usar a função para: cálculos vetorizados, gráficos e procedimentos de otimização numérica. Dessa maneira, foram efetuadas ligeira modificações. Foi trocado o comando de saída de print() (imprime no console e retorna o valor) por return() (retorna o valor e finaliza a função). Foi também eliminada a opção de impressão dos valores no console. Dessa maneira, a função logL.exp() calcula também a log-verossimilhança de uma amostra proveniente de uma população exponencial com parâmetro desconhecido.

logL.exp <- function(lambda){
    ## Função para cálculo da logverossimilhança da exponencial para 
    ## amst (amostra observada deve ser armazenada nesse objeto fora da função).
  
        # valor da logverossimilhança da amostra dada (amst) para um valor de lambda              
    LV <- sum(log(dexp(rate = lambda, x = amst)))
        # valor da saída
      return(LV)
}

Essa função não é vetorizada. Para calcularmos o valor da log-verossimilhança para vários valores possíveis do parâmetro \(\lambda\) podemos usar o comando sapply()

   # Cálculo da verossimilhança para vários valores de lambda
sapply(c(0.1, 0.5, 0.95, 1, 1.1, 1.5, 2), logL.exp)
## [1] -36.36127 -19.50969 -18.08312 -18.22497 -18.61781 -21.25548 -26.05273

Para construir um gráfico da log-verossimilhança devemos ter uma função vetorizada ou podemos usar o comando sapply() diretamente na função gráfica. Vamos criar uma função logL.exp.vet(), vetorizada, para essa finalidade:

  # Função vetorizada para o cálculo da logverossimilhança para vários valores de lambda
logL.exp.vet <- function(x) sapply(x, logL.exp)

Dessa maneira, calculam-se as log-verossimilhança para valores possíveis do parâmetro \(\lambda\), como, por exemplo, valores no intervalo entre 0 e 2.

logL.exp.vet(c(0.1, 0.5, 0.95, 1, 1.1, 1.5, 2))
## [1] -36.36127 -19.50969 -18.08312 -18.22497 -18.61781 -21.25548 -26.05273

Epa. Nossa função não está legal teríamos de impedir que valores menor ou igual a zero sejam processados, porque \(\lambda\) não pode assumir esses valores. Cuidado, pois o espaço paramétrico é o conjunto \(\boldsymbol{\Theta} = \{\forall \lambda \in \mathbb{R}, \lambda > 0\}\). Observe também que o máximo está entre \(0.5\) e \(1.1\), assumindo que não haja máximos locais (sabemos que o parãmetro verdadeiro é \(1\)!). Inspecionando a função de log-verossimilhança no intervalo entre \(o\) e \(2\), obtemos:

  # curva da logverossimilhança para lambda entre 0 e 2, para amostra amst
curve(logL.exp.vet, from = 0, to = 2)

Nossa questão é determinar qual o valor de lambda que oferece a máxima verossimilhança para a amostra amst?

3. Otimização Numérica da Função de Log-verossimilhança:

No contexto de estimação de máxima verossimilhança, o termo otimização significa obter o valor do parãmetro que leva ao maior valor da função de verossimilhança ou a um valor tão grande quanto possível. As funções de otimização numérica incorporadas nao R Base são optimize() e optim(). optimize() é usada ao trabalhar com função com uma única variável (um único parâmetro, no contexto de estimação de máxima verossimilhança) e optim() é usada quando ao trabalhar com funções de várias variáveis (mais de um parâmetro, no contexto de estimação de máxima verossimilhança). Em nosso singelo exemplo, como temos de estimar um único parâmetro (a taxa \(\lambda\)) vamos utilizar a função optimize(). Essa função busca o valor do parâmetro que maximiza (ou minimiza) a função (f) no intervalo indicado (interval), do limite inferior para o superior. Sua sintaxe é:

optimize(f, interval, ..., lower = min(interval), upper = max(interval),
         maximum = FALSE,
         tol = .Machine$double.eps^0.25)

No nosso exemplo, buscaremos o valor de parâmetro, no intervalo \((0, 10)\) que maximiza (maximum = TRUE) a função de log-verossimilhança (logL.exp):

# Determinação de lambda com o máximo da verossimilhança 

  # Usando a função optimize()
obj.mono <- optimize(f = logL.exp, interval = c(0, 10), maximum = TRUE)
obj.mono
## $maximum
## [1] 0.8230463
## 
## $objective
## [1] -17.92114

A saída da função fornece o valor do parâmetro (obj.mono$maximum) que leva ao máximo da log-verossimilhança, ou da verossimilhança (obj.mono$objective), ou seja, logL.exp(obj.mono$maximum) = obj.mono$objective. Perceba também que a estimação numérica do valor do parâmetro, \(0.8230463\) (obj.mono$maximum), é bastante próximo \(0.8230467\), obtido da expressão \(\hat{\lambda} = \frac{1}{\bar{x}}\). Lembre-se que essa expressão é válida para obter a estimação de máxima verossimilhaça do parãmetro \(\lambda\) de qualquer amostra proveniente de população exponencial. O importante é que o verdadeiro (e desconhecido) valor do parâmetro \(\lambda\) é \(1\). Na curva abaixo, pode-se observar os valores estimados obtidos pelo procedimento de otimização numérica.

Por outro lado, em nosso exemplo (busca de um único parâmetro), poderíamos ainda ter usado a função optim(), que busca valores que maximizam funções multivariadas, tomando alguns cuidados importantes. Essa função também pode ser usada para otimização unidimensional. Ela recebe como entrada uma função a ser otimizada (f), um intervalo de valores para a variável independente (lower e upper)e um método de otimização (method). O formato geral da função optim() é:

optim(par, fn, lower = -Inf, upper = Inf, method = "Nelder-Mead", ...),

em que,

O resultado da função optim() é um objeto do tipo list, que contém as seguintes informações:

Os métodos de otimização disponíveis na função optim() são:

Consulte a documentação do R para obter mais informações sobre a função optim().

Por outro lado, no argumento control (lista de controle do parâmetro), é preciso usar fnscale = -1 como a escala global do procedimento, em que a otimização é efetuada em fn(par)/fnscale. Ao utilizar essa escala, o procedimento de otimização passa a buscar o(s) máximo(s) da função fn. Outro argumento que deve ser adotado é method, o qual precisa ser ajustado para "Brent", que é usado para problemas de otimização unidimensional usando optim(). Como argumentos adicionais para o presente problema (otimização unidimensional) é necessário especificar os limites lower e upper do intervalo de busca do parâmetro. (A especificação do intervalo de busca também é necessário nos problemas de otimização com restrição, com method = "L-BFGS-B".). O código para usar a função optim() para obter a estimativa de máxima verossimilhança de nossa amostra (amst) é dado abaixo:

obj.bi <- optim(par = 0, fn = logL.exp, method = "Brent",
            lower = 0, upper = 3, control = list(fnscale = -1))
obj.bi
## $par
## [1] 0.8230467
## 
## $value
## [1] -17.92114
## 
## $counts
## function gradient 
##       NA       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Perceba que foi utilizado o valor par = 0, como nossa aposta inicial no valor do parâmetro, Devem ser tomar cuidado com essa escolha, principalmente quando houver a possibilidade de máximos locais na região de busca do parâmetro. Nesses casos, é interessante verificar se a solução encontrada é robusta a diferentes valores iniciais para esse parâmetro. Outra especificação interessante seria ajustar o intervalo de busca (argumentos lower e upper) , ao intervalo entre o mínimo e o máximo da amostra. Sugiro uma leitura atenta dos posts How to use optim in R e optim() Function in R para verificar mais detalhes sobre uso da função optim() em otimização numérica de funções uni e multi-dimensionais. Lembre-se que grande parte da solução de problemas em Estatística passa pelo uso de procedimentos de otimização ou de integração numéricas.

Os valores encontrados são: \(0.8230467\) (usando a expressão \(1/\bar{x}\)); \(0.8230463\) (usando a função optimize()) e \(0.8230467\) (usando a função optim(()). O valor verdadeiro do parâmetro (usado na simulação da amostra amst) foi \(\lambda = 1\). Na prática esse valor é desconhecido.

4. Próximos Passos na Construção da Função:

Esse foi o primeiro passo na construção de uma função para obter numericamente a estimativa de máxima verossimilhança de uma amostra oriunda de uma população exponencial. Um primeiro aprimoramento é o vetor com os valores da amostra ser usado como argumento de entrada da função. Uma função montada dessa maneira poderia ser usada para simular \(N\) estimativas de máxima verossimilhança usando otimização numérica. O objetivo dessa simulação seria verificar a como os valores das estimativas estão distribuídos (em média o valor das \(N\) estimativas estão próximos do verdadeiro valor do parãmetro?. A dispersão em torno da média das estimativas muda com o tamanho das amostras? Poderíamos sugerir alguma distribuição para essas estimativas?). Outro aprimoramento seria montar uma função que aceite qualquer tipo de distribuição para a população. Há alguma outra sugestão? Compartilhem nesse fórum.

Referências: