In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo para gerar de uma normal truncada entre 0 e infinito pelo metodo da
# rejeicao utilizando proposta Gama(1,1)=Exp(1)
#==============================================================================
set.seed(54345) # Semente
n <- 100000 # Tamanho da amostra
In [2]:
#==============================================================================
# Gerando de uma proposta Exp(1)
# C e' o valor que maximiza a aceitacao para Exp(1)
#==============================================================================
C <- 2*dnorm(1)/dexp(1,1)
x <- rep(NA,n)
k <- 0
v <- 0
while( k < n ){
v <- v+1
w <- rexp(1,1)
prob <- 2*dnorm(w)/( C*dexp(w,1) )
if( runif(1) < prob ){
k <- k+1
x[k] <- w
}
}
In [3]:
#==============================================================================
# Taxa de aceitacao empirica
#==============================================================================
n / v
0,759526359362302
In [4]:
#==============================================================================
# Histograma
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(6,5,0),
mar=c(4.5,5,2,1),bty="n")
hist(x,nclass=50,prob=T,main="",ylim=c(0,0.8),xlim=c(0,5),col="darkgreen",
ylab=expression(f(x)),xlab=expression(x))
xseq <- seq(0,5,length=1000)
yseq <- 2*dnorm(xseq)
lines(xseq,yseq,col="red",lwd=3)
In [5]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================