28 - Simulación, densidad de equilibrio de una distribución gamma

Simulaciones de una densidad de equilibrio gamma.

David J. Santana

16 de abril de 2017

############################################################
## 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. Y=|Z|, ##
## donde Z es una v.a. N(0,(1.2)^2).                      ##
############################################################
library(latex2exp)
m <- 10000
sigma <- 1.2
fe <- function(x){1-pgamma(x,3.5,3.5)}
g <- function(x){2*dnorm(x,0,sigma)}
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,7.2,7.5)
curve(C*g(x),col=2,add = T)
xx <- 7.328
abline(v=xx)

1-pgamma(xx,4.5,3.5)-xx*(1-pgamma(xx,3.5,3.5))

## [1] 2.522786e-09

# Aunque C*g(x)<fe(x) para x>7.328, no es 
# importante porque 1-Fe(7.328) es aprox. 0. 
##############################################
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 <- rnorm(1,0,sigma)
        Y <- abs(Z)
        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] 4918

hist(X,freq = F, border = 'blue',lwd =3 , add= T)
title('Éxitos y fracasos en la simulación',cex=2)

mean(X)

## [1] 0.6464991

var(X)

## [1] 0.2647186

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 <- rnorm(1,0,sigma)
                Y <- abs(Z)
                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)