Vamos fazer aqui os exemplos vistos em sala de aula sobre intervalos de confiança baseado na deviance.
Seja
a.a. de uma distribuição normal de média
e variância 1.
Vimos que:
Vamos considerar que temos uma amostra onde e
.
Neste caso a função deviance é como mostrada na Figura 6
que é obtida com os comandos abaixo onde primeiro definimos uma função para calcular a deviance que depois é mostrada em um gráfico para valores entre 30 e 34.
Para obtermos um intervalo a 95% de confiança escolhemos o quantil correspondente na distribuição
e mostrado pela linha tracejada no gráfico. Os pontos onde esta linha cortam a função são, neste exemplo, determinados analiticamente pela expressão dada acima e indicados pelos setas verticais no gráfico.
> dev.norm.v1 <- function(theta, n, xbar){n * (xbar - theta)^2} > thetaN.vals <- seq(31, 33, l=101) > dev.vals <- dev.norm.v1(thetaN.vals, n=20, xbar=32) > plot(thetaN.vals, dev.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > corte <- qchisq(0.95, df = 1) > abline(h = corte, lty=3) > limites <- 32 + c(-1, 1) * sqrt(corte/20) > limites [1] 31.56174 32.43826 > segments(limites, rep(corte,2), limites, rep(0,2))
Vamos agora examinar o efeito do tamanho da amostra na função.
A Figura 7 mostra as funções para três tamanhos de amostra, e
que são obtidas com os comandos abaixo. A linha horizontal mostra o efeito nas amplitudes dos IC's.
> dev10.vals <- dev.norm.v1(thetaN.vals, n=10, xbar=32) > plot(thetaN.vals, dev10.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > dev20.vals <- dev.norm.v1(thetaN.vals, n=20, xbar=32) > lines(thetaN.vals, dev20.vals, lty=2) > dev50.vals <- dev.norm.v1(thetaN.vals, n=50, xbar=32) > lines(thetaN.vals, dev50.vals, lwd=2) > abline(h = qchisq(0.95, df = 1), lty=3) > legend(31, 2, c('n=10', 'n=20', 'n=50'), lty=c(1,2,1), lwd=c(1,1,2), cex=0.7)
Seja
a.a. de uma distribuição exponencial de parâmetro
com função de densidade
.
Vimos que:
A seguir vamos ilustrar a obtenção destes intervalos no R.
Vamos considerar que temos uma amostra onde e
para a qual
a função deviance é mostrada na Figura 8
e obtida
de forma análoga ao exemplo anterior.
> dev.exp <- function(theta, n, xbar){ + 2*n*(log((1/xbar)/theta) + xbar*(theta-(1/xbar))) + } > thetaE.vals <- seq(0.04,0.20, l=101) > dev.vals <- dev.exp(thetaE.vals, n=20, xbar=10) > plot(thetaE.vals, dev.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta)))
Neste exemplo, diferentemente do anterior, não determinamos a distribuição exata da deviance
e usamos a distribuição assintótica na qual se baseia a linha de corte tracejada
mostrada no gráfico para definir o IC do parâmetro ao nível de 95% de confiança.
Para encontrar os limites do IC precisamos dos valores no eixo dos parâmetros nos pontos onde a linha
de corte toca a função deviance o que corresponde a resolver a equação
onde
é quantil da distribuição da
com 1 grau de liberdade
correspondente ao nível de confiança desejado. Por exemplo, para 95% o valor de
é 3.84.
Como esta equação não tem solução analítica (diferentemente do exemplo anterior)
vamos examinar a seguir duas possíveis soluções para encontrar os limites do intervalo.
Iremos aqui considerar uma solução simples baseada no gráfico da função deviance para encontrar os limites do IC que consiste no seguinte:
Para fazermos o gráfico da deviance criamos uma sequência de valores do parâmetro . A cada um destes valores corresponde um valor de
. Vamos então localizar os valores de
para os quais
é o mais próximo possível do ponto de corte.
Isto é feito com o código abaixo e o resultado exibido na Figura 9.
> plot(thetaE.vals, dev.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > corte <- qchisq(0.95, df = 1) > abline(h = corte, lty=3) > dif <- abs(dev.vals - corte) > linf <- thetaE.vals[thetaE.vals<(1/10)][which.min(dif[thetaE.vals<(1/10)])] > lsup <- thetaE.vals[thetaE.vals>(1/10)][which.min(dif[thetaE.vals>(1/10)])] > limites.dev <- c(linf, lsup) > limites.dev [1] 0.0624 0.1504 > segments(limites.dev, rep(corte,2), limites.dev, rep(0,2))
Note que neste código procuramos primeiro e limite inferior entre os valores menores que a estimativa do parâmetro (1/10) e depois o limite superior entre os valores maiores que esta estimativa.
Embora este procedimento bastante simples e sujeito a imprecisão podemos torná-lo quão preciso quanto quisermos bastando para isto definir um vetor com menor espaçamento para os valores para o parâmetro, por exemplo poderiamos usar thetaE.vals <- seq(0.04,0.20,l=1001)
.
Nesta abordagem aproximamos a função deviance por uma função quadrática obtida pela expansão
por série de Taylor ao redor do estimador de máxima verossimilhança:
> devap.exp <- function(theta, n, xbar){n * (xbar *(theta - (1/xbar)))^2} > devap.vals <- devap.exp(thetaE.vals, n=20, xbar=10) > plot(thetaE.vals, devap.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > corte <- qchisq(0.95, df = 1) > abline(h = corte, lty=3) > limites.devap <- c((1/10)*(1 - sqrt(corte/20)), (1/10)*(1 + sqrt(corte/20))) > limites.devap [1] 0.05617387 0.14382613 > segments(limites.devap, rep(corte,2), limites.devap, rep(0,2))
![]() |
Examinando os limites dos intervalos encontrados anteriormente podemos ver que são diferentes. Vamos agora colocar os resultados pelos dois métodos em um mesmo gráfico (Figura 11) para comparar os resultados.
> plot(thetaE.vals, dev.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > lines(thetaE.vals, devap.vals, lty=2) > abline(h = corte, lty=3) > segments(limites.dev, rep(corte,2), limites.dev, rep(0,2)) > segments(limites.devap, rep(corte,2), limites.devap, rep(0,2), lty=2) > legend(0.07, 12, c('deviance', 'aproximacão quadrática'), lty=c(1,2), cex=0.7)
![]() |
Vamos agora examinar o efeito do tamanho da amostra na função deviance e sua aproximação quadrática.
A Figura 7 mostra as funções para três tamanhos de amostra,
e
que são obtidas com os comandos abaixo onde vemos que a aproximação fica cada vez melhor com o aumento do tamanho da amostra.
> thetaE.vals <- seq(0.04, 0.20, l=101) > dev10.vals <- dev.exp(thetaE.vals, n=10, xbar=10) > plot(thetaE.vals, dev10.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > devap10.vals <- devap.exp(thetaE.vals, n=10, xbar=10) > lines(thetaE.vals, devap10.vals, lty=2) > abline(h = qchisq(0.95, df = 1), lty=3) > > dev30.vals <- dev.exp(thetaE.vals, n=30, xbar=10) > plot(thetaE.vals, dev30.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > devap30.vals <- devap.exp(thetaE.vals, n=30, xbar=10) > lines(thetaE.vals, devap30.vals, lty=2) > abline(h = qchisq(0.95, df = 1), lty=3) > > dev100.vals <- dev.exp(thetaE.vals, n=100, xbar=10) > plot(thetaE.vals, dev100.vals, ty='l', + xlab=expression(theta), ylab=expression(D(theta))) > devap100.vals <- devap.exp(thetaE.vals, n=100, xbar=10) > lines(thetaE.vals, devap100.vals, lty=2) > abline(h = qchisq(0.95, df = 1), lty=3)
![]() |