# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo para gerar de uma normal padrao pelo metodo da rejeicao utilizando
# proposta Laplace(1)
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library("smoothmest"))
#==============================================================================
set.seed(54345) # Semente
n  <- 100000    # Tamanho da amostra

#==============================================================================
# Gerando de uma proposta ExpDupla(1)=Laplace(1)
# C e' o valor que maximiza a aceitacao para Laplace(1)
#==============================================================================
C  <- dnorm(1)/ddoublex(1,0,1)
x  <- rep(NA,n)
k  <- 0
v  <- 0
while( k < n ){
	v    <- v+1
	w    <- rdoublex(1,0,1)
	prob <- dnorm(w)/(C*ddoublex(w,0,1) )
	if( runif(1) < prob ){
		k    <- k+1
		x[k] <- w
	}
}

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

[1] 0,7587599

#==============================================================================
# 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=50,prob=T,main="",ylim=c(0,0.40),xlim=c(-5,5),col="darkgreen",
     ylab=expression(f(x)),xlab=expression(x))
xseq <- seq(-5,5,length=1000)
yseq <- dnorm(xseq)
lines(xseq,yseq,col="red",lwd=3)

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