|
Métodos Computacionais para Inferência Estatística |
Estimação por Máxima Verossimilhança: Distribuição exponencial
Prof. Wagner Hugo Bonat LEG/UFPR
2019
1 Introdução
Como vimos nos slides da aula há pelo menos três formas básicas para obtenção dos intervalos baseados na função de verossimilhança. A intuitivamente mais simples é baseada na verossimilhança relativa, isto é, o intervalo é definido pelo conjunto de valores do parâmetro que tenham uma verossimilhança não inferior a um certo percentual da verossimilhança máxima. O intervalo obtido desta forma não tem interpretação probabilística.
A segunda forma baseia-se na função deviance. Considerando a sua distribuição assintótica, define-se um ponto de corte nesta função e o intervalo é dado pelo conjunto dos valores do parâmetro que possuem deviance inferior a esse ponto de corte.
A terceira utiliza uma aproximação quadrática da função de log-verossimilhança. As últimas duas opções associam uma interpretação frequentista ao intervalo de confiança. Desta forma, os objetivos deste tutorial são: mostrar como obter os intervalos de confiança por cada uma das três abordagens, explorar as relações existentes e discutir algumas possíveis dificuldades computacionais que podem ocorrer quando trabalhamos diretamente com a função de verossimilhança relativa ou deviance.
2 Especificação do problema
Considere o caso onde amostras independentes de uma variável aleatória \(Y_i\) com \(i = 1, 2, \ldots,n\) cada uma tendo distribuição exponencial de parâmetro \(\theta\) está disponível. Vamos denotar esta situação por \(Y_i \sim {\rm Exp}(\theta)\). O objetivo é estimar o parâmetro \(\theta\) e um intervalo de confiança, digamos \(\hat{\theta}_L\) e \(\hat{\theta}_U\) que, sob o enfoque frequentista, tenha probabilidade \(1 - \alpha\) de conter \(\theta\). Note que a estrutura probabilística está em \(\hat{\theta}_L\) e \(\hat{\theta}_U\).
3 Estimação via máxima verossimilhança
O primeiro passo é sempre escrever a função de verossimilhança, \[L(\theta) = \prod_{i=1}^n \theta \exp\{-\theta y_i\} = \theta^n \exp\{-\theta \sum_{i=1}^n y_i\} = \theta^n \exp\{-n \overline{y} \theta \}.\]
Consequentemente, a função de log-verossimilhança é \[l(\theta) = n \log \theta - \theta \sum_{i=1}^n y_i = n \log \theta - \theta n \overline{y}.\]
Derivando a log-verossimilhança em \(\theta\) chegamos a função escore \[U(\theta) = n \theta^{-1} - n \overline{y}.\]
Para encontrar a estimativa de máxima verossimilhança basta igualar a função escore a \(0\) e isolar o \(\theta\) isto resulta \[\hat{\theta} = \frac{1}{\overline{y}}.\]
A segunda derivada da função de log-verossimilhança tem um papel fundamental na construção de intervalos de confiança por aproximação quadrática da função de log-verossimilhança, uma vez que ela mede a curvatura da função no entorno do ponto de máximo.
Em termos estatísticos, trabalhamos com a informação observada ou esperada, uma vez que esta quantidade descreve o comportamento assintótico dos estimadores de máxima verossimilhança, já que sua inversa fornece a variância do estimador. Neste caso a função de log-verossimilhança é côncava e polinomial de ordem infinita. As informações esperada e observada são iguais e dependem do valor do parâmetro e na prática é necessário usar a informação estimada pela amostra. \[I_O(\theta) = I_E(\theta) = n \theta^{-2} \;\; \mbox{ e } \;\; I_O(\hat{\theta}) = I_E(\hat{\theta}) = n \hat{\theta}^{-2}.\]
4 Intervalo deviance
Para encontrar o intervalo de confiança vamos começar pela função deviance, lembrando que de acordo com os resultados apresentados em aula a deviance segue uma distribuição \(\chi^2\) com \(p\) graus de liberdade, no exemplo \(p = 1\). Desta forma o valor de \(c^*\) é obtido definindo algum quantil \(q_{1-\alpha}\) da distribuição qui-quadrado com \(1\) grau de liberdade e \(1-\alpha\) de confiança. Como exemplo, fixando \(\alpha = 0,05\) toma-se o valor do quantil \(c^* = 3,84\). Com isto, temos os elementos necessários para construir o intervalo de confiança baseado na função que é dado pelo conjunto de valores que obedecem: \[\begin{eqnarray}\label{devexponencial} D(\theta) &=& 2 [ l(\hat{\theta}) - l(\theta)] \nonumber \\ &=& 2 [ n \log \hat{\theta} - \hat{\theta} n \overline{y} - ( n \log \theta - \theta n \overline{y}) ] \nonumber \\ &=& 2n [ \log \left(\hat{\theta}/\theta \right) + \overline{y} ( \theta - \hat{\theta}) ] \leq c^* \end{eqnarray}\]Os limites \((\hat{\theta}_L, \hat{\theta}_U)\) são encontrados resolvendo a equação acima, mas mesmo em uma situação simples como a deste exemplo a obtenção das raízes requer a solução de uma equação não-linear que não tem solução analítica. Desta forma, algum método numérico é necessário para encontrar os limites do intervalo.
5 Intervalo de Wald
Uma outra opção para encontrar o intervalo de confiança é fazer uma aproximação quadrática utilizando séries de Taylor para a função \(l(\theta)\) em torno do ponto \(\hat{\theta}\) e usar esta aproximação no lugar da função de log-verossimilhança original para obter o intervalo de confiança. Neste caso, os limites do intervalo são as raízes de um polinômio de segundo grau e os intervalos são da forma \[\hat{\theta} \pm z_{\frac{\alpha}{2}} \sqrt{ V(\hat{\theta})},\] o que corresponde a usar a distribuição assintótica do EMV.
Para a distribuição exponencial, temos que a \(V(\hat{\theta}) = I_O^{-1}(\hat{\theta}) = \hat{\theta}^2/n\), logo o intervalo de confiança é dado por \[\hat{\theta}_L = \hat{\theta} - z_{\frac{\alpha}{2}} \hat{\theta}/\sqrt{n} \quad \text{e} \quad \hat{\theta}_U = \hat{\theta} + z_{\frac{\alpha}{2}} \hat{\theta}/\sqrt{n} . \]
Aqui a construção do intervalo de confiança fica bastante facilitada. Conhecendo o valor de \(z_{\frac{\alpha}{2}}\) o intervalo se resume a uma conta simples, ou seja, não é necessário utilizar métodos numéricos para obtenção dos limites do intervalo. Como a distribuição amostral de \(\hat{\theta}\) é assintóticamente gaussiana podemos escolher o valor de \(z_{\frac{\alpha}{2}}\) como o quantil da distribuição normal, com \((1-\alpha)\) 100% de confiança, escolhendo \(\alpha = 0,05\), temos aproximadamente que \(z_{\frac{\alpha}{2}} = 1,96\). No caso de \(d=1\) o quantil da distribuição \(\chi^2_{(1)}\) é o quadrado de um score \(z\) da distribuição normal e o método também corresponde a obter o intervalo como na função deviance porém utilizando a sua aproximação quadrática dada por: \[ D(\theta) \approx n \left(\frac{\theta-\hat{\theta}}{\hat{\theta}} \right)^2. \] Desta forma as raízes que definem os limites dos intervalos equivalentes aos anteriores são: \[ \left(\hat{\theta}_L \approx \hat{\theta} (1-\sqrt{c^*/n}) \;\; , \;\; \hat{\theta}_U \approx \hat{\theta} (1+\sqrt{c^*/n})\right) . \]
6 Intervalo baseado em perda relativa
A última opção que vamos considerar é a construção do intervalo definindo um valor limite para a verossimilhança relativa. Um intervalo baseado em verossimilhança relativa é da forma \(LR(\theta) = \frac{L(\theta)}{L(\hat{\theta})} \ge r\), onde \(r\) representa um percentual do máximo da função de verossimilhança, para o qual admite-se ainda obter uma estimativa aceitável de \(\theta\). Por exemplo, definindo \(r=0.25\), o intervalo é definido pelo valores que apresentam uma verosimilhança que seja ao menos 25% da máxima para o conjunto de dados. Para escolher um valor de \(r\) comparável com o escolhido para a deviance, note o seguinte: \[\begin{eqnarray*} \frac{L(\theta)}{L(\hat{\theta})} &\ge& r \\ \log \left[ \frac{L(\theta)}{L(\hat{\theta})} \right] &\ge& \log r \\ -2[ l(\theta) - l(\hat{\theta})] &\ge& -2 \cdot \log r\\ 3,84 & \ge& -2 \cdot \log r \\ r &\approx& 0,146 \end{eqnarray*}\]A Tabela abaixo mostra a relação entre alguns valores de corte da verossimilhança relativa e os valores correspondentes em log-verossimilhança e em deviance. Na prática, os níveis de confiança são aproximados e válidos apenas assintóticamente.
| \(r\) | \(c\) | \(c^*\) | Nível de confiança |
|---|---|---|---|
| 50% | 0,693 | 1,386 | 0,761 |
| 26% | 1,661 | 3,321 | 0,899 |
| 15% | 1,897 | 3,794 | 0,942 |
| 3,6% | 3,324 | 6,648 | 0,990 |
Com isso, podemos observar que usar a verossimilhança relativa ou a deviance são formas equivalentes de construir o intervalo de confiança. A diferença está na interpretação que é associada aos intervalos. A vantagem em usar a deviance é que conhecemos a sua distribuição assintótica o que permite a construção de intervalos de confiança com a desejada estrutura probabilística. Da mesma forma que para a deviance a obtenção do intervalo pela verossimilhança relativa também requer resolver uma equação não-linear utilizando algum método numérico para encontrar as raízes da equação \(LR(\theta) \ge r\).
7 Implementação computacional
Para ilustrar a construção dos intervalos de confiança, vamos considerar uma amostra de tamanho \(n = 20\) gerada de uma distribuição exponencial com parâmetro \(\theta = 1\).
set.seed(123)
y <- rexp(20, rate=1)Vamos escrever funções em R para avaliar as funções de verossimilhança e log-verossimilhança. Baseado nestas funções podemos avaliar a verossimilhança relativa e a deviance. Por fim, criamos uma função para avaliar a deviance aproximada usando a aproximação em série de Taylor.
grid.theta <- seq(0.68,2,l=100)
## Verossimilhança
Lik <- function(theta,y){
Lt <- prod(dexp(y,rate=theta,log=FALSE))
return(Lt)
}
L.theta <- sapply(grid.theta, Lik, y=y)
## Verossimilhança relativa
LR.theta <- L.theta/Lik(theta=1/mean(y),y=y)
## Log-verossimilhança
ll <- function(theta,y){
llt <- sum(dexp(y,rate=theta,log=TRUE))
return(llt)
}
ll.theta <- sapply(grid.theta, ll, y=y)
## Deviance
dev.theta <- 2*(ll(1/mean(y),y=y)-ll.theta)
## Deviance aproximada
dev.aprox <- function(theta,theta.hat,y){
dev <- (theta - theta.hat)^2*(length(y)/(theta.hat^2))
return(dev)
}
devA.theta <- sapply(grid.theta, dev.aprox, theta.hat = 1/mean(y), y=y)Vamos plotar alguns gráficos para ver o comportamento das funções.
theta.est <- round(1/mean(y), dig=2)
par(mfrow = c(1,2), mar=c(2.6, 2.8, 1.2, 0.5), mgp = c(1.6, 0.6, 0))
plot(LR.theta ~ grid.theta,type="l",ylab=expression(LR(theta)),xlab=expression(theta))
abline(h=0.146)
abline(v = theta.est)
plot(dev.theta ~ grid.theta,type="l",ylab=expression(D(theta)),xlab=expression(theta))
lines(devA.theta ~ grid.theta,type="l",col="red",lty=1)
abline(h=3.84)
abline(v = theta.est)
legend("top",legend=c("Deviance"," Aproximada"),col=c("black","red"), lty=c(1,1), bty="n")
Como é possível ver na Figura acima a função deviance é razoavelmente aproximada por uma forma quadrática mas esta aproximação vai depender do valor do parâmetro e do tamanho da amostra. Para obter os limites inferior e superior do intervalo de confiança pela aproximação quadrática temos: \[ \hat{\theta}_L = \hat{\theta} - z_{\frac{\alpha}{2}} \sqrt{\hat{\theta}^2/n} \quad \text{e} \quad \hat{\theta}_U = \hat{\theta} + z_{\frac{\alpha}{2}} \sqrt{\hat{\theta}^2/n} \] Para a amostra simulada utilizada para gerar os gráficos temos os seguintes valores: \[ \hat{\theta}_L = 1.23(1 - 1.96 \sqrt{1/20}) \quad \text{e} \quad \hat{\theta}_U = 1.23(1 + 1.96 \sqrt{1/20}) , \] que resulta no seguinte intervalo de confiança: \[ \hat{\theta}_L = 0.6909; \quad \hat{\theta}_U = 1.7690.\]
Usando a função deviance precisamos resolver uma equação não-linear. Em R podemos facilmente resolver este problema usando as funções do pacote rootSolve. O algoritmo implementado na função uniroot() é uma variação do método de Brent que veremos na aula sobre métodos de otimização. De forma geral é um algoritmo de confinamento e por isso precisamos especificar os limites inferior e superior de busca. Lembrando que nos limites a função deve ter sinais trocados.
O primeiro passo é escrever a função deviance,
ICdevExp <- function(theta, theta.hat, y, nivel=0.95){
n <- length(y)
dv <- 2*n*( log( theta.hat/theta) + mean(y)*(theta- theta.hat))
return(dv - qchisq(nivel,df=1))
}Uma vez com a função escrita podemos encontrar suas raízes usando a função uniroot.all().
require(rootSolve)
# Loading required package: rootSolve
uniroot.all(ICdevExp, interval=c(0,10), theta.hat=1/mean(y),y=y)
# [1] 0.7684028 1.8547415A Figura acima mostra o intervalo aproximado pela forma quadrática com um deslocamento para a esquerda quando comparado com o intervalo baseado na função deviance. É importante ter bastante cuidado na interpretação destes intervalos.
De acordo com a interpretação frequentista de probabilidade, se realizarmos o mesmo experimento um grande número de vezes e em cada um destes experimentos calcularmos o respectivo intervalo de confiança esperamos que em \((1-\alpha)\) 100% dos intervalos construídos contenham o verdadeiro valor do parâmetro. Isto pode ser ilustrado com o seguinte estudo de simulação em que geramos 100 amostras cada uma de tamanho \(n=70\), com valor do parâmetro igual a 1. Ao final do experimento verificamos a taxa de cobertura, isto é, construímos o intervalo de confiança e ao final verificamos a proporção dos intervalos que contém o verdadeiro valor do parâmetro.
THETA <- 1
set.seed(12)
ic <- matrix(NA, ncol=2, nrow=100)
for(i in 1:100){
y <- rexp(70, rate=THETA)
est <- 1/mean(y)
ic[i,] <- uniroot.all(ICdevExp, int=c(0,5), theta.hat=est, y=y)
}
mean(apply(ic, 1, function(x) (x[1] < THETA & x[2] > THETA)))
# [1] 0.95No código acima simulamos a cada passo do laço for() uma nova realização da variável aleatória \(Y\), com esta realização calculamos a estimativa de máxima verossimilhança e o seu respectivo intervalo de confiança baseado na deviance e guardamos o resultado em um objeto chamado \(ic\). De acordo com a interpretação frequentista dos \(100\) intervalos de confiança que calculamos esperamos que \(95\) deles contenham o verdadeiro valor do parâmetro neste caso \(\theta = 1\).
O gráfico apresentado abaixo ilustra os resultados. Neste caso, conforme esperado, temos exatamente que \(5\) dos intervalos não contem o valor verdadeiro do parâmetro. É claro que por se tratar de um exemplo de simulação variações podem ocorrer.
plot(c(0.4,1.8)~c(1,100),type="n",ylab=expression(theta),xlab="Ensaio")
abline(h=1)
for(i in 1:100){
arrows(c(i),ic[i,1],c(i),ic[i,2],code=3,angle=90,length=0.03,
col=ifelse(ic[i,1] > 1 | ic[i,2] < 1, "darkred","lightgray"))
# , lty=ifelse(ic[i,1] > 1 | ic[i,2] < 1, 1, 2))
}
A taxa de cobertura de um intervalo é a proporção de vezes em que os intervalos contêm o verdadeiro valor do parâmetro sobre o total de ensaios simulados. Neste caso obtivemos exatamente \(95/100\), ou seja, a taxa de cobertura é de \(95\)% igual ao nível nominal especificado para este intervalo. A taxa de cobertura é uma forma muito utilizada para avaliar e comparar métodos de construção de intervalos de confiança.
| Métodos Computacionais para Inferência Estatística | Prof. Wagner Hugo Bonat |
