#==============================================================================
# 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 unilaterais obervados 
# em cada amostra
#==============================================================================
gammax <- 0.95
ICinf  <- xbar - qt(gammax,n-1)*xdesvpad/sqrt(n)
ICsup  <- xbar + qt(gammax,n-1)*xdesvpad/sqrt(n)

#==============================================================================
# Proporcao dos ICs que contem o valor verdadeiro de mu
#==============================================================================
print(mean(as.numeric(ICinf < mu))) # Para o IC inferior: (a(x),infinito)
print(mean(as.numeric(ICsup > mu))) # Para o IC superior: (-infinito,b(x))

[1] 0.95007
[1] 0.94963

#==============================================================================
# Visualizacao dos diversos IC Inferior
#==============================================================================
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(c(ICinf[1],11),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(c(ICinf[i],11),rep(i,2),lwd=2,col=i)
lines(rep(mu,2),c(1,Mcorte),lwd=3,col=1)

#==============================================================================
# Visualizacao dos diversos IC Superior
#==============================================================================
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(c(-1,ICsup[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(c(-1,ICsup[i]),rep(i,2),lwd=2,col=i)
lines(rep(mu,2),c(1,Mcorte),lwd=3,col=1)