B Código de S-Plus para el Ejemplo B

# ANALISIS DEL MODELO CAUCHY VIA GIBBS SAMPLING


gibbscau_function(N, R, l10, l20)
{
# N: Tamanio de la muestra
# R: Numero de iteraciones
# l10,l20: Valores iniciales
        sal <- list(muestras=0,medias=0)
        l1 <- rgamma(N, l10)
        l2 <- rgamma(N, l20)
        thm <- 1:R
        l1m <- 1:R
        l2m <- 1:R
        for(i in 1:R) {
                mu <- (-5 * l1 + 3 * l2)/(l1 + l2)
                sc <- 1/(l1 + l2)
                th <- rnorm(N, mu, sqrt(sc))
                thm[i] <- mean(th)
                a <- 1
                b1 <- 0.5 * (1 + (th + 5)^2)
                b2 <- 0.5 * (1 + (th - 3)^2)
                l1 <- rgamma(N, a)/b1
                l2 <- rgamma(N, a)/b2
                l1m[i] <- mean(l1)
                l2m[i] <- mean(l2)
                print(i)
        }
        sal$muestras <- cbind(th,l1,l2)
        sal$medias <- cbind(thm,l1m,l2m)
        sal
}


N_5000
R_1100

gc_gibbscau(N,R,0.5627507,0.5363435)


#------------
# Diagnostico

th_gc$muestras[,1]
thm_gc$medias[,1]

l1_gc$muestras[,2]
l2_gc$muestras[,3]
l1m_gc$medias[,2]
l2m_gc$medias[,3]

 
ITC_trunc(R/11)
IT_R
thme_thm[(ITC+1):IT]
l1me_l1m[(ITC+1):IT]
l2me_l2m[(ITC+1):IT]
thmeg_cumsum(thme)/(1:(IT-ITC))
l1meg_cumsum(l1me)/(1:(IT-ITC))
l2meg_cumsum(l2me)/(1:(IT-ITC))

win.graph()
par(mfrow=c(3,2))

plot(thmeg,type="l",xlab="Iteracion",ylab="theta")
title("Promedios ergodicos: theta")
hist(th,xlab="theta",main="Distribucion final de theta",nclass=100,prob=T)

plot(l1meg,type="l",xlab="Iteracion",ylab="lambda_1")
title("Promedios ergodicos: lambda_1")
hist(l1,xlab="l1",main="Distribucion final de lambda_1",nclass=100,prob=T)

plot(l2meg,type="l",xlab="Iteracion",ylab="lambda_2")
title("Promedios ergodicos: lambda_2")
hist(l2,xlab="l2",main="Distribucion final de lambda_2",nclass=100,prob=T)


#---------------------------
# Medias y Varianzas finales

c(mean(th),var(th))
>
> [1]  -1.013873  17.321513
>

c(mean(l1),var(l1))
>
> [1] 0.5581948 1.3310835
>

c(mean(l2),var(l2))
>
> [1] 0.5502849 1.2409282
>


#------------------------
# Densidad final de theta

final_function(tha, l1, l2)
{
        N <- length(l1)
        promedio <- seq(0, 0, length = N)
        suma <- 0
        media <- 0
        var <- 0
        for(i in 1:N) {
                media <- (l1[i] * (-5) + l2[i] * 3)/(l1[i] + l2[i])
                var <- 1/(l1[i] + l2[i])
                suma <- suma + dnorm(tha, media, sqrt(var))
                promedio <- suma/N
        }
        promedio
}


win.graph()
a_-15
b_15

tha_seq(a,b,length=500)
denth_final(tha,l1,l2)


plot(tha,denth,xlab="theta",ylab="p(theta|x)",xlim=c(-15,15),type="l")
title("Densidad final de theta")


#--------------------------------------------------
# Encuentra las modas de la densidad final de theta

locator(2)
>
> [1] -4.876482  2.847174
>


win.graph()
hist(th,xlab="theta",nclass=100,xlim=c(-15,15),ylim=c(0,0.2),prob=T,
main="Distribucion final de theta")

lines(tha,denth)