Zig-zag sampling for discrete structures and non-reversible phylogenetic MCMC

We construct a zig-zag process targeting a posterior distribution defined on a hybrid state space consisting of both discrete and continuous variables. The construction does not require any assumptions on the structure among discrete variables. We demonstrate our method on two examples in genetics based on the Kingman coalescent, showing that the zig-zag process can lead to efficiency gains of up to several orders of magnitude over classical Metropolis-Hastings algorithms, and that it is well suited to parallel computation. Our construction resembles existing techniques for Hamiltonian Monte Carlo on a hybrid state space, which suffers from implementationally and analytically complex boundary crossings when applied to the coalescent. We demonstrate that the continuous-time zig-zag process avoids these complications.


Introduction
The zig-zag process is a non-reversible, piecewise deterministic Markov process introduced by [BR17,BFR19b] for continuous-time MCMC.It has several advantages over reversible methods such as Metropolis-Hastings [Has70] and Gibbs sampling [GS90]: it avoids diffusive backtracking which slows their mixing, and is rejection-free so that no computation is wasted on rejected moves.
In brief, the generator of the zig-zag process (x t , v t ) t≥0 targeting the probability density (with respect to the product of the Lebesgue and counting measures) π(x, v) := π(x)/2 d on a state space where ∂ i is the derivative with respect to x i , and F i flips the sign of v i .The flip rates λ i (x, v) := (−v i ∂ i log(π(x, v))) +  (1) with (x) + := max{x, 0}, ensure that (x t , v t ) t≥0 leaves π(x, v) invariant [BFR19b, Theorem 2.2].
In words, the coordinates of x t move at constant velocities v t until a flip at coordinate i, when the corresponding velocity changes sign.
1 arXiv:2004.08807v4[stat.CO] 18 Jan 2022 To date, the zig-zag processes have been applied to targets such as the Curie-Weiss model [BR17] and logistic regression [BFR19b], whose state spaces have simple geometric structures with natural notions of direction and velocity.Discrete variables (other than the velocities) have been restricted to cases with simple partial orders, such as model selection [CFS20,GD21].We construct a zig-zag process on a general hybrid state space with both continuous and discrete coordinates, without imposing any structure on discrete coordinates.This is done by introducing a separate space of continuous variables for each value of the discrete variable, introducing boundaries into the continuous spaces, and updating the discrete variable when boundaries are hit.The strategy takes advantage of the full generality of piecewise-deterministic Markov processes [Dav93, Section 24], and is resembles similar work for Hamiltonian Monte Carlo (HMC) [DBZM17,NDL20].Our method can also been seen of as a generalization of the zig-zag process on a restricted domain [BBCD + 18] to a union of many restricted sub-domains, with jumps between sub-domains at boundary hitting events.We illustrate our method with an application to the coalescent [Kin82]: a tree-valued target with continuous branch lengths, discrete tree topologies with no natural partial order, and no canonical geometric structure.
The coalescent examples illustrate the need for methods which are implementable on complex state spaces.They are also of interest because existing MCMC algorithms for coalescents tend to mix slowly.The key difficulty lies in designing Metropolis-Hastings proposal distributions which combine efficient exploration with a high acceptance rate [MV05, HDPD08, LVH + 08].Workarounds consist of empirical searches for efficient proposals [HD12,ASR16] or Metropoliscoupled MCMC [Gey92].The former does not scale to problems for which empirical optimization is infeasible.The latter helps mixing between modes, but does not overcome low acceptance rates or the backtracking behavior of reversible MCMC.
The zig-zag process has some similarities with HMC [Nea10], which augments the state space with momentum and uses Hamiltonian dynamics to propose large steps which are accepted with high probability, though they are not rejection-free.Like the zig-zag process, HMC requires gradient information and a suitable geometric embedding of the target.[DBZM17] provided those for the coalescent using an orthant complex construction of phylogenetic tree space [BHV01].
Our examples differ from the method of [DBZM17] in three ways.Firstly, we replace the embedding of [BHV01] with τ -space [GD16], which is better suited to ultrametric trees.Secondly, the zig-zag process is readily implementable on τ -space via Poisson thinning without a numerical integrator such as leap-prog [DBZM17, Algorithm 1].Finally, the zig-zag process has simple boundary behavior between orthants and does not require boundary smoothing [DBZM17, Section 3.3], chiefly because discontinuous gradients are easier to handle on continuous rather than discretized paths.
The rest of the paper is structured as follows.Section 2 defines the zig-zag algorithm with discrete and continuous variables, and proves that it has the desired invariant distribution.Section 3 recalls the coalescent and the τ -space embedding.In Sections 4 and 5 we recall the popular infinite and finite sites models of mutation, derive zig-zag processes for each, and demonstrate their performance via simulation studies.Section 6 concludes with a discussion.The algorithms and data sets used in the simulation studies are available at https://github.com/JereKoskela/tree-zig-zag.

Zig-zag on hybrid spaces
The definition of our zig-zag process follows [Dav93, Section 24].Let F be a countable set.
. Integrals over Ω o and ∂Ω, or subsets thereof, are taken to incorporate discrete sums in the m ∈ F coordinate.
By [Dav93,Theorem 26.14], the zig-zag process with dynamics defined above is a piecewisedeterministic Markov process with extended generator whose domain D(L) consists of measurable functions f (m, x, v) on Ω * satisfying: ) is integrable for each n ∈ N, where {T k } k≥1 are the successive jump times (both velocity flips and jumps due to hitting a boundary) of (m t , x t , v t ) t≥0 .
For a set A, let B(A) and C 1 (A) be the respective spaces of bounded and continuously differentiable functions on A. Let C 1 b (A) := B(A) ∩ C 1 (A).For t > 0, we define Suppose the initial distribution of (m, x) has a density on Ω * and that (2), (3), and (4) hold.
Then the zig-zag process with generator (5) and domain D(L) as described above has stationary distribution π.
Proof.Provided in the appendix.
Remark 1.In addition to having the right invariant distribution, a practical algorithm needs to be ergodic.To discuss ergodicity of the zig-zag process (m t , x t , v t ) t≥0 from Theorem 1, let {(m j t , x j t , v j t ) t≥0 } j∈F be zig-zag processes restricted to respective spaces Ω o j by boundary jump kernels ), each with target proportional to π(j, •, •).When F is finite, a sufficient condition for ergodicity of the global process (m t , x t , v t ) t≥0 is that {(m j t , x j t , v j t ) t≥0 } j∈F are all ergodic, and that for ordered pairs (m, j) ∈ F 2 which form a cycle spanning the support of π. [BRZ19] provide conditions for ergodicity of single-domain zig-zag processes.
We conclude this section with a pseudocode specification of our zig-zag algorithm.

6:
Set ρ to be the solution of Set τ ← ρ and I ← i 9: 14: end while

The coalescent and a geometric embedding
An ultrametric binary tree with n labeled leaves is a rooted binary tree in which each leaf is an equal graph distance away from the root.We follow [GD16] and encode such a tree with leaf labels {1, . . ., n} as the pair (E n , t n ), where E n is a ranked topology and t n ∈ (0, ∞) n−1 .The continuous variables t n encode times between mergers.The time from the leaves to the first merger is t 1 , and subsequent t i variables are times between successive mergers.The ranked topology E n is an (n − 1)-tuple of pairs of labels, where the ith pair specifies the two child nodes of the ith merger.Non-leaf nodes are labeled by the leaves they subtend.For example, the ranked topology encodes the four leaf caterpillar tree depicted in Figures 3 and 7 with nodes labeled left to right.We order the two entries of E n,i by their least element for definiteness.
The coalescent [Kin82] is a seminal model for the genetic ancestry of samples from large populations.Under the coalescent, a tree (E n , t n ) has probability density which arises as the law of a tree constructed by starting a lineage from each leaf, and merging each pair of lineages at rate 1 until the most recent common ancestor (MRCA) is reached.
The success of the coalescent is due to robustness: distributions of ancestries of a large class of individual-based models converge to the coalescent in the infinite population limit under suitable rescalings of time.For details, see e.g.[Wak09].
We define swap operators In words, s i swaps the order of the (i − 1)th and ith mergers.We also define pivot operators p ↓ i and p ↑ i for i ∈ {2, . . ., n − 1} : where ) is the entry of E n,i with the higher (resp.lower) least element, and E s n,i is the sibling: the entry of E n,i that is not E n,i−1 .Pivot p ↑ i is defined by interchanging ↓ and ↑ in (10).The pivots are the two nearest neighbor interchanges between the ith merger and the merger involving its nearest child.Figure 1 illustrates all three operators.
Next we describe τ -space, which gives a geometric structure to the set of pairs (E n , t n ).For fixed E n , the space of t n -vectors is the orthant [0, ∞) n−1 .Each boundary point with t i = 0 corresponds to a degenerate tree in which one of three things happens: 1.The two leaves of E n,1 merge at time 0 if i = 1.

There are two simultaneous mergers if
3. There is a simultaneous merger of three lineages if Type 1 boundaries are boundaries of the whole τ -space.Type 2 boundaries separate orthants corresponding to two ranked topologies separated by an s i -step.Trajectories crossing the boundary move from one ranked topology to the other.Type 3 boundaries separate the three orthants which resolve the triple merger into two binary mergers, which differ by a p ↓ i or p ↑ i -step.A trajectory that crosses the boundary visits a tree with the triple merger.Figure 2 depicts the τ -space with three leaves, and a type 2 boundary in a general τ -space.An example with five leaves is depicted in [GD16, Figure 2], and the t n variables are illustrated in Figures 3 and 7 in Sections 4 and 5. 1 , and s 2 .The horizontal arrangement of leaves is arbitrary throughout this paper; only vertical distance is meaningful.
We will use τ -space to construct zig-zag processes whose state spaces consist of tree topologies, branch lengths, and a scalar parameter introduced in the next section.The discrete variables F will be the ranked topologies, with the boundary crossings described above defining the boundary jump kernel Q.For more on τ -space, e.g.existence and uniqueness of geodesics and Fréchet means, we refer to [GD16].

The infinite sites model
The infinite sites model [Wat75] connects the coalescent tree to DNA sequence data by associating the MRCA with the unit interval (0, 1).Mutations with independent, U (0, 1)-distributed locations accrue along branches of the tree (with branch lengths as specified by t n ) at rate θ/2.The type of a leaf consist of the mutations along branches separating it from the MRCA.We denote the resulting list of types of the leaves by D n .A realization of a coalescent tree and associated D n is shown in Figure 3.The typical task is to sample from the conditional law (E n , t n , θ)|D n corresponding to observing D n , but not the tree which gave rise to it.
To handle mutations, we define F n as the rooted graphical tree with 2n − 1 nodes, the first n of which are leaves labeled 1, . . ., n, while the remaining n − 1 are labeled as in E n .Edges connect children to their parents, and edge lengths are determined by t n .For an edge γ ∈ F n , we denote by c γ and p γ the respective labels of the child and parent nodes of γ, by m γ the number of mutations on γ, and by l γ := t j ∈γ t j the edge length, where we write t i ∈ γ if t i contributes to the length of γ in which case we say γ spans t i .
Given a prior density π 0 (θ) for θ, the posterior distribution of (E n , t n , θ)|D n follows from (9), the fact that the number of mutations on branch γ  and that mutations on distinct branches are independent.In particular, provided E n is consistent D n , and π = 0 otherwise.This distribution can be sampled using a zigzag algorithm by taking F to be the set of ranked topologies on n leaves which are consistent with D n , as well as In the boundary classification of Section 3, θ = 0 is another type 1 boundary.For (t n , θ) ∈ ∂Ω, we define k(t n , θ) as the index of the zero variable, taken to be n in the case of θ.
We augment the state space with n zig-zag velocities v n , of which (v 1 , . . ., v n−1 ) drive t n and v n drives θ.For γ ∈ F n , we also define v γ := j:t j ∈γ v j as the rate of change of l γ .The boundary kernel Q is defined separately on each boundary type: for type 1, At a type 1 boundary the process reflects back into Ω o En via a velocity flip.At a type 2 boundary it undergoes an s i -step and a velocity flip to pass through the boundary.For type 3 boundaries it chooses an adjacent orthant uniformly at random.
In the interiors of orthants, velocity flip rates are Simulating holding times with these rates is difficult due to the time intervals during which they vanish.One strategy is Poisson thinning via dominating rates consisting of only those terms in the round brackets in ( 13) and ( 14) whose sign matches that of the corresponding velocity v i , but these can result in loose bounds and inefficient algorithms.Instead, we define t n := θ for brevity, and for i ∈ {1, . . ., n − 1} define γ(E n , i) := argmin{l γ : p γ = E n,i } as the shorter child branch from parent node E n,i .For i ∈ {1, n} we adopt the short-hands and finally define the time localization T ≡ T (E n , t n , θ; v n ) as where min ∅ = ∞, and K 0 is a maximum increment for the case when all velocities are positive.The indicator functions in the denominator pick out boundaries where (11) vanishes: type 1 or 3 boundaries corresponding to length 0 branches which carry at least one mutation, and the θ = 0 boundary if there is at least one mutation in total.The parameter c > 0 ensures that, when the current process time is t, such boundaries cannot be hit on [t, t + T ], and at most one other boundary can be reached.We found c = 4 gave good performance across our tests in Sections 4 and 5.A larger value results in tighter bounds on (13) and ( 14), but wastes more computation as t + T is hit more often before an accepted velocity flip.
Proof.Provided in the appendix.
We compared the zig-zag process to a Metropolis-Hastings algorithm by reanalyzing the data of [WFDP91] with n = 55, 14 distinct types, and 18 mutations.We used the improper prior π 0 (θ) ∝ 1 (0,∞) (θ), set v i = ±2/[(n + 1 − i)(n − i)], and set v n from trial runs to cross the θmode in unit time.The appendix details the Metropolis-Hastings algorithm, and other tuning parameters in Table 3.We compared both methods to a hybrid combining zig-zag dynamics with continuous time Metropolis-Hastings moves at rate κ = 10.Performance was insensitive to κ provided it was not extreme: small values resemble a zig-zag process, while large values resemble Metropolis-Hastings. Figure 4 shows that the zig-zag and hybrid methods mix visibly better than Metropolis-Hastings over the latent tree, as measured by the tree height H n := t 1 + . . .+ t n−1 .However, they are not as effective at exploring the upper tail of the θ-marginal, likely because they do not stay in regions of short trees for long enough for θ to increase into the tail.
To assess scaling, we simulated two data sets: one of size n = 550 with mutation rate θ = 5.5 (the approximate posterior mean in Figure 4) and one with n = 55 and θ = 55, which models a segment of DNA 10 times longer.Figures 5 and 6 demonstrate that the zig-zag and hybrid processes scale far better than Metropolis-Hastings, particularly when θ = 55.Estimates in Table 1 quantify the improvement to 1-3 orders of magnitude.

The finite sites model
The finite sites model [JC69] is more detailed than the infinite sites model, but has greater computational cost.Consider a finite set of sites S with a finite number of possible types H per site; for example H = {0, 1} or H = {A, T, C, G}.Mutations occur along branches of the coalescent tree at each site with rate θ/(2|S|), and the type of a mutant child is drawn from stochastic matrix P with unique stationary distribution ν.We denote the transition matrix of the H-valued compound Poisson process with jump rate θ/(2|S|) and jump transition matrix P by (Q θ hg (t)) h,g∈H;t≥0 .Figure 7 depicts a realization of the finite sites coalescent.As in Section 4, we denote the configuration of types at the leaves by D n , and seek to sample from the posterior π(E n , t n , θ|D n ), which can be written as a sum over the types of internal nodes: where h(s; η) ∈ H is the type at site s ∈ S on the node with label η ∈ E n , and γ ∈ F n : |c γ | > 1 denotes edges which do not end in a leaf.The target ( 16) can be evaluated efficiently using the pruning algorithm of [Fel81].
The posterior (16) can be sampled using zig-zag dynamics with the same construction as in  Section 4. Velocity flip rates can be written in terms of branch-specific gradients as which can be evaluated using the linear-cost method of [JZH + 20].
We show that events with rates ( 17) and ( 18) can be simulated using the example of [GT94, Section 7.4], in which S = {1, . . ., 20}, H = {0, 1} and P is the 2 × 2 matrix which always changes state, corresponding to ν = (1/2, 1/2) and As ( 19) is not bounded away from 0 when h = g, (17) and (18) lack simple bounds for Poisson thinning.As in (15), bounds can be obtained by time localization using where K 0 is a default increment in case all velocities are positive.The variable T localizes the next zig-zag time step beginning at time t so that at most one branch can shrink to length zero on [t, t + T ], θ can fall by at most 1/(1 + c) of its present value, and t 1 can fall by at most 1/(1 + c) of its present value if the first two leaves to merge are of distinct types.This treatment of θ and t 1 is needed as (17) and (18) diverge in these cases, rendering the θ = 0 and t 1 = 0 boundaries inaccessible.
Given T ∈ (0, ∞), we have the following bounds on (19) on the time interval [t, t + T ]: where h = g.Substituting these bounds into [JZH + 20, Equation ( 9)] provides bounds on flip rates that can be evaluated with O(|S|n) cost.
Figure 8 demonstrates that the zig-zag process mixes over latent trees faster than Metropolis-Hastings again, but struggles to explore the upper tail of the θ-marginal.The hybrid method was run with κ = 100 to compensate for shorter run lengths than in Section 4, and thus resembles Metropolis-Hastings rather than the zig-zag process.

Discussion
We have presented a general method for using zig-zag processes to sample targets defined on hybrid spaces consisting of discrete and continuous variables.This was done by introducing boundaries into the state space of continuous variables and updating discrete components via a Markov jump kernel Q whenever a boundary was hit.The resulting algorithm remains a piecewise-deterministic Markov processes in the sense of [Dav93, Section 24], and generalizes We have demonstrated the method on two examples involving the coalescent, which is a goldstandard model in phylogenetics.It is defined on the space of binary trees consisting of discrete tree topologies and continuous branch lengths, which lacks a simple geometric structure, e.g. a partial order or a tractable norm.We have also shown that the zig-zag process can improve mixing over trees relative to Metropolis-Hastings, particularly under the infinite sites model.This model is widely used to analyze ever larger data sets, and the zig-zag process shows promise for expanding the scope of feasible MCMC computations.
The zig-zag process was more efficient than Metropolis-Hastings under the infinite sites model in terms of effective sample size, but struggled to explore the tails of the θ-marginal.A likely reason is correlation in the target: high mutation rates are only be attainable when branch lengths are short.A Metropolis-Hastings algorithm can jump to a high mutation rate as soon as the latent tree has short branches, while the zig-zag process must traverse all intervening mutation rates before branch lengths grow.The speed up zig-zag method of [VR21] has statedependent velocities, and could provide further improvement.The hybrid method with both zig-zag motion and Metropolis-Hastings updates interpolated between the two algorithms.
All three algorithms exhibited much longer run times under the finite sites model than under infinite sites.For the zig-zag and hybrid methods, that is due to the O(|S|n) cost per evaluation of ( 17) and (18), of which there are O(n).However, flip times for different velocities are conditionally independent given the current state and can be generated in parallel, unlike f ∈ D(L), The change of variable v → F i v combined with the fact that π(m, x, •) is uniform, then (1), and then an integration by parts yield where ∇ x is the gradient operator with respect to x.To conclude, we will show that the first term on the right hand side of (20) cancels with the second.Using (6), then (2), and then × π(j, y, −w)dy.
We now exchange order of integration (justified by the fact that f is bounded), then apply (3) and the fact that π(x, •) is uniform to obtain Substituting ( 21) into (20) concludes the proof.
Proof of Proposition 1.The proof consists of checking the hypotheses of Theorem 1.The set of ranked topologies of binary trees with n leaves is finite, and hence so is its restriction to topologies consistent with D n .The initial condition has a density by assumption.The set of points from which deterministic motion at a constant velocity first hits a boundary at a corner consists of finitely many Lebesgue-null lines.Hence, since the initial condition has a density and the boundary kernel Q does not change (t n , θ), corners are hit with probability 0.
The augmented posterior π corresponding to (11) is C 1 (Ω * ) by inspection, and only vanishes at a boundary when θ = 0 or a branch with mutations has length zero.
It remains to verify (7), which we do by stochastic domination.By construction, the zigzag process crosses boundaries only when a corresponding velocity is negative.The boundary crossing flips the velocity, and a boundary corresponding to the same coordinate cannot be crossed again until a further flip.Hence, the number of jumps in the zig-zag process started at (E n , t n , θ) by time t ∈ (0, ∞) is dominated by n + 2Y (t n , θ, t), where which arises as the volume of the set which is reachable from (E n , t n , θ) by time t multiplied by upper bounds for (13) and ( 14) for positive velocities on the reachable set.Here • 1 is the L 1 -norm, which is invariant for all valid zig-zag velocities v n .In the construction of the dominating random variable, negative velocities are assumed to flip to positive instantaneously which accounts for the initial summand of n (the maximum number of initial negative velocities) and the factor of 2.
The remainder of the appendix details the zig-zag and Metropolis-Hastings algorithms used in Sections 4 and 5. Algorithms 2 and 3 provide pseudocode implementations of the zig-zag methods.Line 7 of Algorithm 2 is given twice, with the first applicable to Section 4 and the second to Section 5. Otherwise, both algorithms apply to both sections.
The for-loop on lines 4-10 of Algorithm 1 implements the time localization and ensures that time steps cannot cause negative entries in t n .The distribution of the increment τ thus has an atom at the value determined by lines 4-10, signifying either a boundary crossing or a need to recompute the time localization and rates {λ * i : i ∈ {1, . . ., n − 1, θ}}.However, the continuoustime motion of the trajectory is not interrupted in either case; see also [BBCD + 18, Section 2.2].
Iterations of the outer while-loop of Algorithm 1 flip a velocity for one of two reasons: a flip event occurs or a boundary is hit.In the former case the reason for the flip is clear by construction.In the latter, a trajectory arriving at a boundary must have v i < 0 for the corresponding velocity.After crossing, the trajectory moves into the interior of an adjacent orthant, which corresponds to v i > 0 in its local coordinates.

8:
Set α ← 1 9: Sample U ∼ U (0, 1) 10: until U < α 11: return ρ The Metropolis-Hastings algorithms in Sections 4 and 5 both used the same proposal mechanisms consisting of three steps.An iteration of the algorithm is one scan through all three steps, with an accept/reject correction after each step.The hybrid method only used steps 1 and 3. 1.A Gaussian random walk update of θ reflected at 0 with proposal variance σ 2 θ tuned to obtain an average acceptance probability α θ ≈ 1/4.
2. An update of holding times t n under a fixed topology.where ξ c i,1 and ξ c i,2 are the perturbed times of the child nodes of the node at time t i .Again, σ tn was tuned so that the average acceptance probability was α tn ≈ 1/4.Under the infinite sites model, moves into topologies incompatible with observed mutations were rejected essentially instantaneously, before costly likelihood evaluations.
For each m ∈ F, let Ω o m be an open subset of R d , Ω o m be its closure, and ∂Ω * m := Ω o m \ Ω o m be its boundary.We assume that {∂Ω * m } m∈F are piecewise Lipschitz and denote by ∂Ω m the restriction of ∂Ω * m to non-corner points.Let Ω o := ∪ m∈F Ω o m = {(m, x) : m ∈ F, x ∈ Ω o m } and ∂Ω := ∪ m∈F ∂Ω m .For a point (m, x) ∈ ∂Ω m , let n(m, x) be the outward unit normal, and let Γ ± (m, x) := {v ∈ {−1, 1} d : ±(v • n(m, x)) > 0} be the sets of velocities with which a zig-zag process can exit (+) or enter (-) Ω o m at x. Zig-zag dynamics imply

Figure 2 :Figure 3 :
Figure2: (Left) τ -space with n = 3 embedded into R 3 .Each square is a copy of [0, ∞) 2 associated with the given topology.The coordinates (t 1 , t 2 ) are the respective time of the first merger, and the time between the first and second merger.The dot is the origin, and the line on which all three orthants intersect is a type 3 boundary consisting of trees in which all three leaves merge simultaneously at time t 1 .The dashed lines are boundaries at ∞. (Right) A segment of τ -space depicting a type 2 boundary, in which each square represents [0, ∞) n−1 .Only the two orthants adjacent to the boundary are shown.

Figure 4 :
Figure 4: Trace plots under the infinite sites model and the data set of [WFDP91].

Figure 5 :
Figure 5: Trace plots for the infinite sites model and the data set with n = 550, θ = 5.5, 30 distinct types, and 38 mutations.

Figure 6 :
Figure 6: Trace plots for the infinite sites model and the data set with n = 55, θ = 55, 40 distinct types, and 252 mutations.

Figure 8 :
Figure 8: Trace plots for the finite sites model and data from [GT94].

Figure 9 :
Figure 9: Trace plots for the finite sites model and data set with n = 500 and S = 20 consisting of five distinct sequences.

Figure 10 :
Figure 10: Trace plots for the finite sites model and data set with n = 50 and S = 200 consisting of 18 distinct sequences.

Table 1 :
[FHVD17]e sample sizes and run times for all three methods and data sets for the infinite sites model.Estimates were computed with the ess method[FHVD17]under default settings for Metropolis-Hastings, and as in [BFR19a, Section 2] for the other two.Stars highlight poorly mixing chains with unreliable ESS estimates.

Table 2 :
Effective sample sizes and run times for all three methods and data sets.Estimates were computed as in Table1.Stars highlight unmixed chains with unreliable ESS estimates.
Type 3 updates sample an ordered pair of edges (γ, γ) ∼ U (F n × [F n ∪ γ M RCA ]), where γ M RCA is an edge connecting the MRCA to ∞. Edge γ is deleted, and its child c γ is reattached into γ .Letting t η be the time of the node with label η ∈ E n in an abuse of notation, the reattachment time t has distributiont ∼ U (t cγ ∧ t c γ , t p γ ) if γ = γ M RCA , t − t cγ M RCA ∼ Exp(1) otherwise.

Table 3
summarizes hyperparameters and quantities of interest used in the simulation.