# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo:
#
# Integracao Monte Carlo
# Rao-Blackwell
#
# Calculando E(exp(X^2)) para X ~ t_nu(0,1)
#
# rm(list=ls())
#==============================================================================
set.seed(12345)         # Semente
nu   <- 5               # Grau de liberdade
r    <- 1000            # Numero de replicoes
n    <- seq(20,1000,20) # Tamanhos de amostra
m    <- length(n)       # Tamanho da grade
I1   <- matrix(NA,m,r)
I3   <- matrix(NA,m,r)
for(i in 1:m){      # Diferentes tamanhos de amostra
    for(j in 1:r){  # Replicações
		x       <- rt(n[i],nu)
		I1[i,j] <- mean(exp(-(x^2)))
		y       <- rgamma(n[i],nu/2,nu/2)
		I3[i,j] <- mean((2/y+1)^(-1/2))
	}
}

#==============================================================================
# Estatisticas do estudo de Monte Carlo
#==============================================================================
media1  <- apply(I1,1,"mean")
media2  <- apply(I3,1,"mean")
ic.inf1 <- apply(I1,1,"quantile",probs=0.025)
ic.inf2 <- apply(I3,1,"quantile",probs=0.025)
ic.sup1 <- apply(I1,1,"quantile",probs=0.975)
ic.sup2 <- apply(I3,1,"quantile",probs=0.975)

#==============================================================================
# Grafico para comparacao
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,
    lab=c(6,8,5),mar=c(4.5,5,1,1),xpd=T,cex.main=2.0,bty="n")
plot(n,media1,type="l",lwd=4,col="black",xlab=expression(n),
  ylab=expression(I),ylim=c(0.39,0.69))
lines(n,media2,lwd=4,col="red")
lines(n,ic.inf1,lwd=4,col="black",lty=2)
lines(n,ic.sup1,lwd=4,col="black",lty=2)
lines(n,ic.inf2,lwd=4,col="red",lty=2)
lines(n,ic.sup2,lwd=4,col="red",lty=2)
legend("topright",legend=c("Simples","Rao-Blackwell"),
  col=c("black","red"),lty=c(1,1),lwd=c(4,4),bty="n",cex=1.5)
legend("bottomright",
       legend=c("IC 95%: Simples","IC 95%: Rao-Blackwell"),
       col=c("black","red"),lty=c(2,2),lwd=c(4,4),bty="n",cex=1.5)

# rm(list=ls()) 
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================