####### Exemplo construção de intervalos de confiança ###### x <- rnorm(100, mean=0, sd=1) ## Construir um intervalo de confiança para a média cont <- qnorm(0.975) mu.u <- mean(x) + (cont/length(x)) mu.l <- mean(x) - (cont/length(x)) intervalo <- c(mu.l, mean(x) , mu.u) ## Interpretação do intervalo de confiança # Se eu retirar repetidas amostras e para cada amostra calcular o mesmo intervalo de confiança espero que em 95% deste intervalos # o verdadeiro valor do parâmetro esteja contido. intervalo <- matrix(ncol=2,nrow=100) cont <- qnorm(0.975) popu <- rnorm(10000, mean=0, sd=1) for(i in 1:100){ amostra <- sample(1:10000, 100) x <- popu[amostra] mu.u <- mean(x) + (cont/sqrt(length(x))) mu.l <- mean(x) - (cont/sqrt(length(x))) intervalo[i,] <- c(mu.l,mu.u) } plot(range(intervalo)~c(1,100),type="n") abline(=0) for(i in 1:100){ Sys.sleep(0.25) arrows(c(i),intervalo[i,1],c(i),intervalo[i,2],code=3,angle=90,length=0.03,col=ifelse(intervalo[i,1]*intervalo[i,2] > 0, "red","black")) } ### Distribuição assintótica do EMV estimativa <- c() for(i in 1:1000){ amostra <- sample(1:10000, 100) x <- popu[amostra] estimativa[i] <- mean(x) } hist(estimativa,prob=T,ylim=c(0,4.5),xlab="Estimativas",ylab="Densidade",main="Distribuição empírica") curve(dnorm(x,mean=mean(estimativa),sd=sd(estimativa)),add=TRUE) curve(dnorm(x,mean=0,sd=1/sqrt(100)),add=TRUE,col="black",lty=2) lines(density(estimativa),col="blue",lty=3) mean(estimativa) sd(estimativa) #### E se fosse Poisson ##### popu = rpois(10000, lambda=1) intervalo <- matrix(ncol=2,nrow=100) cont <- qnorm(0.975) for(i in 1:100){ amostra <- sample(1:10000, 100) x <- popu[amostra] mu.u <- mean(x) + cont * mean(x)/sqrt(sum(x)) mu.l <- mean(x) - cont * mean(x)/sqrt(sum(x)) intervalo[i,] <- c(mu.l,mu.u) } plot(range(intervalo)~c(1,100),type="n") abline(h=1) for(i in 1:100){ Sys.sleep(0.25) arrows(c(i),intervalo[i,1],c(i),intervalo[i,2],code=3,angle=90,length=0.03,col=ifelse(intervalo[i,1] > 1 | intervalo[i,2] < 1, "red","black")) } estimativa <- c() for(i in 1:1000){ amostra <- sample(1:10000, 100) x <- popu[amostra] estimativa[i] <- mean(x) } hist(estimativa,prob=T,xlab="Estimativas",ylab="Densidade",main="Distribuição empírica") curve(dnorm(x,mean=mean(estimativa),sd=sd(estimativa)),add=TRUE) curve(dnorm(x,mean=1,sd=mean(x)/sqrt(sum(x))),add=TRUE,col="black",lty=2) lines(density(estimativa),col="blue",lty=3) mean(estimativa) sd(estimativa)