# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo para gerar de uma normal padrao pelo metodo da rejeicao utilizando
# propostas U(-10,10) e Cauchy
#==============================================================================
set.seed(54345)  # Semente
n      <- 100000 # Tamanho da amostra

#==============================================================================
# Gerando de uma proposta U(-10,10)
# C e' o valor que maximiza a aceitacao para U(-10,10)
#==============================================================================
C  <- 20/sqrt(2*pi) 
x  <- rep(NA,n)
k  <- 0
v  <- 0
while( k < n ){
	v    <- v+1
	w    <- runif(1,-10,10)
	prob <- dnorm(w)/( C*dunif(w,-10,10) )
	if( runif(1) < prob ){
		k    <- k+1
		x[k] <- w
	}
}

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

[1] 0,1253673

#==============================================================================
# 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)

#==============================================================================
# Gerando de uma proposta Cauchy
# D e' o valor que maximiza a aceitacao para Cauchy
#==============================================================================
D  <- dnorm(1)/dt(1,1)
z  <- rep(NA,n)
k  <- 0
u  <- 0
while( k < n ){
	u    <- u+1
	w    <- rt(1,1)
	prob <- dnorm(w)/( D*dt(w,1) )
	if( runif(1) < prob ){
		k    <- k+1
		z[k] <- w
	}
}

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

[1] 0,6592479

#==============================================================================
# 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(z,nclass=50,prob=T,main="",ylim=c(0,0.40),xlim=c(-5,5),col="darkblue",
     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
#==============================================================================