29 - Simulación, densidad de equilibrio de una distribución gamma MEJORADO

Simulaciones de una densidad de equilibrio gamma 2.

David J. Santana

16 de abril de 2017

Una simulación mejorada, gracias a una mejor elección de la

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)