# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo para gerar de f(x) = cst*exp( -(x^4)/2 ) pelo metodo da rejeicao
# utilizando propostas N(0,1)
#==============================================================================
set.seed(54345)  # Semente
n      <- 100000 # Tamanho da amostra

#==============================================================================
# Aproximacao de Riemann da constante de integracao; veremos adiante no curso.
# f(x) = cst*exp( -(x^4)/2 )
#==============================================================================
xseq <- seq(-4,4,0.00001)     # Grade de valores
yseq <- exp(-0.5*(xseq^4))    # Avaliacao da funcao
cst  <- 1/(sum(yseq)*0.00001) # Aproximacao de Riemann

#==============================================================================
# Função de densidade de probabilidade f(x)
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,0),
    mar=c(4.5,5,2,1),bty="n")
plot(xseq,cst*yseq,type="l",xlab=expression(x),ylab=expression(f(x)),lwd=3,
     col="darkblue",xlim=c(-4,4),ylim=c(0,0.5))

#==============================================================================
# Objetivo e' gerar de f(x) = cst*exp( -(x^4)/2 )
# Gerando de uma proposta N(0,1)
# C e' o valor que maximiza a aceitacao para N(0,1)
#==============================================================================
C <- cst*sqrt(2*pi)*exp(1/8)*1.05 # +5% de acrescimo
x  <- rep(NA,n)
k  <- 0
v  <- 0
while( k < n ){
	v    <- v+1
	w    <- rnorm(1)
	prob <- cst*exp(-0.5*(w^4))/( C*dnorm(w) )
	if( runif(1) < prob ){
		k    <- k+1
		x[k] <- w
	}
}

#==============================================================================
# Taxa de aceitacao empirica
#==============================================================================
print(n/v)

[1] 0,7222302

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

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