Subsections

23 Funções de verossimilhança

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.

23.1 Exemplo 1: Distribuição normal com variância conhecida

Seja o vetor uma amostra aleatória de uma distribuição normal de média $\mu$ 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:

Temos que $X_1, \ldots, X_n$ onde, neste exemplo , é uma a.a. de $X \sim N(\mu, 4)$ ,
a densidade para cada observação é dada por $f(x_i) = \frac{1}{2\sqrt{2\pi}} \exp\{-\frac{1}{8}(x_i - \mu)^2\}$ ,
a verossimilhança é dada por $L(\mu) = \prod_1^{10} f(x_i)$ ,
e a log-verossimilhança é dada por

$\displaystyle l(\mu)$ $\textstyle =$ $\displaystyle \sum_1^{10} \log(f(x_i))$

$\textstyle =$ $\displaystyle -5\log(8\pi) - \frac{1}{8} (\sum_1^{10}x_i^2 - 2\mu \sum_1^{10}x_i + 10\mu^2),$ (3)
que é uma função de $\mu$ e portanto devemos fazer um gráfico de $l(\mu)$ versus $\mu$ tomando vários valores de $\mu$ e calculando os valores de $l(\mu)$ .

Vamos ver agora uma primeira possível forma de fazer a função de verossimilhança no R.

Primeiro entramos com os dados que armazenamos no vetor x
```
> x <- c(12, 15, 9, 10, 17, 12, 11, 18, 15, 13)
```
e calculamos as quantidades $\sum_1^{10}x_i^2$ e $\sum_1^{10}x_i$
```
> sx2 <- sum(x^2)
> sx <- sum(x)
```
agora tomamos uma sequência de valores para $\mu$ . Sabemos que o estimador de máxima verossimilhança neste caso é $\hat{\mu} = 13.2$ (este valor pode ser obtido com o comando mean(x)) e portanto vamos definir tomar valores ao redor deste ponto.
```
> mu.vals <- seq(11, 15, l=100)
```

e a seguir calculamos os valores de $l(\mu)$ de acordo com a equação acima

> lmu <- -5 * log(8*pi) - (sx2 - 2*mu.vals*sx + 10*(mu.vals^2))/8

e finalmente fazemos o gráfico visto na Figura

> plot(mu.vals, lmu, type='l', xlab=expression(mu), ylab=expression(l(mu)))

**Figura** : Função de verossimilhança para o parâmetro $\mu$ da distribuição normal com variância $\sigma^2=4$ com os dados do Exemplo 1.
$\begin{figure}\centerline{\includegraphics[width=0.5\textwidth]{figuras/vero01.ps}}\end{figure}$

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.

Primeiro vamos criar uma função que calcula o valor da log-verossimilhança para um certo valor do parâmetro e para um certo conjunto de dados,
```
> logvero <- function(mu, dados){
     sum(dnorm(dados, mean = mu, sd = 2, log = TRUE))
  }
```
a seguir criamos uma sequência adequada de valores de $\mu$ e calculamos $l(\mu)$ para cada um dos valores
```
> mu.vals <- seq(11, 15, l=100)
> mu.vals
> lmu <- sapply(mu.vals, logvero, dados = x)
> lmu
```
Note na sintaxe acima que a função sapply aplica a função logvero anteriormente definida em cada elemento do vetor mu.vals.

Finalmente fazemos o gráfico.

> 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 $\sigma^2$ =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 $\mu$ 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)))

23.2 Exemplo 2: Distribuição Poisson

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 $\lambda$ . 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)

**Figura** : Função de verossimilhança e deviance para o parâmetro $\lambda$ da distribuição Poisson.
$\begin{figure}\centerline{\includegraphics[width=0.9\textwidth]{figuras/veropois.ps}}\end{figure}$

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 $\lambda = \bar{x}$ 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)

**Figura** : Função de verossimilhança e deviance para o parâmetro $\lambda$ da distribuição Poisson.
$\begin{figure}\centerline{\includegraphics{figuras/devpois.ps}}\end{figure}$

23.3 Exercícios

Seja a amostra abaixo obtida de uma distribuição Poisson de parâmetro $\lambda$ .
5 4 6 2 2 4 5 3 3 0 1 7 6 5 3 6 5 3 7 2
Obtenha o gráfico da função de log-verossimilhança.
Seja a amostra abaixo obtida de uma distribuição Binomial de parâmetro e com .
7 5 8 6 9 6 9 7 7 7 8 8 9 9 9
Obtenha o gráfico da função de log-verossimilhança.
Seja a amostra abaixo obtida de uma distribuição $\chi^2$ de parâmetro $\nu$ .
8.9 10.1 12.1 6.4 12.4 16.9 10.5 9.9 10.8 11.4
Obtenha o gráfico da função de log-verossimilhança.

Paulo Justiniano Ribeiro Jr

$\displaystyle l(\mu)$	$\textstyle =$	$\displaystyle \sum_1^{10} \log(f(x_i))$
	$\textstyle =$	$\displaystyle -5\log(8\pi) - \frac{1}{8} (\sum_1^{10}x_i^2 - 2\mu \sum_1^{10}x_i + 10\mu^2),$	(3)