Simple BNPR/MCMC Vignette

This vignette shows a simple and typical use of phylodyn and its MCMC tools. If we increase the number of MCMC samples, we reproduce part of the split HMC results from Lan, et al. 2014.

We start by loading the phylodyn package.

set.seed(8675309)
library(phylodyn)

We need to set the true effective population size trajectory function, and also its reciprocal. For this example, we choose exponential growth (already implemented in phylodyn).

traj = exp_traj

For simplicity's sake, we will use isochronous sampling (taking all samples simultaneously at t=0). We sample 100 individuals at the present time.

samp_times = 0
n_sampled  = 100

We simulate a genealogy based on our sample using the coalescent.

gene = coalsim(samp_times = samp_times, n_sampled = n_sampled, traj = traj, lower_bound = 1/20)

We first use BNPR to calculate approximate marginals.

res_BNPR = BNPR(data = gene, lengthout = 100)

We set the number of samples and burn-in parameters. For expediency in this vignette we set them to be small. To produce more meaningful results, increase both parameters.

nsamp = 500
nburnin = 100

We invoke the mcmc_sampling function with splitHMC to run MCMC on the data.

res_MCMC = mcmc_sampling(data = gene, alg = 'splitHMC', nsamp = nsamp,
                         nburnin = nburnin)

We plot and compare the results from BNPR and splitHMC.

# Plot the results of BNPR
plot_BNPR(BNPR_out = res_BNPR, traj = exp_traj, col = "blue", traj_col = "black")
title("Exponential growth")

# plot the results of splitHMC
lines(res_MCMC$med_fun, pch="", col='red', lwd=1)
lines(res_MCMC$low_fun, pch="", col='red', lwd=1, lty=1)
lines(res_MCMC$hi_fun,  pch="", col='red', lwd=1, lty=1)
legend('topleft',c('Truth','BNPR',"splitHMC"),
       col=c('black','blue','red'),lwd=c(1,2.5,2.5),bty='n', lty=c(2,1,1))

References

  1. J. A. Palacios and V. N. Minin. Integrated nested Laplace approximation for Bayesian nonparametric phylodynamics. In Proceedings of the Twenty-Eighth International Conference on Uncertainty in Artificial Intelligence, pages 726–735, 2012.

  2. Lan, S., Palacios, J. A., Karcher, M., Minin, V. N., & Shahbaba, B. (2014). An Efficient Bayesian Inference Framework for Coalescent-Based Nonparametric Phylodynamics. arXiv preprint arXiv:1412.0158.



Try the phylodyn package in your browser

Any scripts or data that you put into this service are public.

phylodyn documentation built on May 29, 2017, 1:28 p.m.