Probabilistic programming for birth-death models of evolution using analive particle filter with delayed sampling

Jan Kudlicka
Uppsala University
Uppsala, Sweden

Lawrence M. Murray
Uber AI LabsSan Francisco, CA, USA

Abstract
We consider probabilistic programming forbirth-death models of evolution and introduce
a new widely-applicable inference methodthat combines an extension of the alive
particle filter (APF) with automatic RaoBlackwellization via delayed sampling. Birthdeath models of evolution are an importantfamily of phylogenetic models of the diversification processes that lead to evolutionary trees.
Probabilistic programming languages (PPLs)
give phylogeneticists a new and exciting tool:
their models can be implemented as probabilistic programs with just a basic knowledge ofprogramming. The general inference methods
in PPLs reduce the need for external experts,
allow quick prototyping and testing, and accelerate the development and deployment of
new models. We show how these birth-deathmodels can be implemented as simple programs in existing PPLs, and demonstrate theusefulness of the proposed inference method
for such models. For the popular BiSSE model
the method yields an increase of the effective
sample size and the conditional acceptance rateby a factor of 30 in comparison with a standard bootstrap particle filter. Although concentrating on phylogenetics, the extended APF
is a general inference method that shows its
strength in situations where particles are often assigned zero weight. In the case when theweights are always positive, the extra cost ofusing the APF rather than the bootstrap particle filter is negligible, making our method asuitable drop-in replacement for the bootstrap
particle filter in probabilistic programming inference.

Fredrik Ronquist
Swedish Museumof Natural History
Stockholm, Sweden

1

Thomas B. Schön
Uppsala University
Uppsala, Sweden

INTRODUCTION

The development of new probabilistic models of evolution is an important part of statistical phylogenetics.
These models require inference algorithms that are able
to cope with increased model complexity as well as thelarger amount of observational data available today. Experts from several fields typically need to be involved,
both to design bespoke inference algorithms, and to implement the new models and the inference algorithms
in existing software or to develop new software fromscratch. Probabilistic programming languages (PPLs)(e.g., Goodman et al., 2008; Tolpin et al., 2016; Mansinghka et al., 2014; Paige and Wood, 2014) have the potential to accelerate this: generative models are specified
as simple programs and compiled into executable applications that include general inference engines. Writing
models in PPLs requires just basic programming skills,
and thus allows quick prototyping and testing.
Quite a few software applications for statistical phylogenetics exist today, including the popular MrBayes(Huelsenbeck and Ronquist, 2001) and BEAST (Drummond and Rambaut, 2007). They typically take a
Bayesian approach and implement Markov chain Monte
Carlo inference (see review by Nascimento et al., 2017).
Most of these applications do not allow the user to specify models outside of a predefined model space, which
can be quite narrow. Even when adding new models is
possible, it is usually a challenging task requiring not
only good programming skills but also detailed knowledge of the software design and implementation.
Statistical phylogeneticists recognize the benefits of software that supports the addition of new models and inference methods. For example, the design of BEAST 2(Bouckaert et al., 2014) allows users to create and use
custom modules. RevBayes (Höhna et al., 2016) goeseven further: it uses a domain-specific probabilistic programming language for phylogenetics based on probabilistic graphical models (e.g., Koller and Friedman,2009). However, the language is not Turing-complete,
which means it has some limitations. For example it does
not allow unbounded recursion.
In this paper we concentrate on birth-death models of
evolution, an important family of phylogenetic models.
In these models, births correspond to lineage splits (speciation events) and deaths to extinction events. Thesemodels specify probability distributions of evolutionary
trees and the task is to infer model parameters given apart of a complete tree that represents evolution of currently living species.
We take a step toward using PPLs in statistical phylogenetics: our main contribution is a new general inference
algorithm based on an extension of the alive particle filter (APF) (Del Moral et al., 2015) combined with automatic Rao-Blackwellization via delayed sampling (Murray et al., 2018). We also show how to implement birthdeath models of evolution in existing PPLs, and show
the usefulness of our inference algorithm for such models. Interestingly, by using this algorithm we avoid sampling of birth and death rates. We believe that the algorithm may be of interest for other models with highlyinformative observations. Finally, we prove that the estimator of the marginal likelihood in the extended APF isunbiased.
The rest of the paper is organized as follows: in Section 2
we give a brief recapitulation of basic concepts in evolution and introduce probabilistic programming in moredetail. We derive our inference algorithm and show how
phylogenetic birth-death algorithms can be implemented
in PPLs in Section 3. We give implementations of two
well-known phylogenetic birth-death models and compare several general inference algorithms for these models in Section 4. We offer some conclusions and ideas for
future research in Section 5.

2

BACKGROUND

2.1

SPECIATION, EXTINCTION AND
PHYLOGENIES

There are two types of events that play a significant rolein the evolution of any species. Speciation occurs when
the population of one species splits and eventually forms
two new species. Extinction occurs when the whole population of one species dies out. Species that are not extinct, i.e., species with individuals alive at the present
time, are called extant.
In phylogenetics, the before present (BP) time is usuallyused for dating, i.e., if an event happened at time τ itmeans it happened τ time units ago.

The result of an evolutionary process is a binary tree
called the complete phylogeny. A very simple example
of a complete phylogeny is depicted in Figure 1a. The
nodes represent events and species at significant times:• the root node represents the most recent common
ancestor (MRCA) of all species of interest,• an internal node represents a speciation event,• a leaf at τ = 0 (i.e. the present time) represents an
extant species,• a leaf at τ > 0 (i.e. in the past) represents an extinction event.
The length of edges—or branches as they are called inphylogenetics—is the difference between the time of theparent and the child node.
The reconstructed phylogeny is obtained from a complete tree by removing all subtrees that involve only extinct species. We will refer to this as pruning. An example of a reconstructed tree is depicted in Figure 1b.
The reconstructed phylogeny represents the evolution of
the extant species and only contains information that canbe observed directly (the extant species) or reconstructedby statistical analysis of the DNA sequences of extant
species (the topology of the tree and the times of the speciation events).
2.2

PROBABILISTIC PROGRAMMING

The development of new probabilistic models and inference algorithms is a time-consuming and possibly errorprone process that usually requires skilled experts inprobability, statistics and computer science. Probabilistic programming is a relatively new approach to solve
this problem: generative models are expressed as computer programs in probabilistic programming languages(PPLs) with support for random variables and operations
on them. Integral to PPLs are general inference enginesthat perform the inference in such programs. These engines estimate the distribution of all latent random variables conditioned on the observed data and use it to answer the queries of interest.
PPLs allow us to define and initialize random variables
with a given distribution, for example:x ∼ Normal(0, 1)
The program may use random variables as though any
ordinary variable, and control the flow of the execution.
Depending on the PPL, conditioning on the observed
data might be specified explicitly or implicitly. The former means that conditioning on the observed data is apart of the program, for example:

(a) Complete tree
Root
(MRCA)

τMRCA

Time τ

Observed
speciation(b) Reconstructed tree

Unobserved(hidden)speciation

Pruning0

Augmentation

Extinction

Extant species

Hidden (unobserved) subtree

Figure 1: Complete (on the left) vs. reconstructed (on the right) tree. The reconstructed tree shows only the evolutionof the extant species.
x ∼ Normal(0, 1)
observe 0.892 ∼ Normal(x, 1)
The latter implies that observed values of random variables are not part of the program, but instead specified at
run time (e.g. as arguments).
The main general inference methods used in PPLs are
adaptations of various inference algorithms for the universal setting described below, including Markov chain
Monte Carlo (MCMC) (Metropolis et al., 1953; Hastings, 1970), sequential Monte Carlo (SMC) (Gordon
et al., 1993; Del Moral et al., 2006), and Hamiltonian
Monte Carlo (HMC) (Neal, 2011).
There already exist quite a few PPLs today based on different programming paradigms, for example functional
PPLs like Anglican (Tolpin et al., 2016) and Venture(Mansinghka et al., 2014); imperative Probabilistic C(Paige and Wood, 2014), Turing (Ge et al., 2018), Stan
(Carpenter et al., 2017), Edward (Tran et al., 2016) andPyro (Bingham et al., 2019); and object-oriented Birch(Murray and Schön, 2018).
2.3

PROGRAMMATIC MODEL AND SMC

The execution of a probabilistic program can be modeledusing a programmatic model (Murray and Schön, 2018).
Let {Vi }i denote a countable set of all random variables,
both latent and observed, in a probabilistic program. Thisset might be infinite due to loops and recursion. During execution, whenever a random variable Vi is encountered, its realization vi is drawn from a distribution associated with it.
Multiple executions of the program might in general encounter different subsets of the random variables (e.g.,

due to using random variables in conditional expressions) and encounter them in a different order (with the
exception of the first one). For each random variable Vj
not encountered during the execution we set vj = ⊥ (thesymbol ⊥ represents an undefined value). We will however assume that any execution encounters all observed
random variables and that these observations are encountered in the same order.
Let σ denote a sequence of indices into {Vi } specifying
the order in which the random variables are encounteredduring an execution of the program, and let |σ| denotethe length of this sequence. Also, let {vi }i∈σ denote therealizations of the random variables indexed by σ, i.e.,vσ[1] , . . . , vσ[|σ|] . In a similar manner, we will use {Vi =vi }i∈σ to denote Vσ[1] = vσ[1] , . . . , Vσ[|σ|] = vσ[|σ|] .
The index of the k-th encountered variable is given by adeterministic function Ne (for next) of the realizations of
the previously encountered random variables, so that
σ[k] = Ne({vi }i∈σ[1:k−1] ),where σ[1:k−1] denotes the indices of the first k − 1 encountered random variables. The function Ne is uniquely
defined by the probabilistic program. If there are no more
random variables to encounter, Ne returns ⊥.
The k-th encountered random variable, Vσ[k] , is sampled
fromVσ[k] ∼ pσ[k] (·| Pa({vi }i∈σ[1:k−1] )),where pσ[k] is the distribution specified by the program,
Pa (for parents) is a deterministic function returning theparameters of this distribution, and again, it is a functionof the realizations of the previously encountered randomvariables.

Algorithm 1 Bootstrap particle filter (BPF).
for n = 1 to N do
. Initialize(n)(n)
v0 ← ∅; w0 ← 1/N
for t = 1 to T do
for n = 1 to N do P
(m)
(l)
N
a ∼ Cat({wt−1 / l=1 wt−1 }Nm=1 ) . Resample(n)(a)
vt ← P ROPAGATE(vt−1 ). Propagate(n)(n)
wt ← pγ[t] (yt | Pa(vt )). Weigha test function h of interest:
P (n)  (n) w h vTb = n PT
E[h]
,
(n)
n wT
as well as to estimate the marginal likelihood p(y1:T ):
Zb =

(n)

N
T
Y
1 X

N n=1
t=1

(n)

wt .

vt [γ[t]] ← yt
function P ROPAGATE(v)κ ← Ne(v)while κ 6∈ γ dov[κ] ∼ pκ (·| Pa(v))κ ← Ne(v)return v

3

METHODS

3.1

EXTENDED ALIVE PARTICLE FILTER. Run until next observe

The joint distribution function encoded by the program
can be given recursively (starting with σ = []):
p({vi }i6∈σ |σ, {Vj = vj }j∈σ )0 p({vi }i6∈σ0 |σ , {Vj = vj }j∈σ0 )×pκ (vκ | Pa({vj }j∈σ ))if κ 6= ⊥,=
1if κ = ⊥ ∧ ∀i 6∈ σ : vi = ⊥,





0
otherwise,
where κ = Ne({vi }i∈σ ) and σ 0 is obtained from σ byappending κ. The first case is the conditional probability
chain rule, the remaining cases cover the situation where
there are no more random variables to encounter.

In the bootstrap particle filter, each particle is propagatedby simulating the prior, and may make random choicesthat lead to a state with zero weight. In phylogenetic
birth-death models this happens quite often: when simulating the evolution of subtrees that must ultimately become extinct, if any species happen to survive to thepresent time, the particle is assigned zero weight. In extreme cases, all particles have zero weight, and the BPF
degenerates.
Del Moral et al. (2015) considered this problem in a setting with indicator potentials (such as in approximate
Bayesian computation), i.e. all weights being either zeroor one. They proposed a modification of the BPF, where
the resampling and propagation steps are repeated forparticles that have weight zero until all particles haveweight one. Details of the resulting alive particle filter(APF) as well as proofs of some of its theoretical properties can be found in Del Moral et al. (2015).
Although the original APF was designed specifically forindicator potentials, we have extended the algorithm towork with importance weights, see Algorithm 2 (theP ROPAGATE function is the same as in Algorithm 1). The
APF requires N + 1 particles rather than N in order toestimate the marginal likelihood without bias.

We wish to sample from the posterior distributionp({vi }i6∈γ |Vγ[1] = y1 , . . . , Vγ[T ] = yT ), where T denotes the number of observations, yt denotes the t-th observation and γ denotes the sequence of indices of the observed random variables in {Vi }. The sequential natureof the joint distribution allows us to employ Sequential
Monte Carlo methods (Del Moral et al., 2006) to sample
from this posterior distribution, including the bootstrap
particle filter (BPF) summarized in Algorithm 1. For thesake of brevity we have assumed that the last observationis also the last encountered random variable. In the pseudocode, Cat() denotes the categorical distribution withthe given event probabilities. Variables denoted by v areassociative arrays (also known as maps or dictionaries)
used to store the realizations of random variables (v[i]denotes the realization of Vi ). The P ROPAGATE function
runs the program until it encounters an observation.

At the t-th observe statement, if the weight of a particle iszero, the resampling and propagation steps are repeated.
This procedure is repeated until the weights of all N + 1
particles are positive. The APF counts the total numberof propagations Pt for each observation. The algorithm
never uses the states or weights of the N + 1-th particle,
but propagations made using this particle are includedin Pt , and used to calculate the unbiased estimate of themarginal likelihood p(y1:T ):
T PN
(n)
Y
n=1 wtb.
Z=
Pt − 1t=1

Samples from the joint distribution and the corresponding weights can be used to estimate the expected value of

The proof of unbiasedness can be found in Appendix A
in the supplementary material.

Algorithm 2 Alive particle filter (APF).
for n = 1 to N do
. Initialize(n)(n)
v0 ← ∅; w0 ← 1/N
for t = 1 to T do
Pt ← 0for n = 1 to N + 1 do
repeat. Resample(m) PN(l) Na ∼ Cat({wt−1 / l=1 wt−1 }m=1 )(n)(a)
vt ← P ROPAGATE(vt−1 ). Propagate
Pt ← Pt + 1(n)(n)
wt ← pγ[t] (yt | Pa(vt )). Weighwhere ∆r is the length of the branch between the noder and its parent. The factor 2 occurs in the formulas
because the orientation of the two child subtrees does
not matter. If r is a speciation event, no extinction occurs along the branch (factor e−µ∆r ) and the speciation
happens after a waiting time ∆r (factor λ exp−λ∆r ). If
r is an extinction event, no speciation occurs along thebranch (factor e−λ∆r ) and the extinction occurs after
waiting time ∆r (factor µe−µ∆r ). Finally, if r is anextant species, neither extinction nor speciation occurs
along the branch. The likelihood of the complete phylogeny T is given by p(Troot ).

(n)until wt > 0(n)vt [γ[t]] ← yt

Unbiasedness of the marginal likelihood estimate opens
for the possibility to use the APF within particle Markov
chain Monte Carlo methods.
3.2

BIRTH-DEATH MODELS AS
PROBABILISTIC PROGRAMS

Phylogentic birth-death models constitute a family ofmodels where speciation (birth) events and extinction(death) events occur along the branches of a phylogenetic tree. Typically, the waiting times between eventsare exponentially distributed. In general, the rates ofthese exponential distributions do not remain constant
but rather change continuously, discontinuously, or both.
Some models assume that these rates further depend ona state variable that itself evolves discontinuously along
the tree; in some cases the value of this state variable isgiven for the extant species.
The constant-rate birth death (CRBD) model (Kendall,1948) is the simplest birth-death model, where the speciation rate λ and extinction rate µ remain constant overtime. Pseudocode for generating phylogenetic trees using the CRBD model can be found in Appendix B in the
supplementary material.
Let Tr denote the subtree rooted at the node r, andCh(r) denote the children of this node. The likelihoodof the subtree Tr can be expressed recursively (we have
dropped conditioning on λ and µ in the notation forbrevity): Q
2p(Tc )
if r is the root node,c∈Ch(r)Qp(Tc ) if r is a speciation,
 2λe−(λ+µ)∆rc∈Ch(r)p(Tr ) =µe−(λ+µ)∆r
if r is an extinction,






 −(λ+µ)∆re
if r is an extant species,

The likelihood in other birth-death models that admit
varying rates and/or include the state can be derived ina similar way.
The likelihood of the complete phylogeny for the CRBD
model can be conveniently written asp(T ) = 2S+1 λS µX e−(λ+µ)L ,where L is the sum of all branch lengths, S the number of
speciation events (excluding the root) and X the numberof extinction events.
The task of interest is to infer model parameters givena reconstructed tree. Recall that this tree is a part ofthe complete tree corresponding to the extant species andtheir ancestors. A naive approach to inference is to simulate trees from the generative model, prune back the extinct subtrees, and reject those for which the pruned tree
does not equal the observed tree. Such an approach always results in rejection.
Instead, we turn the problem upside down: starting withthe observed tree and augmenting it with unobserved information to obtain a complete tree. Recalling Figure 1,
the observed tree is traversed in depth-first order. At theroot node the importance weight is doubled (since theorientation of the children does not matter). Along eachbranch, the generative model is used to simulate:• changes to the state (in models with state),• changes to the speciation and extinction rates,• hidden speciation events.
For each of the hidden speciation events, the model simulates the evolution of the new species (i.e. a hidden subtree). If any portion of a hidden subtree survives to thepresent time, the weight is set to zero. If not, the current
weight is doubled, since the children created by a hidden
speciation event on an observed branch are interchangeable: either can correspond to the observed subtree, theother must correspond to the hidden subtree and goes extinct.

If the examined branch ends with a speciation event,the algorithm observes 0 ∼ Exponential(λ) and doubles the weight. Finally, as there were no extinction
events along the processed branch, the algorithm observes 0 ∼ Poisson(µ∆). In models with state, if the
branch ends at τ = 0 (i.e. the present time) we also condition on the simulated state being equal to the observed
state. We trigger resampling at the end of each branch.
Let us return to the CRBD model in light of the discussion above. The likelihood of a proposed complete tree
T 0 is given by000

the parameters and sample them after running the particle filter. Exploiting the conjugacy in a probabilistic
program can be automated by an algorithm known as delayed sampling (Murray et al., 2018).
Consider the following prior:ν ∼ Gamma(k, θ),with k ∈ N. When the program needs to make a draw
from a Poisson distributionn ∼ Poisson(ν∆),

0

q(T 0 ) = λH e−λLobs × (2λ)S µX e−(λ+µ)L ,where H 0 denotes the number of all simulated hidden
speciation events along the observed tree, Lobs the sumof the branch lengths in the observed tree, S 0 the numberof speciation events in all hidden subtrees, X 0 the number of extinction events and finally L0 denotes the sumof the branch lengths in the hidden subtrees. The factor0
λH e−λLobs is related to the hidden speciation events, andthe rest is the combined likelihood of all hidden subtrees.
The weight of the proposal T 0 that is compatible with theobservation (i.e. without any extant species in the hidden
subtrees) is given by

it can instead make a draw from the marginalized distribution:1,n ∼ NegativeBinomial k,
1 + ∆θ
where NegativeBinomial(k, p) is the negative binomial
distribution counting the number of failures given thenumber of successes k and the probability of success pin each trial. The distribution for ν is then updated to theposterior distribution according toθν ∼ Gamma k + n,.
1 + ∆θ0w(T 0 ) = 2H × 2Sobs +1 × λSobs × e−µLobs ,where Sobs is the number of observed speciation events.
0
The factor 2H corresponds to doubling the weight foreach hidden subtree, the factor 2Sobs +1 corresponds todoubling the weight at the root node and at all observed
speciations, the factor λSobs is due to observing 0 ∼Exponential(λ) at all observed speciations, and finallythe factor e−µLobs is due to observing 0 ∼ Poisson(µ∆)
for all observed branches.
Multiplying the likelihood q(T 0 ) and the weight w(T 0 )of the proposal and summing the event numbers and thebranch lengths we getq(T 0 )w(T 0 ) = p(T 0 ).
3.3

DELAYED SAMPLING OF THE RATESIn a Bayesian setting, the parameters are associated witha prior distribution. Using the gamma distribution (orthe exponential distribution as its special case) as a priorfor the rates of speciation, extinction and state change ismathematically convenient since the gamma distribution
is a conjugate prior for both the Poisson and the exponential likelihood. Instead of sampling these parametersfrom the prior distribution before running the particle filter (which we refer to as immediate sampling), we can
exploit the conjugacy, which allows us to marginalize out

Similarly for variables distributed according to the exponential distribution, instead of drawing∆ ∼ Exponential(ν),
the program makes a draw from the marginalized distribution:


1
,k ,
∆ ∼ Lomax
θwhere the first parameter denotes the scale and the second the shape of the Lomax distribution, and the distribution for ν is then updated toθν ∼ Gamma k + 1,.
1 + ∆θ
Using this strategy there is actually no need to sample
the rates at all; all draws involving these rates are replaced by draws from the negative binomial and the Lomax distributions with a consequent update of the rate
distribution. Details of these conjugacy relationships canbe found in Appendix C in the supplementary material.

4

EXPERIMENTS

We implemented the inference algorithms describedabove in the probabilistic programming language

Birch (Murray and Schön, 2018) and added support forthe conjugacy relationships described in the previous
section, so that Birch can provide automated delayedsampling for these. We also implemented two phylogenetic birth-death models, described in Sections 4.1 and4.2 below.
We ran the inference for these models using different combinations of the inference method, the samplingstrategy (immediate or delayed) and different number ofparticles N . For each combination we executed the program M times, collected the estimates {Zbm }M
m=1 of themarginal likelihood and calculated the relative effective
sample size (RESS):
P
2
M
bm
Z
m=1
1
RESS =
PM b2 ,
MZ
m=1

m

as well as the conditional acceptance ratio (CAR) (seeMurray et al., 2013 for more detail):
!
M
X
1
CAR =
2
ci − 1 ,M
i=1bm }m .
where ci is the sum of the i smallest elements in {Zb
We also calculated the sample variance var log Z.
For the experiments with the APF we also compared thetotal number of propagations with the number of propagations in the BPF by calculatingPMPmρ = m=1,
MNT
where Pm is the number of all propagations made during
the m-th execution and T is the number of branches in
the observation. Note that N T is the number of propagations in the BPF.
4.1

CONSTANT-RATE BIRTH DEATH MODEL

Pseudocode for the probabilistic program implementing
the constant-rate birth death (CRBD) model is listed as
Algorithm 3. To sample speciation events along a branch
we first sample a number of events from a Poisson distribution and then sample the time of each event from a
uniform distribution. The implementation in Birch canbe found in the supplementary material.
We used the phylogeny of cetaceans (Steeman et al.,2009) as the observation. This phylogeny (Figure 2 inthe supplementary material) represents the evolution ofwhales, dolphins and porpoises and contains 87 extant
species. We used Gamma(1, 1) as the prior for boththe speciation and extinction rate. The results of experiments comparing BPF and APF with immediate or delayed sampling for different number of particles N , andAlgorithm 3 CRBD model as a probabilistic program.
Input:• T – a pre-ordered list of nodes in the observation• kλ , θλ – the shape and scale of the prior Gamma
distribution for λ (kλ ∈ N)
• kµ , θµ – the shape and scale of the prior Gamma
distribution for µ (kµ ∈ N)
λ ∼ Gamma(kλ , θλ )µ ∼ Gamma(kµ , θµ )for all r ∈ T do
if r is the root thendouble the weightcontinuechs ∼ Poisson(λ∆r )for i ← 1 to chs doτ ∼ Uniform(τr , τr + ∆r )if B RANCH S URVIVES(τ ) then
set the weight to 0 and returnelse
double the weight
if r has children thenobserve 0 ∼ Exponential(λ)
double the weightobserve 0 ∼ Poisson(µ∆r )function B RANCH S URVIVES(τ, λ, µ)∆ ∼ Exponential(µ)
if ∆ ≥ τ thenreturn truecb ∼ Poisson (λ∆)
for i ← 1 to cb doτ 0 ∼ Uniform(τ − ∆, τ )if B RANCH S URVIVES(τ 0 , λ, µ) thenreturn true
return false
running M = 200 executions for each combination, aresummarized in Table 1 and Figure 3a.
When using delayed sampling, the speciation and extinction rates are never sampled; the rates are instead represented by gamma distributions with parameters that are(m) (m) (m)updated during the execution. Let kλ , θλ , kµ and(m)θµ denote the final values of these parameters for aparticle drawn from all particles in the m-th run withthe probabilities proportional to their final weights. The
posterior distributions for λ and µ can be estimated by
mixtures of gamma distributions:λ ∼ PM

M
X

1

m=1 Zm m=1b(m) (m)Zbm Gamma kλ , θλ,Our implementation of the BiSSE model can be found in
the supplementary material.

λ
µ

40

20

5

0
0

0.05

0.1

0.15

0.2

Figure 2: The posterior distributions for the speciationand extinction rates for the cetaceans using the CRBD
model.

µ ∼ PM

M
X

1

m=1 Zm m=1


Zbm Gamma kµ(m) , θµ(m) .

b

Figure 2 depicts the posterior distributions for the speciation and extinction rates estimated using M = 1000
runs of the APF with N = 4096 particles.
4.2

BINARY-STATE SPECIATION AND
EXTINCTION MODEL

The binary-state speciation and extinction (BiSSE)
model (Maddison et al., 2007) introduces a binary state
for species, denoted by s ∈ {0, 1}. Each state has its own(but constant) speciation and extinction rates, denoted by
λs and µs . The waiting time between switching stateis exponentially distributed with rates q01 for switchingfrom state 0 to state 1, and q10 from state 1 to state 0. In
our experiments we made a (common) assumption thatq01 = q10 = ς.
We used the same observation as Rabosky and Goldberg(2015), i.e., we extended the cetacean phylogeny with thestate variable related to the body length of cetaceans obtained from Slater et al. (2010). Data are available for 74of the 87 extant species. The binary state 0 and 1 refers
to the length being below and above the median. Again
we used Gamma(1, 1) as the prior for λ0 , λ1 , µ0 and µ1 ,and Gamma(1, 10/820.28) as the prior for ς (the number in the denominator is the sum of all branch lengths).
The initial state at the root is drawn from {0, 1} withequal probabilities. The results for experiments comparing the BPF and the APF with immediate or delayedsampling for different number of particles N , and running M = 200 executions for each combination, aresummarized in Table 2 and Figure 3b. When runningthe experiments using the BPF and immediate sampling,a certain fraction of the executions degenerated—from
25% of the executions with 1024 particles down to 1.5%of the executions with 4096 particles. These executionsbwere excluded when calculating var log Z.

DISCUSSION AND CONCLUSION

In this paper we introduced a new general inference
method for probabilistic programming combining an extended alive particle filter (APF) with delayed sampling,
and proved that the resulting estimate of the marginal
likelihood is unbiased. We showed how phylogenetic
birth-death models can be implemented in probabilistic programming languages, in particular, we considered
two models—CRBD and BiSSE and their implementation in the probabilistic programming language Birch.
We showed the strength of this inference method forthese models compared to the standard bootstrap particle filter (BPF) (Tables 1 and 2, and Figures 3a and 3b):for the BiSSE model using 8192 particles we increased
RESS approximately 29 times, CAR approximately 30
times and lowered var log Zb more than 1150 times at thecost of running 3 times more propagations.
The extended APF is a suitable drop-in replacement for
the BPF for black-box probabilistic programs. If a program produces only positive weights, the APF produces
the same result as the BPF at the overhead of just one extra particle, used to estimate the marginal likelihood. On
the other hand, if the program can produce zero weights,
the APF behaves much more reasonably than the BPF:
resampling and propagation are repeated until all particles have positive weight. This may seem equivalent tousing the BPF with a higher number of particles (ρ times
more to be precise), but this is not the case: the numberof propagations Pt is not the same throughout the execution, but rather adapted dynamically for each t. This
simplifies the tuning of the number of particles for suchmodels.
The learning of rates in birth-death models sits in thecontext of broader problems in phylogenetics, such as
the learning of trees. Interesting future work would be toconsider whether models and methods for learning rates
can be combined with models and methods for learning
tree structures for joint inference.
Acknowledgements
The authors wish to thank Johannes Borgström, Viktor
Senderov and Andreas Lindholm for their useful feedback. This research was financially supported by the
Swedish Foundation for Strategic Research (SSF) via theproject ASSEMBLE and by the Swedish Research Council grants 2013-4853, 2014-05901 and 2017-03807.

-250

-225
-275

b
log Z

log Zb

-230
-300

-235
-325

MethodMethodBPF immediate

BPF immediate

APF immediate-240

APF immediate-350

BPF delayed

BPF delayed

APF delayed
512

1024

2048

APF delayed

4096

1024

2048

N

4096

8192

N

(a) CRBD(b) BiSSE

Figure 3: Box plot of log Zb for the CRBD (on the left) and BiSSE (on the right) model for different number ofparticles (N ) and methods. The lower and upper hinges correspond to the first and third quartile, whereas the midhingecorresponds to the median. Values outside of the interquartile range are shown as dots. The horizontal dashed line
shows the true value of log Z for the CRBD model.

Table 1: Summary of the results of the experiments with the CRBD model using the cetacean phylogeny as theobservation, priors λ, µ ∼ Gamma(1, 1), and M = 200.
Bootstrap particle filter (BPF)
Immediate sampling
N
512
1024
2048
4096

RESS CAR
0.02 0.04
0.11 0.12
0.14 0.17
0.18 0.23

var
334.2
117.5
52.2
17.2

Alive particle filter (APF)

Delayed sampling
RESS CAR
0.13 0.23
0.28 0.35
0.47 0.47
0.67 0.63

var
31.6
12.9
8.7
0.7

Immediate samplingρ RESS CAR
1.8 0.11 0.15
1.8 0.14 0.18
1.7 0.29 0.32
1.7 0.36 0.42

var
50.7
20.2
7.6
2.7

Delayed samplingρ RESS CAR
1.7 0.40 0.46
1.7 0.54 0.55
1.7 0.73 0.69
1.7 0.84 0.76

var
2.7
0.8
0.3
0.2

Table 2: Summary of the results of the experiments with the BiSSE model using the cetacean phylogeny extended with information about the body length as the observation, priors λ0 , λ1 , µ0 , µ1 ∼ Gamma(1, 1) andς ∼ Gamma(1, 10/820.28), and M = 200.
Bootstrap particle filter (BPF)
Immediate sampling
N
1024
2048
4096
8192

RESS CAR
0.01 0.01
0.01 0.01
0.01 0.01
0.02 0.02

var3382.22954.01894.1
968.4

Alive particle filter (APF)

Delayed sampling
RESS CAR
0.06 0.09
0.09 0.15
0.22 0.27
0.28 0.35

var
72.4
22.2
7.6
6.1

Immediate samplingρ RESS CAR
10.0 0.01 0.01
6.6 0.02 0.02
5.9 0.01 0.01
3.9 0.02 0.03

var2294.91044.5614.3
160.9

Delayed samplingρ RESS CAR
3.1 0.10 0.21
3.1 0.14 0.27
3.1 0.34 0.43
3.0 0.54 0.55

var
4.8
2.9
1.3
0.8

References
E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer,
N. Pradhan, T. Karaletsos, R. Singh, P. Szerlip,
P. Horsfall, and N. D. Goodman. Pyro: Deep universal probabilistic programming. Journal of Machine
Learning Research, 20(28):1–6, 2019.

J. P. Huelsenbeck and F. Ronquist. MRBAYES: Bayesianinference of phylogenetic trees. Bioinformatics, 17(8):754–755, 2001.
D. G. Kendall. On the generalized “birth-and-death” process. The Annals of Mathematical Statistics, 19(1):1–15, 1948.

R. Bouckaert, J. Heled, D. Kühnert, T. Vaughan, C.H. Wu, D. Xie, M. A. Suchard, A. Rambaut, andA. J. Drummond. BEAST 2: A software platform forBayesian evolutionary analysis. PLOS Computational
Biology, 10(4):e1003537, 2014.

D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT Press, 2009.

B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee,
B. Goodrich, M. Betancourt, M. Brubaker, J. Guo,
P. Li, and A. Riddell. Stan: A probabilistic programming language. Journal of Statistical Software, 76(1),2017.

V. Mansinghka, D. Selsam, and Y. Perov.
Venture: a higher-order probabilistic programming platform with programmable inference. arXiv preprint
arXiv:1404.0099, 2014.

P. Del Moral, A. Doucet, and A. Jasra. Sequential Monte
Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3):411–
436, 2006.
P. Del Moral, A. Jasra, A. Lee, C. Yau, and X. Zhang.
The alive particle filter and its use in particle Markov
chain Monte Carlo. Stochastic Analysis and Applications, 33(6):943–974, 2015.
A. J. Drummond and A. Rambaut. BEAST: Bayesianevolutionary analysis by sampling trees. BMC Evolutionary Biology, 7(1):214, 2007.
H. Ge, K. Xu, and Z. Ghahramani. Turing: A language
for flexible probabilistic inference. In Proceedingsof the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84, pages
1682–1690, 2018.
N. D. Goodman, V. K. Mansinghka, D. M. Roy,K. Bonawitz, and J. B. Tenenbaum. Church: a language for generative models. In Proceedings of the24th Conference in Uncertainty in Artificial Intelligence, pages 220–229, 2008.
N. J. Gordon, D. J. Salmond, and A. F. M. Smith.
Novel approach to nonlinear/non-Gaussian Bayesianstate estimation. In IEE Proceedings on Radar and
Signal Processing, volume 140, pages 107–113, 1993.
W. K. Hastings. Monte Carlo sampling methods using
Markov chains and their applications. Biometrika, 57(1):97–109, 1970.
S. Höhna, M. J. Landis, T. A. Heath, B. Boussau,
N. Lartillot, B. R. Moore, J. P. Huelsenbeck, and
F. Ronquist. RevBayes: Bayesian phylogenetic inference using graphical models and an interactive modelspecification language. Systematic Biology, 65(4):726–736, 2016.

W. P. Maddison, P. E. Midford, and S. P. Otto. Estimating a binary character’s effect on speciation and extinction. Systematic Biology, 56(5):701–710, 2007.

N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth,A. H. Teller, and E. Teller. Equation of state calculations by fast computing machines. The Journal ofChemical Physics, 21(6):1087–1092, 1953.
L. M. Murray and T. B. Schön. Automated learning witha probabilistic programming language: Birch. Annual
Reviews in Control, 46:29–43, 2018.
L. M. Murray, E. M. Jones, and J. Parslow. On disturbance state-space models and the particle marginal
metropolis-hastings sampler. SIAM/ASA Journal on
Uncertainty Quantification, 1(1):494–521, 2013.
L. M. Murray, D. Lundén, J. Kudlicka, D. Broman, andT. B. Schön. Delayed sampling and automatic RaoBlackwellization of probabilistic programs. In International Conference on Artificial Intelligence andStatistics, pages 1037–1046, 2018.
F. F. Nascimento, M. dos Reis, and Z. Yang. A biologist’s guide to Bayesian phylogenetic analysis. Nature
Ecology & Evolution, 1(10):1446–1454, 2017.
R. M. Neal. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11):2, 2011.
B. Paige and F. Wood. A compilation target for probabilistic programming languages. In International
Conference on Machine Learning, pages 1935–1943,
2014.
D. L. Rabosky and E. E. Goldberg. Model inadequacy
and mistaken inferences of trait-dependent speciation.
Systematic Biology, 64(2):340–355, 2015.
G. J. Slater, S. A. Price, F. Santini, and M. E. Alfaro.
Diversity versus disparity and the radiation of modern
cetaceans. Proceedings of the Royal Society of London
B: Biological Sciences, 277(1697):3097–3104, 2010.
M. E. Steeman, M. B. Hebsgaard, R. E. Fordyce, S. Y.
Ho, D. L. Rabosky, R. Nielsen, C. Rahbek, H. Glenner, M. V. Sørensen, and E. Willerslev. Radiation ofextant cetaceans driven by restructuring of the oceans.
Systematic Biology, 58(6):573–585, 2009.
D. Tolpin, J.-W. van de Meent, H. Yang, and F. Wood.
Design and implementation of probabilistic programming language Anglican. In Proceedings of the 28th
Symposium on the Implementation and Application of
Functional programming Languages, pages 6:1–6:12.
ACM, 2016.
D. Tran, A. Kucukelbir, A. B. Dieng, M. Rudolph,
D. Liang, and D. M. Blei. Edward: A library forprobabilistic modeling, inference, and criticism. arXiv
preprint arXiv:1610.09787, 2016.