#==============================================================================
# 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 media amostral das M amostras
#==============================================================================
xbar <- apply(x,1,"mean") # vetor de tamanho M
print(length(xbar))

[1] 100000

#==============================================================================
# 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 o desvio padrao amostral das M amostras
#==============================================================================
xdesvpad <- sqrt(xvar) # vetor de tamanho M
print(length(xdesvpad))

[1] 100000

#==============================================================================
# Calculando os intervalos de confiancas obervados em cada amostra
#==============================================================================
gammax <- 0.95
liminf <- xbar - qt((1+gammax)/2,n-1)*xdesvpad/sqrt(n)
limsup <- xbar + qt((1+gammax)/2,n-1)*xdesvpad/sqrt(n)
IC     <- cbind(liminf,limsup)

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

[1] 0.94918

#==============================================================================
# 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(-1,11),
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(mu,2),c(1,Mcorte),lwd=3,col=1)