In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo para gerar de uma normal padrao pelo metodo do envelope utilizando
# proposta U(-10,10)
#==============================================================================
set.seed(54345) # Semente
n <- 100000 # Tamanho da amostra
In [2]:
#==============================================================================
# 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)
k1 <- 0 # Contagem de valores aceito na Etapa 1
k2 <- 0 # Contagem de valores aceito na Etapa 2
v <- 0 # Numero total de iteracoes
m <- 0 # Numero de avaliacoes da Etapa 2
while( (k1+k2) < n ){
v <- v+1
w <- runif(1,-10,10)
u <- runif(1)
if( u < (1-w^2/2) ){
k1 <- k1+1
x[k1+k2] <- w
}else{
m <- m+1
# Esse passo e' utilizado um numero menor de vezes
prob <- dnorm(w)/( C*dunif(w,-10,10) )
if ( u < prob){
k2 <- k2+1
x[k1+k2] <- w
}
}
}
In [3]:
#==============================================================================
# Numero total de iteracoes
# Mesmo numero de avaliacoes de 1-x^2/2
#==============================================================================
v
797656
In [4]:
#==============================================================================
# Taxa de aceitacao empirica da primeira etapa
#==============================================================================
print(k1/v)
[1] 0,09422232
In [5]:
#==============================================================================
# Taxa de aceitacao empirica da segunda etapa
#==============================================================================
print(k2/v)
[1] 0,031145
In [6]:
#==============================================================================
# Taxa de aceitacao empirica
#==============================================================================
print(n/v)
print((k1+k2)/v)
[1] 0,1253673 [1] 0,1253673
In [7]:
#==============================================================================
# Numero de avaliacoes de dnorm(w)/(C*dunif(w,-10,10))
#==============================================================================
print(m)
print(v)
print(m/v) # Taxa de avaliacao de dnorm(w)/(C*dunif(w,-10,10))
[1] 722499 [1] 797656 [1] 0,9057777
In [8]:
#==============================================================================
# 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)
In [9]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================