# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo para gerar de 
# f(x) = cst*exp(-x^2/2)*(sin(6*x)^2+3*cos(x)^2*sin(4*x)^2+1)
# com cst desconhecida
#==============================================================================
f  <- function(x){
	exp(-x^2/2)*(sin(6*x)^2+3*cos(x)^2*sin(4*x)^2+1)
}

#==============================================================================
# Analise grafica do comportamento da f.d.p.
#==============================================================================
xseq <- seq(-5,5,0.001)
yseq <- f(xseq)
cst  <- sum(yseq)*0.001  # Aproximacao de Riemann
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,6,0),
    mar=c(4.5,5,2,1),bty="n")
plot(xseq,yseq/cst,lwd=3,type="l",xlab=expression(x),xlim=c(-5,5),
     ylim=c(0,0.7),ylab=expression(f(x)),col="darkblue")

#==============================================================================
# Reamostragem ponderada
#==============================================================================
set.seed(54345)  # Semente
n      <- 500000 # Tamanho da amostra das propostas
m      <- 100000 # Reamostragem (tamanho da amostra)

#==============================================================================
# Utilizando como proposta a normal padrao: N(0,1)
#==============================================================================
xprop <- rnorm(n)

#==============================================================================
# Note que calcularemos os pesos abaixo sem a necessidade da constante de
# normalizacao da f.d.p.
#==============================================================================
pesos  <- f(xprop)/dnorm(xprop)
prob   <- pesos/sum(pesos)	
print(sum(prob))
x      <- sample(xprop,m,replace=T,prob=prob)

[1] 1

#==============================================================================
# Histograma
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,6,0),
    mar=c(4.5,5,2,1),bty="n")
hist(x,nclass=100,prob=T,main="",xlab=expression(x),ylab=expression(f(x)),
     xlim=c(-5,5),ylim=c(0,0.7),col="darkgreen")
lines(xseq,yseq/cst,lwd=3,col="red")

#==============================================================================
# Utilizando como proposta Cauchy
#==============================================================================
xprop  <- rt(n,df=1)

#==============================================================================
# Note que calcularemos os pesos abaixo sem a necessidade da 
# constante de normalizacao da f.d.p.
#==============================================================================
pesos  <- f(xprop)/dt(xprop,df=1)
prob   <- pesos/sum(pesos)	
print(sum(prob))
y      <- sample(xprop,m,replace=T,prob=prob)

[1] 1

#==============================================================================
# Histograma
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,6,0),
    mar=c(4.5,5,2,1),bty="n")
hist(y,nclass=100,prob=T,main="",xlab=expression(x),ylab=expression(f(x)),
     xlim=c(-5,5),ylim=c(0,0.7),col="darkgreen")
lines(xseq,yseq/cst,lwd=3,col="red")

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