A função de verossimilhança é central na inferência estatística. Nesta sessão vamos ver como traçar funções de verossimilhança utilizando o programa R.
Seja o vetor uma amostra aleatória de uma distribuição normal de média e variância conhecida e igual a . O objetivo é fazer um gráfico da função de log-verossimilhança.
Solução:
Vejamos primeiro os passos da solução analítica:
(3) |
Vamos ver agora uma primeira possível forma de fazer a função de verossimilhança no R.
> x <- c(12, 15, 9, 10, 17, 12, 11, 18, 15, 13)
> sx2 <- sum(x^2) > sx <- sum(x)
mean(x)
) e portanto vamos definir tomar valores ao redor deste ponto.
> mu.vals <- seq(11, 15, l=100)
> lmu <- -5 * log(8*pi) - (sx2 - 2*mu.vals*sx + 10*(mu.vals^2))/8
> plot(mu.vals, lmu, type='l', xlab=expression(mu), ylab=expression(l(mu)))
Entretanto podemos obter a função de verossimilhança no R de outras forma mais geral e menos trabalhosas.
Sabemos que a função dnorm
calcula a densidade da distribuição normal e
podemos usar este fato para evitar a digitação da expressão acima.
> logvero <- function(mu, dados){ sum(dnorm(dados, mean = mu, sd = 2, log = TRUE)) }
> mu.vals <- seq(11, 15, l=100) > mu.vals > lmu <- sapply(mu.vals, logvero, dados = x) > lmuNote na sintaxe acima que a função
sapply
aplica a função logvero
anteriormente definida em cada elemento do vetor mu.vals
.
> plot(mu.vals, lmu, type='l', xlab=expression(mu), ylab=expression(l(mu)))
Para encerrar este exemplo vamos apresentar uma solução ainda mais genérica que consiste em criar uma função que vamos chamar de vero.norm.v4
para cálculo da verossimilhança de distribuições normais com =4.
Esta função engloba os comandos acima e pode ser utilizada para obter o gráfico da log-verossimilhança para o parâmetro para qualquer amostra obtida desta distribuição.
> vero.normal.v4 <- function(mu, dados){ logvero <- function(mu, dados) sum(dnorm(dados, mean = mu, sd = 2, log = TRUE)) sapply(mu, logvero, dados = dados) } > curve(vero.normal.v4(x, dados = x), 11, 15, xlab=expression(mu), ylab=expression(l(mu)))
Considere agora a amostra armazenada no vetor y
:
> y <- c(5, 0, 3, 2, 1, 2, 1, 1, 2, 1)de uma distribuição de Poisson de parâmetro . A função de verossimilhança pode ser definida por:
lik.pois <- function(lambda, dados){ loglik <- function(l, dados){sum(dpois(dados, lambda = l, log = TRUE))} sapply(lambda, loglik, dados = dados) }
E podemos usar esta função para fazer o gráfico da função de verossimilhança como visto à esquerda da Figura
> lambda.vals <- seq(0, 10, l=101) > loglik <- sapply(lambda.vals, lik.pois, dados=y) > plot(lambda.vals, loglik, ty = "l") ## ou mudando o texto do eixos > plot(lambda.vals, loglik, type = "l", xlab=expression(lambda), > ylab=expression(l(lambda))) ## ou > curve(lik.pois(x, dados=y), 0, 10)
Alternativamente pode-se fazer um gráfico da função deviance, como nos comandos abaixo.
> dev.pois <- function(lambda, dados){ > lambda.est <- mean(dados) > lik.lambda.est <- lik.pois(lambda.est, dados = dados) > lik.lambda <- lik.pois(lambda, dados = dados) > return(-2 * (lik.lambda - lik.lambda.est)) > } > curve(dev.pois(x, dados=y), 0, 10) ## fazendo novamente em um intervalo menor > curve(dev.pois(x, dados=y), 0.5, 5)
O estimador de máxima verossimilhança é o valor que maximiza a função de verossimilhança que é o mesmo que minimiza a função deviance. Neste caso sabemos que o estimador tem expressão analítica fechada e portanto calculado com o comando.
> lambda.est [1] 1.8
Caso o estimador não tenha expressão fechada pode-se usar maximização
(ou minimização) numérica.
Para ilustrar isto vamos encontrar a estimativa do parâmetro da Poisson
e verificar que o valor obtido coincide com o valor dado pela expressão fechada do estimador.
Usamos o função optimise
para encontrar o ponto de mínimo da função deviance.
> optimise(dev.pois, int=c(0, 10), dados=y) $minimum [1] 1.800004 $objective [1] 1.075406e-10
A função optimise()
é adequada para minimizações envolvendo um único parâmetro.
Para dois ou mais parâmetros deve-se usar a função optim()
Finalmente os comandos abaixo são usados para obter graficamente o intervalo de confiança baseado na verossimilhança.
corte <- qchisq(0.95, df=1) abline(h=corte) ## obtendo os limites (aproximados) do IC l.vals <- seq(0.5,5,l=1001) dev.l <- dev.pois(l.vals, dados=y) dif <- abs(dev.l - corte) ind <- l.vals < lambda.est ic2.lambda <- c(l.vals[ind][which.min(dif[ind])], l.vals[!ind][which.min(dif[!ind])]) ic2.lambda ## adicionando ao gráfico curve(dev.pois(x, dados=y), 1, 3.5, xlab=expression(lambda), ylab=expression(l(lambda))) segments(ic2.lambda, 0, ic2.lambda, corte)