# Simulating trees from the fossilized birth-death process

#### Walker Pett

#### 2026-02-26

## `sim.fbd` functions

`FossilSim` can be used to simulate trees directly from a fossilized-birth death process using the `sim.fbd` family of functions, whose interface is similar to the `sim.bd` functions in `TreeSim`. These functions return a list of `SAtree` objects, since the simulated trees may contain sampled ancestors represented as tips with zero-length branches.

## Simulating trees from a time-homogeneous fossilized-birth death process

As in `TreeSim`, trees can be simulated while conditioning either on the number of sampled extant taxa (using the `sim.fbd.taxa` function), or on the total age of the process (using the `sim.fbd.age` function). These functions require the fossil sampling rate parameter `psi` as an argument, but otherwise take the same arguments as the `sim.bd.taxa` and `sim.bd.age` functions in `TreeSim`.

    trees <- sim.fbd.taxa(n = 10, numbsim = 10, lambda = 3, mu = 2, psi = 2)

    plot(trees[[1]])

A stratigraphic range representation of the tree and fossil samples can be plotted using the `rangeplot.asymmetric` function, which assumes that all speciation events are asymmetric.

    rangeplot.asymmetric(trees[[1]])

Note that when simulating a complete tree (`complete = TRUE`), tips in the resulting `SAtree` object represent extinction events, whereas they represent fossil sampling events in the tree without unsampled lineages (`complete = FALSE`). The `complete` field of resulting the `SAtree` object is set to the value of the `complete` function argument.

    trees <- sim.fbd.taxa(n = 10, numbsim = 10, lambda = 3, mu = 2, psi = 2, complete = T)

    print(trees[[1]]$complete)

    ## [1] TRUE

    rangeplot.asymmetric(trees[[1]])

## Simulating trees from an episodic fossilized-birth death process

Trees can also be simulated from a time-varying, piecewise constant fossilized birth-death process using the `sim.fbd.rateshift.taxa` function. Apart from a vector of interval specific `psi` parameters, this function also takes most of the same arguments as the corresponding `sim.rateshift.taxa` function in `TreeSim`. At the moment, this function assumes the sampling fraction of extant taxa (`frac`) is equal to 1.0 in all time intervals.

    trees = sim.fbd.rateshift.taxa(n = 10, numbsim = 1, lambda = c(2,1), mu = c(0,0.3), psi = c(1,0.1), times = c(0,0.3))

    rangeplot.asymmetric(trees[[1]])