|
Métodos Computacionais para Inferência Estatística |
Estimação por Máxima Verossimilhança: Distribuição gaussiana
Prof. Wagner Hugo Bonat LEG/UFPR
2019
1 Introdução
Suponha que \(Y_1, Y_2, \ldots, Y_n\) são variáveis aleatórias independentes e identicamente distribuídas com distribuição gaussiana de média \(\mu\) e variância \(\sigma^2\). Denote isto por \(Y_i \sim N(\mu, \sigma^2)\). Note que neste caso o vetor de parâmetros é \(\boldsymbol{\theta} = (\mu, \sigma)^{\top}\), onde \(\mu \in \Re\) e \(\sigma \in \Re^+\) são os respectivos espaços paramétricos. O objetivo é estimar \(\mu\) e \(\sigma\) além de encontrar intervalos ou regiões de confiança. Note que agora temos dois parâmetros e a função de log-verossimilhança é uma superfície.
Os princípios vistos no caso uniparamétrico (Tutorial I) se mantém, mas a construção de gráficos e a obtenção das estimativas são mais trabalhosas. Vejamos alguns fatos relevantes deste exemplo.
2 Estimação
Como sempre, começamos escrevendo a função de verossimilhança, \[ L(\mu, \sigma) = (2\pi)^{-n/2} \sigma^{-n} \exp\{ - \frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i - \mu)^2\}. \] A log-verossimilhança é dada por, \[ l(\mu,\sigma) = -\frac{n}{2} \log 2 \pi - n \log \sigma - \frac{1}{2 \sigma^2} \sum_{i=1}^n ( y_i - \mu)^2. \] A função escore toma a forma de um sistema de equações, \[\begin{align*} U(\mu) &= \frac{\partial l(\mu, \sigma)}{\partial \mu} = \frac{\sum_{i=1}^n y_i }{\sigma^2} - \frac{n \mu}{\sigma^2} \\ U(\sigma) &= -\frac{n}{\sigma} + \frac{1}{\sigma^3} \sum_{i=1}^n (y_i - \mu)^2 . \end{align*}\] Neste caso podemos facilmente resolver este sistema chegando as estimativas de máxima verossimilhança, \[ \hat{\mu} = \frac{\sum_{i=1}^n y_i}{n} \quad \text{e} \quad \hat{\sigma}^2 = \frac{\sum_{i=1}^n (y_i - \mu)^2}{n}. \] A matriz de informação observada fica da seguinte forma, \[ I_O (\mu, \sigma) = \left[\begin{array}{cc} -\frac{\partial^2 l(\mu, \sigma)}{\partial \mu^2} & - \frac{\partial^2 l(\mu, \sigma)}{\partial \mu \partial \sigma} \\ - \frac{\partial^2 l(\mu, \sigma)}{\partial \mu \partial \sigma} & -\frac{\partial^2 l(\mu, \sigma)}{\partial \sigma^2} \end{array}\right] . \] Tomando a esperança, temos \[\begin{align*} \frac{\partial^2 l(\mu, \sigma)}{\partial \mu^2} &= \frac{\partial U(\mu)}{\partial \mu} = - \frac{n}{\sigma^2} \\ \frac{\partial^2 l(\mu, \sigma)}{\partial \sigma^2} &= \frac{\partial U(\sigma)}{\partial \sigma} = - \frac{2n}{\sigma^2} \\ \frac{\partial^2 l(\mu, \sigma)}{\partial \mu \partial \sigma} &= \frac{\partial U(\sigma)}{\partial \sigma} = E\left(-\frac{2}{\sigma^3} \sum_{i=1}^n (y_i - \overline{y})\right) = 0. \end{align*}\]Logo, \[ I_E(\hat{\mu}, \hat{\sigma}) = \left[\begin{array}{cc} \frac{n}{\hat{\sigma}^2} & 0 \\ 0 & \frac{2n}{ \hat{\sigma}^2} \end{array}\right] . \]
Neste caso a matriz de informação observada coincide com a matriz de informação esperada. Além disso, note a importante propriedade de ortogonalidade entre os parâmetros, indicada pelos termos fora da diagonal da matriz de informação serem zero. A derivada cruzada entre dois parâmetros ser zero, é condição suficiente para que estes parâmetros sejam ortogonais. A ortogonalidade é uma propriedade muito conveniente e simplifica as inferências, uma vez que podemos fazer inferência para um parâmetro sem nos preocupar com os valores do outro.
3 Intervalos de confiança
Para construção dos intervalos de confiança, a maneira mais direta é usar a distribuição assintótica dos estimadores de máxima verossimilhança, neste caso temos que a distribuição assintótica de \(\hat{\underline{\theta}} = (\hat{\mu}, \hat{\sigma})^\top\) é \[ \begin{bmatrix} \hat{\mu} \\ \hat{\sigma} \end{bmatrix} \sim NM_2\left (\begin{bmatrix} \mu \\ \sigma \end{bmatrix} , \begin{bmatrix} \hat{\sigma}^2/n & 0 \\ 0 & \hat{\sigma}^2/2n \end{bmatrix} \right ) \] Intervalos de confiança de Wald podem ser obtidos por: \[ \hat{\mu} \pm z_{\alpha/2} \sqrt{\hat{\sigma}^2/n} \] e para \(\sigma\) temos \[ \hat{\sigma} \pm z_{\alpha/2} \sqrt{\hat{\sigma}^2/2n}. \]
A função de verossimilhança deste exemplo é simétrica e quadrática na direção de \(\mu\), e portanto a aproximação quadrática coincide com a forma exata na direção deste parâmetro.
Porém, a verossimilhança é assimétrica na direção de \(\sigma\). Destacamos ainda que a assimetria é maior \(\sigma^2\), um pouco menos acentuada em \(\sigma\) e ainda menos acentuada em uma transformação não linear como \(\psi = \log(\sigma)\). Na prática, se intervalos baseados na aproximação quadrática serão utilizados (por vezes a partir de hessianos numéricos), o mais recomendado então é reparametrizar a verossimilhança como função de \(\psi\) para uma obter uma forma mais próxima à simetria. Pode-se então obter intervalos assintóticos para \(\psi\) e depois transformá-los para escala original do parâmetro, por transformação direta se verossimilhança em \(\psi\) for muito próxima à simetria ou, caso contrário, pelo método delta.
Outra opção é obter uma região de confiança baseada na deviance, \[\begin{align*} D(\mu, \sigma) &= 2 [ l(\hat{\mu}, \hat{\sigma}) - l(\mu,\sigma)] \\ &= 2[ n \log \left( \frac{\sigma}{\hat{\sigma}} \right) + \frac{1}{2\sigma^2}\sum_{i=1}^n(y_i - \mu)^2 - \frac{1}{2\hat{\sigma}^2} \sum_{i=1}^n (y_i - \hat{\mu})]. \end{align*}\]A deviance aproximada tem a seguinte forma \[ D(\mu, \sigma) \approx ( \underline{\theta} - \underline{\hat{\theta}})^\top I_o(\underline{\hat{\theta}}) ( \underline{\theta} - \underline{\hat{\theta}}). \]
Note que neste caso a superfície de log-verossimilhança em duas dimensões esta sendo aproximada por uma elipse. É bastante intuitivo pensar que aproximar uma função em duas dimensões é mais difícil que em uma. Sabemos também que esta aproximação tenderá a ser melhor quanto maior for o tamanho da amostra. Para exemplificar esta ideia a Figura abaixo apresenta o gráfico da função deviance bidimensional em \((\mu, \sigma)\) para o caso do modelo gaussiano para quatro tamanhos de amostra, n=10, 50, 100 e 1000.
set.seed(123)
y10 <- rnorm(10,10,1.5)
y50 <- rep(y10,5)
y100 <- rep(y10,10)
y1000 <- rep(y10,100)
ll.gauss <- function(theta, y, logsigma=FALSE){
if(logsigma) theta[2] <- exp(theta[2])
ll <- sum(dnorm(y,theta[1], theta[2],log=TRUE))
return(ll)
}
n.m.sd10 <- c(length(y10), mean(y10), sd(y1000))
n.m.sd50 <- c(length(y50), mean(y50), sd(y1000))
n.m.sd100 <- c(length(y100), mean(y100), sd(y1000))
n.m.sd1000 <- c(length(y1000), mean(y1000), sd(y1000))
devi.approx <- function(theta, amostra, logsigma=FALSE){
if(logsigma) theta[2] <- exp(theta[2])
n <- amostra[1]
theta.est <- c(amostra[2],amostra[3]*(n-1)/n)
Io <- diag(c(n/theta.est[2]^2,2*n/theta.est[2]^2))
dev.app <- crossprod(theta - theta.est, Io) %*% (theta - theta.est)
return(dev.app)
}
Nseq <- 401
mu <- seq(8.2,11.8,l= Nseq)
sigma <- seq(0.4,3.5,l= Nseq)
grid.theta <- expand.grid(mu,sigma)
vero10 <- apply(as.matrix(grid.theta),1,ll.gauss,y=y10)
vero50 <- apply(as.matrix(grid.theta),1,ll.gauss,y=y50)
vero100 <- apply(as.matrix(grid.theta),1,ll.gauss,y=y100)
vero1000 <- apply(as.matrix(grid.theta),1,ll.gauss,y=y1000)
de.ap10 <- apply(as.matrix(grid.theta),1,devi.approx,amostra=n.m.sd10)
de.ap50 <- apply(as.matrix(grid.theta),1,devi.approx,amostra=n.m.sd50)
de.ap100 <- apply(as.matrix(grid.theta),1,devi.approx,amostra=n.m.sd100)
de.ap1000 <- apply(as.matrix(grid.theta),1,devi.approx,amostra=n.m.sd1000)
devi10 <- 2*(ll.gauss(theta=c(mean(y10),sd(y10)*0.9),y=y10) - vero10)
devi50 <- 2*(ll.gauss(theta=c(mean(y50),sd(y50)*0.98),y=y50) - vero50)
devi100 <- 2*(ll.gauss(theta=c(mean(y100),sd(y100)*0.99),y=y100) - vero100)
devi1000 <- 2*(ll.gauss(theta=c(mean(y1000),sd(y1000)*0.999),y=y1000) - vero1000)
par(mfrow=c(2,2),mar=c(2.8, 2.8, 1.7, 1),mgp = c(1.8, 0.7, 0))
LEVELS <- c(0.99,0.95,0.9,0.7,0.5,0.3,0.1,0.05)
contour(mu,sigma,matrix(devi10,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,
xlab=expression(mu),ylab=expression(sigma),
main="n=10", lwd=1.7)
par(new=TRUE)
contour(mu,sigma,matrix(de.ap10,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,
xlab=expression(mu),ylab=expression(sigma),axes=FALSE,
lty=1,col="red",lwd=1)
#
contour(mu,sigma,matrix(devi50,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,xlab=expression(mu),ylab=expression(sigma),
main="n=50",ylim=c(1,2),xlim=c(9.5,10.7), lwd=1.7)
par(new=TRUE)
contour(mu,sigma,matrix(de.ap50,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,
xlab=expression(mu),ylab=expression(sigma),axes=FALSE,
lty=1,col="red",lwd=1,ylim=c(1,2),xlim=c(9.5,10.7))
#
contour(mu,sigma,matrix(devi100,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,xlab=expression(mu),ylab=expression(sigma),
main="n=100",ylim=c(1.05,1.7),xlim=c(9.65,10.55), lwd=1.7)
par(new=TRUE)
contour(mu,sigma,matrix(de.ap100,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,xlab=expression(mu),ylab=expression(sigma),axes=FALSE,
lty=1,col="red",lwd=1,ylim=c(1.05,1.7),xlim=c(9.65,10.55))
#
contour(mu,sigma,matrix(devi1000,Nseq,Nseq),levels=qchisq(LEVELS,df=2),
labels=LEVELS,xlab=expression(mu),ylab=expression(sigma),
main="n=1000",xlim=c(9.98,10.25),ylim=c(1.25,1.48), lwd=1.7)
par(new=TRUE)
contour(mu,sigma,matrix(de.ap1000,Nseq,Nseq),levels=qchisq(LEVELS,df=2),
labels=LEVELS,xlab=expression(mu),ylab=expression(sigma),axes=FALSE,
lty=1,col="red",lwd=1,xlim=c(9.98,10.25),ylim=c(1.25,1.48))
Na Figura acima vemos que com \(n=10\) a verossimilhança exibe forte assimetria na direção do parâmetro \(\sigma\) e a aproximação quadrática é claramente insatisfatória. Com o aumento do tamanho da amostra a aproximação quadrática vai ficando cada vez mais próxima da deviance exata, mais uma vez ilustrando o comportamento assintótico da verossimilhança. É importante notar também que em modelos com dois ou mais parâmetros a aproximação pode melhorar mais rapidamente para um do que outro parâmetro. No exemplo a aproximação é exata para \(\mu\) para qualquer tamanho de amostra. Já para \(\sigma\) é necessário um tamanho de amostra relativamente grande para que a deviance em sua direção tome um comportamento próximo do simétrico. A função se aproxima da simetria mais rapidamente se parametrizada com \(\log(\sigma)\).
Em outros modelos, como no caso dos modelos lineares generalizados (MLG) a intuição é a mesma ainda que a aproximação para parâmetros de média deixe de ser exata. De forma geral, para parâmetros de média a aproximação quadrática tende a apresentar bons resultados mesmo com amostras reduzidas. O mesmo não pode ser dito para parâmetros de dispersão ou mesmo de correlação. Em outras palavras, estimar a média é mais simples que estimar a variabilidade, que por sua vez é mais simples do que estimar correlações.
4 Verossimilhança perfilhada
As regiões de confiança são as curvas de nível na superfície de deviance. Note que, apenas com a deviance não temos intervalos marginais como os obtidos pela aproximação quadrática. Uma possível solução é projetar a superfície na direção do parâmetro de interesse na maior amplitude, que é obtida quando fixamos o outro parâmetro na sua estimativa de máxima verossimilhança. Porém esta prática só produz bons resultados quando os parâmetros são ao menos aproximadamente ortogonais. Uma solução mais genérica, ainda que computacionalmente mais trabalhosa é o obtenção das verossimilhanças perfilhadas. No caso particular da distribuição gaussiana, que tem a propriedade de ortogonalidade entre \(\mu\) e \(\sigma\), a verossimilhança condicionada na estimativa de máxima verossimilhança coincide com a verossimilhança perfilhada. Para ilustrar este fato considere a obtenção da verossimilhança perfilhada para \(\mu\) e \(\sigma\) pelas funções a seguir:
## Perfil para mu
pl.mu <- function(sigma, mu, dados){
pll <- sum(dnorm(dados, mean=mu, sd=sigma, log=TRUE))
return(pll)}
## Perfil para sigma
pl.sigma <- function(mu, sigma, dados){
pll <- sum(dnorm(dados, mean=mu, sd=sigma, log=TRUE))
return(pll)}Vamos criar uma malha de valores nos quais a função será avaliada para a construção dos gráficos. Também vamos obter a log-verossimilhança condicionada na estimativa de máxima verossimilhança, que consiste em avaliar apenas a função para um dos parâmetros com o outro fixado em sua estimativa.
set.seed(123)
y10 <- rnorm(10,10,1.5)
grid.mu <- seq(9, 11.3, length=200)
grid.sigma <- seq(0.65, 2.7, length=200)
## Condicional para mu:
mu.cond <- sapply(grid.mu, pl.sigma, sigma=sqrt(var(y10)*9/10), dados=y10)
## Condicional para sigma:
sigma.cond <- sapply(grid.sigma, pl.mu, mu=mean(y10), dados=y10)Para obter o perfil de verossimilhança, por exemplo para \(\sigma\) precisamos de uma malha de valores de \(\sigma\) e para cada valor nesta malha encontrar o valor digamos \(\hat{\mu}_{\sigma}\) que maximiza a verossimilhança perfilhada. Para este exemplo existem formas fechadas para os estimadores, portanto basta aplicar a expressão do estimador de um parâmetro para cada valor na malha de valores do parâmetro sendo perfilhado. Entretanto para ilustração utilizamos uma forma mais geral, adequada para casos onde não há expressões fechadas, na qual maximizamos a função utilizando procedimentos numéricos, o que implica em uma otimização numérica de um parâmetro para cada valor na grade do parâmetro que está sendo perfilhado. Para a maximização usamos a função optimize() própria para maximização em apenas uma dimensão como é o caso neste exemplo. O código abaixo ilustra o procedimento para o conjunto de 10 dados.
O gráfico da esquerda da Figura abaixo mostra a superfície de deviance com as linhas tracejadas indicando os cortes para obtenção das deviances perfilhadas dos gráficos do centro e a direita. Nestes gráficos são também mostradas as funções deviance condicionadas no MLE (linha sólida).
mu.perf <- matrix(0, nrow=length(mu), ncol=2)
for(i in 1:length(mu)){
mu.perf[i,] <- unlist(optimize(pl.mu,c(0,200),
mu=mu[i],dados=y10,maximum=TRUE))}
sigma.perf <- matrix(0, nrow=length(sigma), ncol=2)
for(i in 1:length(sigma)){
sigma.perf[i,] <- unlist(optimize(pl.sigma,c(0,1000),
sigma=sigma[i],dados=y10,maximum=TRUE))}
LEVELS <- c(0.99,0.95,0.9,0.7,0.5,0.3,0.1,0.05)
par(mfrow=c(1,3), mar=c(3,3,1.5,1.5), mgp=c(1.7,0.8, 0))
contour(mu,sigma,matrix(devi10,Nseq,Nseq),
levels=qchisq(LEVELS,df=2),
labels=LEVELS,
xlab=expression(mu),ylab=expression(sigma))
lines(mu, mu.perf[,1], lty=2)
lines(sigma.perf[,1], sigma, lty=2)
plot(grid.mu,-2*(mu.cond - max(mu.cond)),type="l",ylab=expression(D(mu[sigma])),xlab=expression(mu),ylim=c(0,5), lwd=1.5,col="red")
lines(mu,-2*(mu.perf[,2] - max(mu.perf[,2])),lty=1,lwd=1)
abline(h=qchisq(0.95,df=1))
legend("top", c("Condicional", "Perfilhada"), lwd=c(1.5,1),col=c("black","red"),lty=c(1,1),cex=0.85)
plot(grid.sigma,-2*(sigma.cond - max(sigma.cond)),type="l",ylab=expression(D(sigma[mu])),xlab=expression(sigma),ylim=c(0,5), lwd=1.5,col="red")
lines(sigma,-2*(sigma.perf[,2] - max(sigma.perf[,2])),lty=1,lwd=1)
abline(h=qchisq(0.95,df=1))
legend("top", c("Condicional", "Perfilhada"), lwd=c(1.5,1),col=c("red","black"),lty=c(1,1),cex=0.85)
A Figura acima ilustra que a deviance perfilhada e a deviance condicional coincidem para o parâmetro \(\sigma\) porém não para o parâmetro de média \(\mu\). Isto reflete o fato de que ao perfilhar \(\sigma\) o valor maximizado \(\hat{\mu}_{\sigma} = \hat{\mu}\) não depende de \(\sigma\). Já no perfil de \(\mu\) cada um dos valores maximizados \(\hat{\sigma}_{\mu}\) dependem dos valores de \(\mu\). Para obter intervalos de confiança basta definir o corte nestas funções seja por valor relativo da verossimilhança ou usando o quantil da distribuição \(\chi^2\) para o nível de confiança desejado e encontrar as raízes da equação, assim como nos exemplos uniparamétricos. A verossimilhança perfilhada permite tratar um problema multiparamétrico como um problema uniparamétrico levando em consideração a forma da verossimilhança na direção de todos os parâmetros do modelo. Porém, esta abordagem pode ser extremamente cara computacionalmente, uma vez que para cada avaliação da verossimilhança perfilhada pode ser necessário uma maximização, que em geral vai requerer algum método numérico.
| Métodos Computacionais para Inferência Estatística | Prof. Wagner Hugo Bonat |
