#==============================================================================
# A distribuicao amostral da media quando a variancia e' desconhecida 
# (um pequeno estudo de Monte Carlo).
#
# A distribuicao t-Student com n-1 graus de liberade.
#==============================================================================
# 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

#==============================================================================
# Analise da distribuicao de X barra dado sigma^2
#==============================================================================
aux    <- 5*sigma/sqrt(n)
xseq   <- seq(mu-aux,mu+aux,length=1001)
yseq   <- dnorm(xseq,mu,sigma/sqrt(n))

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
hist(xbar,prob=T,lwd=2,nclass=30,main="",xlim=c(mu-aux,mu+aux),
  xlab=expression(bar(x)),ylab=expression(f(bar(x))),col=0)
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(mu,2),c(0,max(yseq)),lwd=4,col="red",lty=2)

#==============================================================================
# Analise da distribuicao de S^2 dado sigma^2
#==============================================================================
aux    <- (n-1)+6*sqrt(2*(n-1))  
wseq   <- seq(0,1.1*aux,length=1001)
useq   <- dchisq(wseq,n-1)

chisq  <- (n-1)*xvar/sigma2 # alterando a distribuicao de S^2

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
hist(chisq,prob=T,lwd=2,nclass=30,main="",xlim=c(0,1.1*aux),
  xlab=expression(y),ylab=expression(f(y)),col=0)
lines(wseq,useq,lwd=4,col="blue")

#==============================================================================
# Analise da distribuicao de X barra (sigma^2 desconhecido)
#==============================================================================
aux    <- 4*sigma
aseq   <- seq(-aux,aux,length=10001)
bseq   <- dt(aseq,n-1)

tpad   <- sqrt(n)*(xbar-mu)/xdesvpad

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
hist(tpad,prob=T,lwd=2,nclass=30,main="",xlim=c(-aux,aux),
  xlab=expression(t),ylab=expression(f(t)),col=0)
lines(aseq,bseq,lwd=4,col="blue")

#==============================================================================
# Analise da distribuicao de S^2 dado sigma^2
# O interesse e' so' mostrar S^2 com distribuicao gama.
#==============================================================================
aux    <- sigma2+6*sqrt(2*(sigma2^2)/(n-1))  
wseq   <- seq(0,1.1*aux,length=1001)
a      <- (n-1)/2
b      <- a/sigma2
useq   <- dgamma(wseq,a,b)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
hist(xvar,prob=T,lwd=2,nclass=30,main="",xlim=c(0,1.1*aux),
  xlab=expression(s^2),ylab=expression(f(s^2)),col=0)
lines(wseq,useq,lwd=4,col="blue")
lines(rep(sigma2,2),c(0,max(useq)),lwd=4,col="red",lty=2)