Aula de 15/06/2009

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                    #### 
######################################################

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)
##