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))
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.
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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.