# Simple BNPR/MCMC Vignette In phylodyn: Statistical Tools for Phylodynamics

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.