# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo:
#
# Integracao Monte Carlo 
#
# Calculando a preditiva de (Y|theta) ~ N(theta,1) e theta ~ N(0,4)
#
# Calculamos a preditiva em uma grade de pontos y
# Para cada y, tomaremos uma amostra de theta
#
#==============================================================================
set.seed(12345) # Semente

m      <- 100                    # Tamanho da grade de valores
y      <- seq(-8,8,length=m)     # Grade de valores y
n      <- 1000                   # Tamanho da amostra 

#------------------------------------------------------------------------------
f1     <- matrix(NA,m,n)         # Primeira maneira
theta1 <- rnorm(n,0,2)           # Primeira maneira
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
f2     <- matrix(NA,m,n) # Segunda maneira
#------------------------------------------------------------------------------

for (i in 1:100){
	# Valores iguais de z para cada y[i]
	f1[i,] <- dnorm(y[i],theta1,1) # Primeira maneira
	#----------------------------------------------------------------------
	# Valores diferentes de z para cada y[i]
	theta2 <- rnorm(n,0,2)         # Segunda maneira
	f2[i,] <- dnorm(y[i],theta2,1) # Segunda maneira
}

#==============================================================================
# Obtendo as estimativas
#==============================================================================
f1y.chapeu <- apply(f1,1,"mean")
f2y.chapeu <- apply(f2,1,"mean")

#==============================================================================
# Ilustracao das aproximacoes
#==============================================================================
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,bty="n")
plot(y,f1y.chapeu,t="l",lwd=3,xlab=expression(y),ylab=expression(f(y)),
    ylim=c(0,0.20))
lines(y,dnorm(y,0,sd=sqrt(5)),lty=2,lwd=3,col="red")
lines(y,f2y.chapeu,col="blue",lwd=3)
legend("topleft",c("Verdadeiro","Monte Carlo 1","Monte Carlo 2"),    
    lty=c(2,1,2),lwd=3,col=c("red","black","blue"),bty="n",cex=1.4)

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