#==============================================================================
# Intervalo de confianca (bilateral) para a media MU quando a variancia 
# SIGMA^2 e' desconhecida (um pequeno estudo de Monte Carlo)
#
#==============================================================================
# Declarando valores para o estudo
#==============================================================================
set.seed(5438)
M      <- 100000         # numero de replicas
n      <- 15             # tamanho da amostra em cada replica
x      <- matrix(NA,M,n) # para salvar as amostras
mu     <- 5              # valor verdadeiro (neste estudo)
sigma2 <- 9              # valor verdadeiro (neste estudo)
sigma  <- sqrt(sigma2)   # valor verdadeiro (neste estudo)

#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rnorm(n,mu,sigma)
}

#==============================================================================
# Calculando a variancia amostral das M amostras
# Esse estimador "var" e' o S^2
#==============================================================================
xvar <- apply(x,1,"var") # vetor de tamanho M
print(length(xvar))

[1] 100000

#==============================================================================
# Calculando os intervalos de confiancas obervados em cada amostra
#==============================================================================
gammax  <- 0.95
gamma2x <- (1+gammax)/2
gamma1x <- 1-gamma2x
liminf  <- (n-1)*xvar/qchisq(gamma2x,n-1) 
limsup  <- (n-1)*xvar/qchisq(gamma1x,n-1)
IC      <- cbind(liminf,limsup)

#==============================================================================
# Proporcao dos ICs que contem o valor verdadeiro de mu
#==============================================================================
print(mean(as.numeric(IC[,1] < sigma2 & sigma2 < IC[,2])))

[1] 0.94937

#==============================================================================
# Visualizacao dos diversos ICs.
#==============================================================================
Mcorte <- 100

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(IC[1,],rep(1,2),type="l",lwd=2,col=2,xlim=c(0,50),
ylim=c(0,Mcorte),xlab="Intervalo de Confiança",ylab="Amostra")
for(i in 2:Mcorte)
    lines(IC[i,],rep(i,2),lwd=2,col=i)
lines(rep(sigma2,2),c(1,Mcorte),lwd=3,col=1)