función de densidad proposal, disminuyendo rechazos.
############################################################
## Si X es una v.a. Gamma(3.5,3.5), el problema consiste ##
## en simular su v.a. de equilibrio, la cual tiene como ##
## función de densidad la función f_e(x)=1-F(x). ##
## Se usará como proposal, la densidad de una v.a. ##
## mezcla de erlangs 1/4(erl(1,3.5) + ... + erl(4,3.5) ##
############################################################
library(latex2exp)
m <- 10000
fe <- function(x){1-pgamma(x,3.5,3.5)}
g <- function(x){
1/4*(dgamma(x,1,3.5)+dgamma(x,2,3.5)+dgamma(x,3,3.5)+dgamma(x,4,3.5))
}
C <- 0.01+1/g(0)
curve(fe,0,6,lwd=3)
curve(C*g(x),col=2,add = T,lwd=3)
abline(h=0,v=0)
lines(c(3.5,4.5),c(0.8,0.8),lwd=3)
lines(c(3.5,4.5),c(0.6,0.6),col=2,lwd=3)
text(5,0.8,TeX('$f_e(x)$'),cex=2.5)
text(5,0.6,TeX('$g(x)$'),cex=2.5)
##############################################
curve(fe,10,16)
curve(C*g(x),col=2,add = T)
# Siempre fe(x)<=C*g(x)
##############################################
curve(fe,0,6)
curve(C*g(x),col=2,add = T)
abline(h=0,v=0)
exitos <- 0
fracasos <- 0
X <- numeric(m)
set.seed(123)
while(exitos < m){
U <- runif(1)
Z <- sample(1:4,1)
Y <- rgamma(1,Z,3.5)
if(U<=fe(Y)/(C*g(Y))){
exitos <- exitos + 1
X[exitos] <- Y
points(Y,U*C*g(Y),pch = 16,cex = 1.5,col = 'green')
}else{
fracasos <- fracasos + 1
points(Y,U*C*g(Y),pch = 16,cex = 1.5,col = 'red')
}
}
fracasos
## [1] 1506
hist(X,freq = F, border = 'blue',lwd =3 , add= T)
title('Éxitos y rechazos en la simulación',cex=2)
mean(X)
## [1] 0.6384351
var(X)
## [1] 0.2531555
mu_real <- 9/14
var_real <- 51/196
mu_real
## [1] 0.6428571
var_real
## [1] 0.2602041
requilibrio <- function(m,seed){
exitos <- 0
X <- numeric(m)
set.seed(seed)
while(exitos < m){
U <- runif(1)
Z <- sample(1:4,1)
Y <- rgamma(1,Z,3.5)
if(U<=fe(Y)/(C*g(Y))){
exitos <- exitos + 1
X[exitos] <- Y
}
}
X
}
a <- requilibrio(50000,123)
hist(a,20,freq = F,col = 'yellow',border = 2,
main = 'Histograma vs densidad',ylab = 'fe(x)')
curve(fe,add=T,lwd=3)
