====== Aula de 15/06/2009 ======

===== Ilustração computacional de conceitos de inferência =====

==== Exemplo 1: Distribuição Amostral de uma proporção====

<code R>
set.seed(37)

p <- runif(1)
n <- 1000

#p <- 0.025
#n <- 200

## retirando uma amostra de tamanho n e calculando a estimativa de p
am1 <- sample(c(0,1), size=n, prob=c(1-p, p), rep=T)
am1
p1 <- mean(am1)
p1
## IC (assintótico)
p1 + qnorm(c(0.025, 0.975))*sqrt(p1*(1-p1)/n)

## outra amostra...
am2 <- sample(c(0,1), size=n, prob=c(1-p, p), rep=T)
am2
p2 <- mean(am2)
p2
p2 + qnorm(c(0.025, 0.975))*sqrt(p2*(1-p2)/n)

## e mais uma amostra
am3 <- sample(c(0,1), size=n, prob=c(1-p, p), rep=T)
am3
p3 <- mean(am3)
p3
p3 + qnorm(c(0.025, 0.975))*sqrt(p3*(1-p3)/n)

## agora retirando 10,000 amostras de tamanho n
ams <- lapply(1:10000, function(x) sample(c(0,1), size=n, prob=c(1-p, p), rep=T))

## calculando as estimativas das proporções para cada uma das 10,000 amostras
p.sim <- sapply(ams, mean)
length(p.sim)
p.sim
summary(p.sim)

## comparando as distribuições empírica e teórica:
## histogramas das estimativas
hist(p.sim, prob=T)
lines(density(p.sim))
## distribuição amostral (teórica) do estimador
vals <- seq(0.48, 0.62, leng=200)
lines(vals, dnorm(vals, m=p, sd=sqrt(p*(1-p)/n)), col=2, lwd=2)

## "margem de erro" teórica
quantile(p.sim, prob=c(0.025, 0.975))
diff(quantile(p.sim, prob=c(0.025, 0.975)))/2

## IC's para todas as 10,000 amostras
ic.sim <- sapply(ams, function(x){p.est <-  mean(x);
                                  p.est + qnorm(c(0.025, 0.975))*sqrt(p.est*(1-p.est)/n)}
                 )
dim(ic.sim)
ic.sim[,1:8]
sum(ic.sim[1,] <= 0.5)
sum(ic.sim[2,] <= 0.5)

dentro <- apply(ic.sim, 2, function(x) (p> x[1] & p< x[2]))
mean(dentro)

## e qual mesmo é o valor de p verdadeiro?
p
abline(v=p, col=4, lwd=2)

######################################################
####          FIM DO EXEMPLO 1                    #### 
######################################################
</code>

====  Exemplo 2: teste de permutação como uma aproximação ou alternativa a expressões teóricas====

<code R>
am1 <- c(6.6, 10.3, 10.8, 12.9, 9.2, 12.3, 7)
length(am1)
mean(am1)

am2 <- c(8.1, 9.8, 8.7, 10, 8.2, 8.7, 10.1)
length(am2)
mean(am2)

## comparação de variancias por simulação
var.q <- var(am2)/var(am1)
var.q

am1
am2
ams <- lapply(1:10000,
              function(x) matrix(sample(c(am1, am2)), nc=2))

var.sim <- sapply(ams, function(x) var(x[,2])/var(x[,1]))
hist(var.sim, prob=T, ylim=c(0, 0.65))

abline(v=quantile(var.sim, c(0.025, 0.975)))

vals <- seq(0,40, l=200)
lines(vals, df(vals, df1=length(am2)-1, df2=length(am2)-1), col=2, lwd=2)

abline(v=var.q, col=4)
mu.dif <- mean(am2) - mean(am1)
mu.dif

mu.dif.sim <- sapply(ams, function(x) mean(x[,2])-mean(x[,1]))

hist(mu.dif.sim, prob=T)
vals <- seq(-6,6, l=200)
lines(vals, dt(vals, df=9.137), col=2, lwd=2)


abline(v=quantile(mu.dif.sim, prob=c(0.025, 0.975)))
abline(v=mu.dif, col=4)

## resultado usando uma função já pronta do R
t.test(am1, am2, alt="two", mu =0, pair=F, var.equal=F,
       conf=0.95)
##
</code>