#============================================================
# Inferência Estatística: Distribuição amostral da média
#                                                (10/05/2011)
#
#                                 Professor Walmes M. Zeviani
#                                      www.leg.ufpr.br/ce003b
#============================================================

#------------------------------------------------------------
# o que é parâmetro da população?
# o que é estimador de um parâmetro?
# o que é estimativa de uma parâmetro?

require(MASS)

#------------------------------------------------------------
# distribuição amostral de uma estatística
#  fica completamente descrita quando conhecemos os
#  resultados assumiados pela estatística em ***todas as
#  amostras possíveis*** e suas probabilidades de ocorrência
#  associadas

#------------------------------------------------------------
# experimento aleatório de retirar e papeis com os números
# 0, 3, e 12 grafados. A v.a. é o número do papel sorteado

x <- c(0,3,12)       # valores assumidos pela v.a.
px <- c(1/3,1/3,1/3) # suas probabilidades
plot(px~x, type="h", pch=19, ylim=c(0,0.5))
points(x,px,pch=19)

#------------------------------------------------------------
# valor esperado da variável aleatória

x*px
mu <- sum(x*px); mu    # E(X) = mu
si <- sum((x-mu)^2*px) # Var(X) = sigma^2
abline(h=0)
points(sum(x*px), 0, cex=2, pch=17)

#------------------------------------------------------------
# espaço amostral para a retirada de 3 papeis

universo <- expand.grid(e1=x, e2=x, e3=x)
universo

#------------------------------------------------------------
# valores do estimador A (média amostral) em cada amostra

universo$A <- apply(universo, 1,
                    function(amostra) mean(amostra))
universo
sort(unique(universo$A))

#------------------------------------------------------------
# valores do estimador B (mediana amostral) em cada amostra

universo$B <- apply(universo, 1,
                    function(amostra) median(amostra[1:3]))
universo
sort(unique(universo$B))

#------------------------------------------------------------
# qual a distribuição do estimador A?

A <- sort(unique(universo$A))             # valores assumidos
pA <- c(table(universo$A)/nrow(universo)) # probabilidades
fractions(pA)

plot(pA~A, type="h", ylim=c(0,max(pA))); points(A, pA, pch=19)

#------------------------------------------------------------
# qual a distribuição do estimador B?

B <- sort(unique(universo$B))             # valores assumidos
pB <- c(table(universo$B)/nrow(universo)) # probabilidades
fractions(pB)

plot(pB~B, type="h", ylim=c(0,max(pB))); points(B, pB, pch=19)

#------------------------------------------------------------
# baseado no experimento e na distribuição amostral de A
# qual P(A)=0, P(A)=1, P(A)=5?
 
pA

#------------------------------------------------------------
# qual estimador é melhor para a média, A ou B?
# que características um estimador deve ter pra ser bom?

mA <- sum(A*pA)
mA                       # não viciado

mB <- sum(B*pB)
mB                       # viciado viciado

vA <- sum((A-mA)^2*pA)
vA                       # variancia pequena

vB <- sum((B-mB)^2*pB)
vB                       # variancia grande

#------------------------------------------------------------
# e se não conseguissemos escrever o espaço amostral do
# experimento para atribuir probabilidade e chegar na
# distribuição amostral de A ou B?
# vamos aumentar o número de bolas na caixa e retiradas?

#------------------------------------------------------------
# e se o experimento possuir características que permitam
# o uso de uma função de distribuição de probabilidades, ou
# seja, um modelo conhecido como binomial, poisson, normal?

#------------------------------------------------------------
# vamos estudar a distribuição amostral dos estimadores
# da variância de uma população: com denominador "n" e "n-1"

#------------------------------------------------------------
# abaixo uma amostra da altura (m) de 12 pessoas de uma
# população com mu=1.75 e sigma^2=0.1 (normal)

amostra <- rnorm(12, mean=1.75, sd=sqrt(0.1))
amostra

#------------------------------------------------------------
# variância com denominador "n" (populacional)

sum((amostra-1.75)^2)/(12)

#------------------------------------------------------------
# variância com denominador "n-1" (amostral)

sum((amostra-1.75)^2)/(12-1)

#------------------------------------------------------------
# não temos como saber quem é melhor olhando apenas uma
# amostra pois sabemos que de amostra para amostra há
# diferenças. então temos que avaliar o comportamento
# considerando infinitas amostras

nexp <- 30
s.pop <- c()
s.amo <- c()
k <- 0.5; op <- par(mar=c(3,3,1,1), fig=c(0,1,0,k))
plot(c(0.01,0.23), c(0,nexp), xlim=c(0.02,0.28), type="n")
for(i in 1:nexp){
  amostra <- rnorm(12, mean=1.75, sd=sqrt(0.1))
  s.pop[i] <- sum((amostra-1.75)^2)/(12)
  s.amo[i] <- sum((amostra-1.75)^2)/(12-1)
  points(c(i,i)~c(s.pop[i], s.amo[i]), col=1:2)
  if(i<15) Sys.sleep(0.5)
}
par(fig=c(0,1,k,1), new=TRUE)
plot(density(s.pop), xlim=c(0.02,0.28), main=NA)
lines(density(s.amo), col=2)
abline(v=c(mean(s.pop), mean(s.amo)), col=c(1,2),
       lwd=c(1,1))
par(op)

#------------------------------------------------------------
# é possivel que a distribuição de probabilidades um
# estimador seja uma distribuição de probabilidades
# conhecida? qual a distribuição da média amostral?

#------------------------------------------------------------
# vamos amostrar dados da seguinte distribuição

dev.off()

x <- 0:25                  # valores assumidos
px <- dpois(x, lambda=10)  # probabilidades
plot(px~x, type="h"); points(x, px, pch=19)

rpois(11, lambda=10)

nexp <- 500
xb <- c()
for(i in 1:nexp){
  amostra <- rpois(11, lambda=10)
  #amostra <- rexp(11, rate=1)
  xb[i] <- mean(amostra)
}
plot(density(xb), col=1)
hist(xb, freq=FALSE)

#------------------------------------------------------------
# existe uma distribuição de probabilidades que descreva bem
# a distribuição da média amostral? TLC!

xbb <- mean(xb)
vxbb <- var(xb)

curve(dnorm(x, mean=xbb, sd=sqrt(vxbb)), add=TRUE,
      col="red", lwd=2)

#------------------------------------------------------------
# qual seria a média a a variância da distribuição das
# médias?

xbb  # quando o estimador é não viciado é igual ao parâmetro
vxbb # é a variância dos dados/número de elementos da amostra
     # a variância das estimativas é chamada de erro padrão!

#------------------------------------------------------------
# de que depende a boa aproximação da distribuição normal
# para representar a distribuição das médias?


par(mfrow=c(3,3), mar=c(3,2,2,1))
for(n in c(1,2,3,5,10,50,90,120,200)){
# xb <- apply(matrix(runif(2000*n), ncol=n), 1, mean); hist(xb, main=n, freq=FALSE, col="gray90")
# xb <- apply(matrix(rexp(2000*n, rate=1), ncol=n), 1, mean); hist(xb, main=n, freq=FALSE, col="gray90")
 xb <- apply(matrix(rbeta(2000*n, 1/5.1, 1/5.1), ncol=n), 1, mean); hist(xb, main=n, freq=FALSE, col="gray90")
 curve(dnorm(x, mean(xb), sd(xb)), col=3, add=TRUE, lwd=2)
 lines(density(xb), col=2, lwd=2); box()
}

layout(1)
curve(dbeta(x, 1/2.1,1/2.1))

#------------------------------------------------------------
# a aproximação não depende só do tamanho da amostra mas da
# forma da distribuição dos dados originais

#------------------------------------------------------------
# estimação por intervalo
#  o que é? é melhor porque? o que o intervalo representa?

layout(1)
n <- 7
nexp <- 100
int <- matrix(0, ncol=3, nrow=nexp)
plot(c(1.2,2.3)~c(0,nexp), xlim=c(0,nexp), type="n",
     xlab="número do experimeto", ylab="Estimativa intervalar (95%)")
abline(h=1.75)
for(i in 1:nexp){
  amostra <- rnorm(n, mean=1.75, sd=sqrt(0.1))
  int[i,] <- mean(amostra)+c(-1,0,1)*1.96*sqrt(0.1)/sqrt(n)
  if(i<15) { Sys.sleep(0.5) }
  col <- ifelse(int[i,1]>1.75 | int[i,3]<1.75, 2, 1)
  points(int[i,2]~i, col=col)
  arrows(i,int[i,1],i,int[i,3], angle=90, code=3,
         length=0.05, col=col)
}

#------------------------------------------------------------
# como obter o intervalo de confiança para a média?
#  usando percentis da distribuição amostral do estimador!