Verossimilhança da Poisson

Dislaimer: eu tenho a formação em estatística de uma batata, não me leve muito a sério

A distribuição de Poisson descreve a probabilidade de que \(k\) eventos discretos ocorram em um espaço ou período de tempo em que \(\lambda\) eventos eram esperados. A densidade é:

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

Existe um bom número de situações em que essa distribuição pode não ser adequada para modelar o que se propõe, mas vamos passar por cima disso. Duas propriedades que eu acho divertidinhas:

• O primeiro e o segundo momento são iguais, valem \(\lambda\).
• Seja \(X\) um conjunto de variáveis aleatórias onde vale que \(X_i \sim \text{Pois}(\lambda_i)\), então\(\sum_{i=1}^n X_i \sim \text{Pois}(\sum_{i=1}^n \lambda_i)\). A soma de várias Poisson, cada uma com \(\lambda_i\) gera uma Poisson cujo parâmetro é a soma dos \(\lambda_i\) individuais.

A pergunta é muito simples. Dado que tenho várias observações \(k_1, k_2, ..., k_n\) que suponho serem tiradas de um processo i.i.d de uma Poisson com parâmetro \(\lambda\), como estimo \(\lambda\)?

Note que para cada observação é possível calcular \(P(k_i \, | \, \lambda)\). Como estamos falando de um processo i.i.d a probabilidade de que todas as observações \(k_i\) tenham sido tiradas de uma mesma Poisson com parâmetro \(\lambda\) é apenas o produto:

\[P(k_1, k_2, ..., k_n \, | \, \lambda) = \prod_{i=1}^n \frac{\lambda^{k_i} e^{-\lambda}}{k_i!}\] Se você tomar as \(n\) observações como dadas então temos uma curva relacionando cada potencial \(\lambda\) à probabilidade de que os dados foram coletadas de uma Poisson com este parâmetro. Essa é a verossimilhança. Faz sentido então escolher o \(\lambda^*\) que maximiza a verossimilhança. Os dados indicam que aí é onde está o nosso melhor chute educado sobre a verdadeira taxa de ocorrência do fenômeno.

\[L(\lambda) = \prod_{i=1}^n \frac{\lambda^{k_i} e^{-\lambda}}{k_i!}\]

Precisamos otimizar \(L(\lambda)\), que tem uma forma funcional meio estranha. Uma boa ideia é otimizar o seu logaritmo, que chamaremos de \(l(\lambda)\). Se você fez inferência estatística na faculdade provavelmente já fez essa conta aqui (ou foi cobrado fazer):

\[l(\lambda) = \sum_{i+1}^n k_i\log{\lambda} - \lambda - \log{k_i!}\]

\[\frac{d l}{d \lambda} = \sum_{i+1}^n \frac{k_i}{\lambda} - 1 = 0 \] \[\sum_{i+1}^n \frac{k_i}{\lambda} = n \] \[\lambda^* = \frac{1}{n} \sum_{i+1}^n k_i \]

Esse trabalho todo para chegar na média. É isso, o melhor chute para a taxa de ocorrência verdadeira é a média das taxas de ocorrência observadas. Faz sentido… Agora vamos ver um pouco da magia acontecendo. Primeiro simular uma amostra com um \(\lambda\) positivo aleatório que precisaremos descobrir.

lambda_original <- runif(1, 25, 100)

amostra <- rpois(1000, lambda_original)

Definir a função de fatorial porque eu nunca entendi onde está a fatorial implementada em R:

fatorial <- function(n) {
  
  if(n == 0L | n == 1L) return(1)
  if(n > 1L)            return(n*fatorial(n - 1))
  
}

fatorial(4) # teste

## [1] 24

A agora a brincar:

(simulacao <- tibble(lambda = seq(25, 100, by = 0.1)) %>%
  mutate(L = map_dbl(
    lambda, 
    function(.x) sum(amostra*log(.x) - .x - log(map_dbl(amostra, fatorial)))
    )))

## # A tibble: 751 x 2
##    lambda       L
##     <dbl>   <dbl>
##  1   25   -17939.
##  2   25.1 -17814.
##  3   25.2 -17689.
##  4   25.3 -17566.
##  5   25.4 -17443.
##  6   25.5 -17321.
##  7   25.6 -17200.
##  8   25.7 -17080.
##  9   25.8 -16961.
## 10   25.9 -16843.
## # … with 741 more rows

(results <- simulacao %>%
  filter(L == max(L)))

## # A tibble: 1 x 2
##   lambda      L
##    <dbl>  <dbl>
## 1   56.5 -3410.

Já o \(\lambda\) verdadeiro era:

lambda_original

## [1] 56.47351

O que dá um erro absoluto de:

(simulacao %>%
  filter(L == max(L)) %>%
  pull(lambda) -
  lambda_original %>%
  abs() ->
  erro_abs)

## [1] 0.02648981

100 * erro_abs / lambda_original %>%
  round(2) # erro percentual

## [1] 0.04690953

Algo próximo de 1,13%

A curva que nós calculamos, por sinal, foi essa aqui:

simulacao %>%
  ggplot(aes(x = lambda, y = L)) +
  geom_line(size = 1.2, col = "red") +
  geom_hline(yintercept = pull(results, L)) +
  geom_vline(xintercept = pull(results, lambda)) +
  labs(title = "Curva da Log-Verossimilhança calculada a partir da amostra",
       x = "Lambda",
       y = "l(.)") +
  theme_minimal()