Nesta sessão vamos examinar um pouco mais a teoria de intervalos de confiança. São ilustrados os conceitos de:
Serão utilizados conceitos do método da quantidade pivotal, a propriedade de normalidade assintótica dos estimadores de máxima verossimilhança e a distribuição limite da função deviance.
Os dados abaixo são uma amostra aleatória da distribuição \(Bernoulli(p)\).
0 0 0 1 1 0 1 1 1 1 0 1 1 0 1 1 1 1 0 1 1 1 1 1 1
Desejamos obter:
o gráfico da função de verossimilhança para \(p\) com base nestes dados
o estimador de máxima verossimilhança de \(p\), a informação observada e a informação de Fisher
um intervalo de confiança de 95% para \(p\) baseado na normalidade assintótica de \(\hat {p}\)
compare o intervalo obtido em (b) com um intervalo de confiança de 95% obtido com base na distribuição limite da função deviance
a probabilidade de cobertura dos intervalos obtidos em (c) e (d). (O verdadeiro valor de \(p\) é 0.8)
Primeiramente vamos entrar com os dados na forma de um vetor.
> y <- c(0,0,0,1,1,0,1,1,1,1,0,1,1,0,1,1,1,1,0,1,1,1,1,1,1)
(a)
Vamos escrever uma função em R para obter a função de verossimilhança usando a função de
densidade da distribuição binomial com argumento log=TRUE para obter a log-verossimilhança.
> vero.binom <- function(p, dados){ + n <- length(dados) + x <- sum(dados) + return(dbinom(x, size = n, prob = p, log = TRUE)) + }
Esta função exige dados do tipo \(0\) ou \(1\) da distribuição Bernoulli. Entretanto às vezes temos dados binomiais do tipo \(n\) e \(x\) (número \(x\) de sucessos em \(n\) observações). Por exemplo, para os dados acima teríamos \(n=25\) e \(x=18\). Vamos então escrever a função acima de forma mais geral de forma a poder utilizar dados disponíveis tanto em um formato quanto em outro.
> vero.binom <- function(p, dados, n = length(dados), x = sum(dados)){ + return(dbinom(x, size = n, prob = p, log = TRUE)) + }
Para obter o gráfico da função de verossimilhança de um conjunto de dados cria-se uma sequência de valores para o parâmetro \(p\) e calcula-se o respectivo valor da (log)verossimilhança. O gráfico da função é obtido com os valores fixados dos parâmetros no eixo-\(x\) e o respectivos valores da função no eixo-\(y\) e unindo-se os pontos assim obtidos. No R isto pode ser feito com os comandos abaixo que produzem o gráfico mostrado na Figura 38. Evitamos os valores nos extremos do espaço paramétrico (\(p=0\) ou \(p=1\)) pois nestes casos a verossimilhança é zero e portanto a log-verossimilhança retornada por dbinom() é -Inf.
> p.vals <- seq(0.01,0.99,l=99) > logvero <- sapply(p.vals, vero.binom, dados=y) > plot(p.vals, logvero, type="l", xlab="p", ylab="l(p)")
Note que os três comandos acima podem ser substituídos por um único que produz o mesmo resultado:
> curve(vero.binom(x, dados=y), from = 0, to = 1)
(b)
Dos resultados para distribuição Bernoulli sabemos que o estimador de máxima verossimilhança é
dado por
\[\hat {p} = \frac {\sum _{i=1}^n y_i}{n}\]
e que a informação esperada coincide com a informação observada e sendo iguais a:
\[I(\hat {p}) = \frac {n}{\hat {p}(1-\hat {p})}.\]
Para indicar o estimador de MV o gráfico poderíamos usar arrows() e para obter os valores
numéricos para a amostra dada utilizamos os comandos a seguir.
> p.est <- mean(y) > io <- ie <- length(y)/(p.est * (1 - p.est)) > io
[1] 124.0079
> ie
[1] 124.0079
(c)
O intervalo de confiança baseado na normalidade assintótica do estimador de máxima
verossimilhança é dado por:
\[\left (\hat {p} - z_{\alpha /2} \sqrt {I(\hat {p})}\;\; , \;\;\hat {p} + z_{\alpha /2} \sqrt {I(\hat {p})}\right )\]
e para obter o intervalo no R usamos os comandos a seguir.
> ic1.p <- p.est + qnorm(c(0.025, 0.975)) * sqrt(1/ie) > ic1.p
[1] 0.5439957 0.8960043
(d)
Vamos agora obter e mostrar em um gráfico o intervalo baseado na função deviance. Lembrando que
a deviance é definida pela expressão
\[D(p) = 2\{l(\hat {p}) - l(p)\} , \]
definimos uma função dev.binom() para calcular a deviance. Com o comando curve() podemos
obter o gráfico de função deviance.
> dev.binom <- function(p, dados, n = length(dados), x =sum(dados)){ + p.est <- x/n + vero.p.est <- vero.binom(p.est, n = n, x = x) + dev <- 2 * (vero.p.est - vero.binom(p, n = n, x = x)) + dev + } > curve(dev.binom(x, dados=y), 0.35, 0.95, xlab="p", ylab="D(p)")
A função deviance \(D(p)\) tem distribuição assintótica \(\chi ^2_{n-1}\) e o intervalo de confiança é dado pelos pontos de intersecção entre a função deviance e o valor de quantil da distribuição \(\chi ^2\) para o nível de significância desejado como ilustrado na Figura 39. Nos comandos a seguir primeiro encontramos o ponto de corte para o nível de confiança de 95%. Depois traçamos a linha de corte com abline() e os pontos de corte que definem o intervalo são as raízes de uma função definida como a diferença entre a função deviance e o valor do ponto de corte.
> L.95 <- qchisq(0.95, df=1) > #dif <- abs(dev.p - L.95) > #inf <- ifelse(p.est==0, 0, p.vals[p.vals<p.est][which.min(dif[p.vals<p.est])]) > #sup <- ifelse(p.est==1, 1, p.vals[p.vals>p.est][which.min(dif[p.vals>p.est])]) > #c(inf,sup) > lim.fc <- function(x) dev.binom(x, dados=y) - L.95 > ICdev <- c(inf=uniroot(lim.fc, c(0, p.est))$root, + sup=uniroot(lim.fc, c(p.est, 1))$root) > ICdev
inf sup 0.5283461 0.8686755
Agora que já vimos as duas formas de obter o IC passo a passo vamos usar os comandos acima para criar uma função geral para encontrar IC para qualquer conjunto de dados e com opções para os dois métodos.
> ic.binom <- function(dados, n=length(dados), x=sum(dados), + nivel = 0.95, tipo=c("assintotico", "deviance")){ + tipo <- match.arg(tipo) + alfa <- 1 - nivel + p.est <- x/n + if(tipo == "assintotico"){ + se.p.est <- sqrt((p.est * (1 - p.est))/n) + ic <- p.est + qnorm(c(alfa/2, 1-(alfa/2))) * se.p.est + } + if(tipo == "deviance"){ + lim.fc <- function(y, ...) dev.binom(y, ...) - qchisq(nivel, df=1) + inf <- ifelse(identical(p.est, 0), 0, uniroot(lim.fc, c(0, p.est), n=n, x=x)$root) + sup <- ifelse(identical(p.est, 1), 1, uniroot(lim.fc, c(p.est, 1), n=n, x=x)$root) + ic <- c(inf, sup) + } + names(ic) <- c("lim.inf", "lim.sup") + ic + }
E agora vamos utilizar a função, primeiro com a aproximação assintótica e depois pela deviance. Note que os intervalos são diferentes!
> ic.binom(dados=y)
lim.inf lim.sup 0.5439957 0.8960043
> ic.binom(dados=y, tipo = "dev")
lim.inf lim.sup 0.5283461 0.8686755
(e)
O cálculo do intervalo de cobertura consiste em:
simular dados com o valor especificado do parâmetro;
obter o intervalo de confiança;
verificar se o valor está dentro do intervalo
repetir (1) a (3) e verificar a proporção de simulações onde o valor está no intervalo.
Espera-se que a proporção obtida seja o mais próximo possível do nível de confiança definido para o intervalo.
Para isto vamos escrever uma função implementando estes passos e que utiliza internamente ic.binom() definida acima.
> cobertura.binom <- function(n, p, nsim, ...){ + conta <- 0 + for(i in 1:nsim){ + ysim <- rbinom(1, size = n, prob = p) + ic <- ic.binom(n = n, x = ysim, ...) + if(p > ic[1] & p < ic[2]) conta <- conta+1 + } + return(conta/nsim) + }
E agora vamos utilizar esta função para cada um dos métodos de obtenção dos intervalos.
> set.seed(3214) > cobertura.binom(n=length(y), p=0.8, nsim=1000)
[1] 0.897
> set.seed(3214) > cobertura.binom(n=length(y), p=0.8, nsim=1000, tipo = "dev")
[1] 0.96
Note que a cobertura do método baseado na deviance é muito mais próxima do nível de 95%, o que pode ser explicado pelo tamanho da amostra. O IC assintótico tende a se aproximar do nível nominal de confiança na medida que aumenta o tamanho da amostra.
Refaça o item (e) do exemplo acima com \(n = 10\), \(n = 50\) e \(n=200\). Discuta os resultados.
Seja \(X_1, X_2, \cdots , X_n\) uma amostra aleatória da distribuição \(U(0,\theta )\). Encontre uma quantidade pivotal e:
construa um intervalo de confiança de 90% para \(\theta \)
construa um intervalo de confiança de 90% para \(\log \theta \)
gere uma amostra de tamanho \(n=10\) da distribuição \(U(0,\theta )\) com \(\theta =1\) e obtenha o intervalo de confiança de 90% para \(\theta \). Verifique se o intervalo cobre o verdadeiro valor de \(\theta \).
verifique se a probabilidade de cobertura do intervalo é consistente com o valor declarado de 90%. Para isto gere 1000 amostras de tamanho \(n=10\). Calcule intervalos de confiança de 90% para cada uma das amostras geradas e finalmente, obtenha a proporção dos intervalos que cobrem o verdadeiro valor de \(\theta \). Espera-se que este valor seja próximo do nível de confiança fixado de 90%.
repita o item (d) para amostras de tamanho \(n=100\). Houve alguma mudança na probabilidade de cobertura?
Note que se \(-\sum _i^n \log F(x_i;\theta ) \sim \Gamma (n,1)\) então \(-2 \sum _i^n \log F(x_i;\theta ) \sim \chi ^2_{2n}\).
Acredita-se que o número de trens atrasados para uma certa estação de trem por dia segue uma distribuição Poisson(\(\theta \)), além disso acredita-se que o número de trens atrasados em cada dia seja independente do valor de todos os outros dias. Em 10 dias sucessivos, o número de trens atrasados foi registrado em:
5 0 3 2 1 2 1 1 2 1
Obtenha:
o gráfico da função de verossimilhança para \(\theta \) com base nestes dados
o estimador de máxima verossimilhança de \(\theta \), a informação observada e a informação de Fisher
um intervalo de confiança de 95% para o número médio de trens atrasados por dia baseando-se na normalidade assintótica de \(\hat {\theta }\)
compare o intervalo obtido em (c) com um intervalo de confiança obtido com base na distribuição limite da função deviance
o estimador de máxima verossimilhança de \(\phi \), onde \(\phi \) é a probabilidade de que não haja trens atrasados num particular dia. Construa intervalos de confiança de 95% para \(\phi \) como nos itens (c) e (d).
Encontre intervalos de confiança de 95% para a média de uma distribuição Normal com variância 1 dada a amostra
9.5 10.8 9.3 10.7 10.9 10.5 10.7 9.0 11.0 8.4
10.9 9.8 11.4 10.6 9.2 9.7 8.3 10.8 9.8 9.0
baseando-se:
na distribuição assintótica de \(\hat {\mu }\)
na distribuição limite da função deviance
Acredita-se que a produção de trigo, \(X_i\), da área \(i\) é normalmente distribuída com média \(\theta z_i\), onde \(z_i\) é quantidade (conhecida) de fertilizante utilizado na área. Assumindo que as produções em diferentes áreas são independentes, e que a variância é conhecida e igual a 1, ou seja, \(X_i \sim N(\theta z_i, 1)\), para \(i=1,\cdots ,n\):
simule dados sob esta distribuição assumindo que \(\theta =1.5\), e \(z=(1,2,3,4,5)\). Visualize os dados simulados através de um gráfico de \((z \times x)\)
encontre o EMV de \(\theta \), \(\hat {\theta }\)
mostre que \(\hat {\theta }\) é um estimador não viciado para \(\theta \) (lembre-se que os valores de \(z_i\) são constantes)
obtenha um intervalo de aproximadamente 95% de confiança para \(\theta \) baseado na distribuição assintótica de \(\hat {\theta }\)