1) EMVs para \(\alpha\) e \(\beta\) da Distribuição Gama

A Função Densidade de Probabilidade (FDP) de uma distribuição Gama(\(\alpha, \beta\)) é dada por:

\[ f(y; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{\alpha-1} e^{-\beta y}, \quad y > 0, \alpha > 0, \beta > 0 \]

Para uma amostra aleatória \(Y_1, \dots, Y_n\) de uma distribuição Gama, a função de verossimilhança é:

\[ L(\alpha, \beta | \mathbf{y}) = \prod_{i=1}^n \frac{\beta^\alpha}{\Gamma(\alpha)} y_i^{\alpha-1} e^{-\beta y_i} \]

A função de log-verossimilhança, \(l(\alpha, \beta) = \ln(L(\alpha, \beta | \mathbf{y}))\), é:

\[ l(\alpha, \beta) = n \alpha \ln(\beta) - n \ln(\Gamma(\alpha)) + (\alpha-1) \sum_{i=1}^n \ln(y_i) - \beta \sum_{i=1}^n y_i \]

(a) Maximização Numérica da Função Log-Verossimilhança (\(l(\alpha, \beta)\))

Nesta abordagem, utilizaremos a função optim do R para maximizar a função de log-verossimilhança em relação a \(\alpha\) e \(\beta\) simultaneamente. A função optim realiza minimização por padrão, então maximizaremos o negativo da função de log-verossimilhança.

Primeiro, vamos gerar uma amostra de dados Gama para exemplificar.

set.seed(123) # Para reprodutibilidade
n_amostra <- 100
alpha_verdadeiro <- 2
beta_verdadeiro <- 0.5
dados_amostra <- rgamma(n_amostra, shape = alpha_verdadeiro, rate = beta_verdadeiro)

# Exibir os primeiros dados da amostra
head(dados_amostra)
## [1] 1.7841871 6.6236948 0.2887499 3.3250466 8.6717580 4.2352314

Agora, definimos a função de log-verossimilhança negativa.

log_verossimilhanca_negativa <- function(parametros, dados) {
  alpha <- parametros[1]
  beta <- parametros[2]
  n <- length(dados)
  
  # Condições para os parâmetros da Gama
  if (alpha <= 0 || beta <= 0) {
    return(NA) # Retorna NA para valores inválidos dos parâmetros
  }
  
  # Função de log-verossimilhança
  # l(alpha, beta) = n * alpha * ln(beta) - n * ln(Gamma(alpha)) + (alpha-1) * sum(ln(yi)) - beta * sum(yi)
  
  termo_alpha_ln_beta <- n * alpha * log(beta)
  termo_gamma_alpha <- n * lgamma(alpha) # lgamma é o log da função Gama
  termo_alpha_menos_1_ln_y <- (alpha - 1) * sum(log(dados))
  termo_beta_y <- beta * sum(dados)
  
  log_lik <- termo_alpha_ln_beta - termo_gamma_alpha + termo_alpha_menos_1_ln_y - termo_beta_y
  
  # Como optim minimiza, retornamos o negativo da log-verossimilhança
  return(-log_lik)
}

Vamos usar a função optim para encontrar os EMVs. Precisamos fornecer valores iniciais para \(\alpha\) e \(\beta\).

# Valores iniciais (podem ser baseados na média e variância da amostra)
# E[Y] = alpha/beta, Var[Y] = alpha/beta^2
# alpha_hat = media^2 / variancia
# beta_hat = media / variancia

media_amostra <- mean(dados_amostra)
variancia_amostra <- var(dados_amostra)

alpha_inicial <- media_amostra^2 / variancia_amostra
beta_inicial <- media_amostra / variancia_amostra

valores_iniciais <- c(alpha_inicial, beta_inicial)

# Realizar a otimização
resultado_optim <- optim(
  par = valores_iniciais,
  fn = log_verossimilhanca_negativa,
  dados = dados_amostra,
  method = "L-BFGS-B", # Método que permite restrições de caixa
  lower = c(0.001, 0.001) # alpha e beta devem ser positivos
)

# Extrair os EMVs
emv_alpha_a <- resultado_optim$par[1]
emv_beta_a <- resultado_optim$par[2]

cat("EMV de alpha (abordagem a):", emv_alpha_a, "\n")
## EMV de alpha (abordagem a): 2.193219
cat("EMV de beta (abordagem a):", emv_beta_a, "\n")
## EMV de beta (abordagem a): 0.6370087
cat("Valor da log-verossimilhança maximizada (negativo):", resultado_optim$value, "\n")
## Valor da log-verossimilhança maximizada (negativo): 209.2646

(b) Redução do Problema de Maximização para uma Função Univariada

Nesta abordagem, primeiro encontramos o EMV de \(\beta\) assumindo \(\alpha\) conhecido. Derivando a função de log-verossimilhança em relação a \(\beta\) e igualando a zero, obtemos:

\[ \frac{\partial l(\alpha, \beta)}{\partial \beta} = \frac{n \alpha}{\beta} - \sum_{i=1}^n y_i = 0 \]

\[ \hat{\beta}(\alpha) = \frac{n \alpha}{\sum_{i=1}^n y_i} = \frac{\alpha}{\bar{y}} \]

Substituindo \(\hat{\beta}(\alpha)\) de volta na função de log-verossimilhança, obtemos \(l(\alpha, \hat{\beta}(\alpha))\), que é uma função apenas de \(\alpha\):

\[ l(\alpha, \hat{\beta}(\alpha)) = n \alpha \ln\left(\frac{\alpha}{\bar{y}}\right) - n \ln(\Gamma(\alpha)) + (\alpha-1) \sum_{i=1}^n \ln(y_i) - \alpha n \]

\[ l(\alpha, \hat{\beta}(\alpha)) = n \alpha \ln(\alpha) - n \alpha \ln(\bar{y}) - n \ln(\Gamma(\alpha)) + (\alpha-1) \sum_{i=1}^n \ln(y_i) - n \alpha \]

(b.1) Abordagem do Grid de Valores para \(\alpha\)

Nesta sub-abordagem, avaliaremos a função \(l(\alpha, \hat{\beta}(\alpha))\) para uma série de valores de \(\alpha\) em um grid predefinido e identificaremos o \(\alpha\) que maximiza essa função. Faremos também um gráfico para visualizar o comportamento da função.

# Função de log-verossimilhança apenas em função de alpha
log_verossimilhanca_alpha <- function(alpha, dados) {
  n <- length(dados)
  media_y <- mean(dados)
  soma_ln_y <- sum(log(dados))
  
  if (alpha <= 0) {
    return(NA)
  }
  
  termo1 <- n * alpha * log(alpha)
  termo2 <- n * alpha * log(media_y)
  termo3 <- n * lgamma(alpha)
  termo4 <- (alpha - 1) * soma_ln_y
  termo5 <- n * alpha
  
  log_lik_alpha <- termo1 - termo2 - termo3 + termo4 - termo5
  
  return(log_lik_alpha)
}

Definindo o grid de valores para \(\alpha\) e avaliando a função.

alpha_grid <- seq(0.1, 600, by = 0.1)

# Calcular a log-verossimilhança para cada alpha no grid
valores_log_lik <- sapply(alpha_grid, log_verossimilhanca_alpha, dados = dados_amostra)

# Encontrar o alpha que maximiza a log-verossimilhança
indice_max <- which.max(valores_log_lik)
emv_alpha_b1 <- alpha_grid[indice_max]

# Calcular o beta correspondente
emv_beta_b1 <- emv_alpha_b1 / mean(dados_amostra)

cat("EMV de alpha (abordagem b.1):", emv_alpha_b1, "\n")
## EMV de alpha (abordagem b.1): 2.2
cat("EMV de beta (abordagem b.1):", emv_beta_b1, "\n")
## EMV de beta (abordagem b.1): 0.6389778
cat("Valor da log-verossimilhança maximizada:", valores_log_lik[indice_max], "\n")
## Valor da log-verossimilhança maximizada: -209.2649

Gerando o gráfico de \(l(\alpha, \hat{\beta}(\alpha))\) vs \(\alpha\).

plot(alpha_grid, valores_log_lik, type = "l", 
     xlab = expression(alpha), 
     ylab = expression(l(alpha, hat(beta)(alpha))), 
     main = expression("Função de Log-Verossimilhança em função de " * alpha))
points(emv_alpha_b1, valores_log_lik[indice_max], col = "red", pch = 19)
text(emv_alpha_b1, valores_log_lik[indice_max], 
     labels = paste0("EMV alpha = ", round(emv_alpha_b1, 2)), 
     pos = 4, col = "red")

(b.2) Abordagem Iterativa (Newton-Raphson ou similar)

Nesta sub-abordagem, utilizaremos um método iterativo para encontrar o \(\alpha\) que maximiza \(l(\alpha, \hat{\beta}(\alpha))\). A equação para \(\alpha\) é dada por:

\[ \ln(\alpha) - \psi(\alpha) = \ln(\bar{y}) - \frac{1}{n} \sum_{i=1}^n \ln(y_i) \]

Onde \(\psi(\alpha) = \frac{d}{d\alpha} \ln(\Gamma(\alpha))\) é a função digama. Não há solução analítica para \(\alpha\), então um método numérico é necessário. Usaremos a função optim novamente, mas agora para uma função univariada.

# A função a ser minimizada é o negativo da log-verossimilhança_alpha
resultado_optim_alpha <- optim(
  par = emv_alpha_b1, # Usamos o resultado do grid como valor inicial
  fn = function(alpha) -log_verossimilhanca_alpha(alpha, dados_amostra),
  method = "Brent", # Método para otimização univariada com restrições
  lower = 0.001, 
  upper = 600 # Definimos um limite superior razoável
)

emv_alpha_b2 <- resultado_optim_alpha$par[1]
emv_beta_b2 <- emv_alpha_b2 / mean(dados_amostra)

cat("EMV de alpha (abordagem b.2):", emv_alpha_b2, "\n")
## EMV de alpha (abordagem b.2): 2.193211
cat("EMV de beta (abordagem b.2):", emv_beta_b2, "\n")
## EMV de beta (abordagem b.2): 0.6370061
cat("Valor da log-verossimilhança maximizada (negativo):", resultado_optim_alpha$value, "\n")
## Valor da log-verossimilhança maximizada (negativo): 209.2646

Comparação dos Resultados

resultados <- data.frame(
  Abordagem = c("a (optim bivariado)", "b.1 (grid)", "b.2 (optim univariado)"),
  EMV_alpha = c(emv_alpha_a, emv_alpha_b1, emv_alpha_b2),
  EMV_beta = c(emv_beta_a, emv_beta_b1, emv_beta_b2)
)
print(resultados)
##                Abordagem EMV_alpha  EMV_beta
## 1    a (optim bivariado)  2.193219 0.6370087
## 2             b.1 (grid)  2.200000 0.6389778
## 3 b.2 (optim univariado)  2.193211 0.6370061

As três abordagens produzem resultados muito próximos para os Estimadores de Máxima Verossimilhança de \(\alpha\) e \(\beta\), o que demonstra a robustez dos métodos numéricos. O valor verdadeiro de \(\alpha\) era 2 e \(\beta\) era 0.5, e os EMVs obtidos estão razoavelmente próximos desses valores, como esperado para uma amostra de tamanho 100.