Capítulo 11: Estimação

Exemplo 11.9

tab11_1<- data.frame(
  x=c(1.2,1.5,1.7,2,2.6),
  y=c(3.9,4.7,5.6,5.8,7.0))

mean(tab11_1$x)

[1] 1.8

mean(tab11_1$y)

[1] 5.4

plot(tab11_1$x,tab11_1$y,pch=16,col= "darkblue")

tab11_1$x_3=3*tab11_1$x
tab11_1$y_3x=tab11_1$y - tab11_1$x_3
tab11_1$y_3x_2=tab11_1$y_3x^2
print(kable(tab11_1,col.names = c("$X$","$Y$", "$3X$","$Y-3X$","$(Y-3X)^2$"), caption = "**Tabela 11.1**: Análise do modelo $Ŷ=3X$"))

**Tabela 11.1**: Análise do modelo $Ŷ=3X$
$X$	$Y$	$3X$	$Y-3X$	$(Y-3X)^2$
1.2	3.9	3.6	0.3	0.09
1.5	4.7	4.5	0.2	0.04
1.7	5.6	5.1	0.5	0.25
2.0	5.8	6.0	-0.2	0.04
2.6	7.0	7.8	-0.8	0.64

## [1] "Total da soma de quadrados dos resíduos(sum(tab11_1$y_3x_2)): 1.06"

Para calcular $\hat{\theta}_{MQ}$, basta calcularmos a soma do produto de X e Y e dividir pela soma de quadrados de X:

sum(tab11_1$x*tab11_1$y)/sum(tab11_1$x^2)

## [1] 2.9441

Exemplo 11.11

Precisamos estimar o valor de $\alpha$ no modelo exponencial. A partir da maximização da função de verossimilhança, a função mle da biblioteca stats4 estima o parâmetro do modelo desejado utilizando os dados de uma amostra coletada.

No contexto do Exemplo 11.11, temos uma distribuição exponencial de parâmetro desconhecido. Vamos então simular uma amostra de uma distribuição exponencial(1/2) e avaliar se o método da função mle estima o parâmetro utilizado na simulação:

x_sim <- rexp(1000, rate = 2) # 1000 observações amostradas de um modelo exponencial de média 1/2
mean(x_sim)

## [1] 0.51367

hist(x_sim, col="lightblue3", border="white", freq = FALSE)

Agora, precisamos da função que será utilizada na otimização da verossimilhança. Como estamos supondo um modelo exmponencial, teremos:

MV <- function(rate) {
     R = dnorm(x_sim, rate)
     return(-sum(log(R))) # devolvendo '-log(verossimilhanca)'
 }

Finalmente, para executar a otimização, façamos:

library(stats4)
emv_rate_x <- mle(MV,  # MV é a função que se deseja otimizar
                  start = list(rate=1)) # start é uma lista com 'chutes' iniciais para o processo de otimização.
emv_rate_x

## 
## Call:
## mle(minuslogl = MV, start = list(rate = 1))
## 
## Coefficients:
##    rate 
## 0.51367

Logo, a estimativa do estimador de máxima verossimilhança de $\alpha$ na amostra selecionada é 0.51367, o que é bastante próximo do valor simulado, 1/2. Note também que este valor é equivalente ao valor da média dos $x_n$ simulados: $\frac{\sum_{n=1}^1000 x_n}{1000}$ = 0.51367.

Exemplo 11.12

Vamos simular 20 amostras de uma $N(5,3^2)$ e construir o Intervalo de confiança(IC) para a média em cada uma delas:

samples=20 # Numero de amostras
N=25 # Numero de observações por amostra

x_sim<-matrix(NA,N,samples) # x_sim armazenara as amostras simuladas
for( i in 1:samples){
  x_sim[,i]<-rnorm(n = N, mean= 5, sd=3) # Simulando 20 N(5,3^2)
}

Os ICs da média para estas amostras são:

cont=0
for( i in 1:samples){
  llim=mean(x_sim[,i])-1.96*sd(x_sim[,i])/sqrt(N) # limite inferior do IC da i-ésima amostra
  ulim=mean(x_sim[,i])+1.96*sd(x_sim[,i])/sqrt(N) # limite inferior do IC da i-ésima amostra
  if (llim>5 || ulim<5){
    line_col="red"
    cont=cont+1
  }
  else line_col="black"
  if(i!=1)
    par(new=T)
  plot(x=i,y=mean(x_sim[,i]), xlim=c(0,21), ylim=c(0,10), pch=20,cex=1, col=line_col, ylab="Médias", xlab="Amostras") # plotando as médias de cada amostra
  lines(x = c(i,i), y =c(llim,ulim), lwd=1, col=line_col) # Adicionando as linhas dos ICs de cada amostra
  abline(h=5,lty = 'dashed', col="blue") # Valor real da média 
}

Figura 11.4: Intervalos de confiança para a média de uma N(5,9), para 20 amostras de tamanho n=25.

Note que 2 intervalos não contiveram o valor real da média, isto é 8% dos casos.

Para os exemplos 11.13,11.15 e 11.16, vamos construir uma função que, dados a média, desvio padrão e tamanho da amostra coletada e o nível de confiança, retorna o intervalo de confiança desejado:

IC<-function(n, media, desvpad, gama){
  left_tail<-(1-gama)/2
  gamma_adj<-left_tail+gama # ajuste para somar a cauda inferior da distribuição somente para encontrar o quantil gama usando a função dist. acumulada
  llim<-media-qnorm(gamma_adj)*desvpad/sqrt(n) # calcula o limite inferior
  ulim<-media+qnorm(gamma_adj)*desvpad/sqrt(n) # calcula o limite inferior
  writeLines(paste("IC(",gama,"):[",round(llim,3),";",round(ulim,3),"]\n",sep=""))
  return(ic_calc=list(inf.lim=llim,upper.lim=ulim))
}

Exemplo 11.13

Usando a função que acabamos de construir, teremos:

IC(n=25, media = 485, desvpad=sqrt(100), gama=0.95)

## IC(0.95):[481.08;488.92]

## $inf.lim
## [1] 481.08
## 
## $upper.lim
## [1] 488.92

Exemplo 11.15

Note que construímos a função IC para calcular a variância do estimador internamente. Em outras palavras, a variância que devemos passar para a função é a variância amostral $p(1-p)=1/4$(no pior caso), e não a variância do estimador. Assim:

IC(n=400, media = 0.6, desvpad=sqrt(1/4), gama=0.95)

## IC(0.95):[0.551;0.649]

## $inf.lim
## [1] 0.551
## 
## $upper.lim
## [1] 0.649

Exemplo 11.16

Novamente, utilizando a função IC:

p=80/400
dp=sqrt(p*(1-p))
IC(n=400, media = p, desvpad=dp, gama=0.9)

## IC(0.9):[0.167;0.233]

## $inf.lim
## [1] 0.1671
## 
## $upper.lim
## [1] 0.2329

Se quisermos o intervalo conservador, temos que p=1/2 para o cáculo da variância:

p=1/2
dp=sqrt(p*(1-p))
IC(n=400, media = 80/400, desvpad=dp, gama=0.9)

## IC(0.9):[0.159;0.241]

## $inf.lim
## [1] 0.15888
## 
## $upper.lim
## [1] 0.24112

11.9 Exemplos Computacionais

Bootstrap

O bootstrap consiste em reamostrar com reposição a amostra que temos um número b de vezes. Calcula-se o valor do estimador em cada uma destas amostras e assim obtém-se uma distribuição empírica do estimador de interesse, como fizemos isso no capítulo anterior. Portanto, teremos:

amostra_original<-rexp(50, rate=3) # Vamos simular uma exponencial de média 1/3
n=length(amostra_original)
md=numeric() # vetor que armazenará as medianas calculadas ao longo das amostras bootstrap
b=30 # Número de amostras bootstrap a serem resorteadas baseadas na amostra original.
for( i in 1:b){ # Vamos calcular 30 valores da mediana para construir sua dist. empirica
  boot<-sample(amostra_original, size = n, replace = TRUE) # Tomando amostras de tamanho n, com reposicao
  md[i]<-median(boot)
}
hist(md, col="lightblue3", border= "white")

Agora vamos estudar o erro padrão da mediana neste exemplo. Note que o erro padrão do estimador depende do número de amostras bootstrap calcumamos. Assim, vamos calcular o erro padrão do estimador quando b={30,50,100,500,1000} e verificar qual o seu comportamento:

ep_md=numeric() # vetor que armazenará os erros padrões da mediana para cada rodada do bootstrap
md=numeric() # vetor que armazenará as medianas calculadas ao longo das amostras bootstrap
b=c(30,50,100,500,1000) # Número de amostras bootstrap a serem resorteadas baseadas na amostra original.
for (j in 1:5){ # neste for, percorreremos os valores de b
  for( i in 1:b[j]){ # Vamos calcular 30 valores da mediana para construir sua dist. empirica
    boot<-sample(amostra_original, size = n, replace = TRUE) # Tomando amostras de tamanho n, com reposicao
    md[i]<-median(boot)
  }
  ep_md[j]=sd(md)/b[j]
}

plot(1:5, ep_md, cex=2, type="b", ylab="Erro Padrão", xlab="b: tamanho do bootstrap")

Como era esperado, o erro padrão diminui a medida que aumentamos o número de amostras bootstrap.

Exemplo 11.19

Utilizando o mesmo procedimento, vamos recalcular o resultado deste exemplo:

amostra_original<-c(2,5,3,4,6)
n=length(amostra_original)
b=5 

boot<-matrix(NA,n,b) # armazena as amostras geradas
md=numeric() 
xbarr=numeric()

for( i in 1:b){ 
  boot[,i]<-sample(amostra_original, size = n, replace = TRUE) 
  md[i]<-median(boot[,i])
  xbarr[i]<-mean(boot[,i])
}

**Tabela 11.2**: Procedimento *bootstrap*.
Bootstrap:x1*	Bootstrap:x2*	Bootstrap:x3*	Bootstrap:x4*	Bootstrap:x5*	md(x*)	média(x*)
4	5	5	5	2	4	4.4
6	5	5	2	3	4	3.8
2	2	6	5	5	5	4.6
4	3	3	5	5	5	3.8
6	4	4	2	3	3	3.6

O erro padrão da mediana é, portanto:

### Calculando agora o EP para a mediana a partir das amostras geradas:  
ep=0
for( i in 1:b){
  ep=ep+(md[i]-mean(md))^2/(b-1)
}
ep=sqrt(ep)
ep

## [1] 0.83666

Exemplo 11.20

Da mesma forma, o erro padrão da média é:

### Calculando agora o EP para a mediana a partir das amostras geradas:  
ep_media=0
for( i in 1:b){
  ep_media=ep_media+(xbarr[i]-mean(xbarr))^2/(b-1)
}
ep_media=sqrt(ep_media)
ep_media

## [1] 0.43359

O IC para a mediana com coeficiente de confiança 95% é:

IC(n = 1,media = 5.2, desvpad = ep, gama = 0.95)

## IC(0.95):[3.56;6.84]

## $inf.lim
## [1] 3.5602
## 
## $upper.lim
## [1] 6.8398

Note que, como já fizemos a conta do erro padrão fora da função IC, inserimos n=1, pois a divisão por $\sqrt{n}$ já foi feita.

OBS: Eventualmente, como b=5, os procedimento apresentados nos Exemplos 11.19 e 11.20 podem resultar em valores diferentes dos valores apresentados no livro. A medida que aumentamos o valor deb, alcançamos valores mais consistentes. Se, por exemplo, utilizarmos b=50 com a mesma amostra original, obteríamos:

amostra_original<-c(2,5,3,4,6)
n=length(amostra_original)
b=50 

boot<-matrix(NA,n,b) # armazena as amostras geradas
md=numeric() 
xbarr=numeric()

for( i in 1:b){ 
  boot[,i]<-sample(amostra_original, size = n, replace = TRUE) 
  md[i]<-median(boot[,i])
  xbarr[i]<-mean(boot[,i])
}

### Calculando o EP para a mediana a partir das 50 amostras geradas:  
ep=0
for( i in 1:b){
  ep=ep+(md[i]-mean(md))^2/(b-1)
}
ep=sqrt(ep)
ep

## [1] 0.88086

### Calculando o EP para a media a partir das50 amostras geradas:  
ep_media=0
for( i in 1:b){
  ep_media=ep_media+(xbarr[i]-mean(xbarr))^2/(b-1)
}
ep_media=sqrt(ep_media)
ep_media

## [1] 0.60117

Há inúmeros pacotes que contêm versões mais sofisticadas do algoritmo de bootstrap. Dois pacotes bastante difundidos são o pacote bootstrap e o pacote boot. Apresentaremos neste texto a função bootstrap() do pacote bootstrap, que constrói amostras de bootstrap de acordo com a seguinte sintaxe:

bootstrap(x,nboot,theta,..., func=NULL)
Em que:
x:     Vetor contendo os dados a serem reamostrados
nboot: numero de amostras bootstrap desejadas  
theta: Função a ser calculada em cada uma das amostras. Deve ter o argumento x.
func:  (opcional) argumento adicional indicando a distribuição funcional do estimador de theta 
       desejada. É utilizada pelo algoritmo para gerar as estimatias de jackknife do erro-padrão.

Teríamos então, nos Exemplos 11.19 e 11.20 os seguintes resultados:

library(bootstrap)
amostra_original<-c(2,5,3,4,6)
### Exemplo 11.19: Mediana
boot_11_19<-bootstrap(x=amostra_original, nboot=5, theta=function(x){median(x)})
boot_11_19$thetastar

## [1] 3 5 6 3 3

### Exemplo 11.19: Mediana
boot_11_20<-bootstrap(x=amostra_original, nboot=5, theta=function(x){mean(x)})
boot_11_20$thetastar

## [1] 3.8 4.4 3.6 2.6 4.2

Vemos então que os objetos boot_11_19$thetastar e boot_11_20$thetastar armazenam os resultados equivalentes às colunas “md(x*)” e “média(x*)” da Tabela 11.2.

Se desejarmos, por exemplo, a estimativa do erro padrão baseada no método jackknife:

#Se desejarmos, como nos exemplos, a estimativa do erro padrão, faríamos:
err_pad<- function(x){sqrt(var(x)/length(x))}
bootstrap(x=amostra_original, nboot=15, theta=function(x){median(x)}, func=err_pad)

## $thetastar
##  [1] 4 5 3 3 5 5 3 3 5 3 4 4 5 5 4
## 
## $func.thetastar
## [1] 0.22817
## 
## $jack.boot.val
## [1] 0.25000 0.28868 0.33333 0.66667 0.50000
## 
## $jack.boot.se
## [1] 0.31003
## 
## $call
## bootstrap(x = amostra_original, nboot = 15, theta = function(x) {
##     median(x)
## }, func = err_pad)

bootstrap(x=amostra_original, nboot=15, theta=function(x){mean(x)}, func=err_pad)

## $thetastar
##  [1] 2.8 4.0 4.0 4.4 4.0 4.0 3.6 4.8 4.0 3.8 4.0 3.6 4.2 3.2 4.0
## 
## $func.thetastar
## [1] 0.12168
## 
## $jack.boot.val
## [1] 0.149666 0.100000 0.068313 0.233238 0.271293
## 
## $jack.boot.se
## [1] 0.1543
## 
## $call
## bootstrap(x = amostra_original, nboot = 15, theta = function(x) {
##     mean(x)
## }, func = err_pad)

Para opções mais sofisticadas de bootstrap, veja o pacote boot, que contém opções de reamostragem ponderada, estratificada, cálculo de intervalos de confiança, dentre outras. Veja help.search("boot") para maiores detalhes sobre este pacote. Referências complementares sobre este assunto pode ser:

Efron, B. and Tibshirani, R. J. (1993). An Introduction to the Bootstrap, New York: Chapman & Hall.
Davison, A. C. and Hinkley, D. V. (1997). Bootstrap Methods and their Application, Cambridge University Press.

Capítulo Anterior | Introdução | Próximo Capítulo

Estatística Básica - Bussab e Morettin (2017)

Comandos R para análises estatísticas