Page 88 - Tequio 11

86 Métodos de reconstrucción filogenética II/Duchen/81-89

##-------Metropolis-Hastings para inferencia bayesiana-------## cat("\nEjemplo Metropolis-Hastings para
inferencia Bayesiana\n")
#Se comienza con la topología inicial (paso 1).
indice_actual <- 1
T_actual <- T[ [ indice_actual ] ]

for (i in 2:length(T) ) {
#La topología candidata es la siguiente en la lista de T (paso 2). T_candidata <- T[[i]]

#Se calcula la verosimilitud de la topología actual.
logV_actual <- pml(T_actual,D)

#Se calcula la verosimilitud de la topología candidata.
logV_candidata <- pml(T_candidata,D)

#Se calcula la relacion de verosimilitudes A (paso 3).
A <- logV_actual$logLik/logV_candidata$logLik

if (A>1) { #Si A>1 (paso 4 parte 1).
T_actual <- T_candidata
indice_actual <- i
print(paste("T candidata =",i,", log V = ",logV_candidata$logLik))
} else {
#Si A<1 (paso 4 parte 2).
nuevoIndice <- sample(c(indice_actual,i),1,prob=c(A,1-A)) T_actual <- T[[nuevoIndice]]

if (nuevoIndice==indice_actual) {
print(paste("T actual =",indice_actual,", log V = ",logV_actual$logLik))
} else {
print(paste("T candidata =",indice_actual,", log V = ",logV_candidata$logLik))
}
indice_actual <- nuevoIndice
}
}

#######-FIN DEL PROGRAMA-#######

El output de este programa es el siguiente:

Ejemplo Metropolis-Hastings para inferencia Bayesiana
[1] "T actual = 1 , log Verosimilitud = -87.1232112880638"
[1] "T candidata = 3 , log Verosimilitud = -87.0763445725479"

Tequio, enero-abril 2021, vol. 4, no. 11