Using rewards

Using rewards#

In order to compute more complex moments such as higher order (cross)-moments that are not directly made available as cached properties of PhaseTypeDistribution, we can specify our own rewards. A Reward is a means of rewarding or weighting each state so as to obtain the moments of the quantity of interest. Examples of common rewards are TreeHeightReward and TotalBranchLengthReward or UnfoldedSFSReward. We can use PhaseTypeDistribution’s moment(), which requires a tuple of rewards to be specified, whose length equals the order of the moment to be computed.

import phasegen as pg

coal = pg.Coalescent(n=10)

We can explicitly compute the moments made available as cached properties. Note that moment() provides centered moments by default.

# mean tree height
assert coal.tree_height.mean == coal.moment(1, (pg.TreeHeightReward(),))

# variance of tree height
assert coal.moment(2, (pg.TreeHeightReward(),) * 2) == coal.tree_height.var

# second non-central moment of total branch length
assert coal.total_branch_length.m2 == coal.moment(2, (pg.TotalBranchLengthReward(),) * 2, center=False)

# mean of the 2nd unfolded SFS entry
assert coal.sfs.mean.data[2] == coal.moment(1, (pg.UnfoldedSFSReward(2),))

If necessary, we can also compute much higher order (cross)-moments.

# 5th non-central moment of tree height
coal.moment(5, (pg.TreeHeightReward(),) * 5)

48.26441927101328

# 3rd non-central moment of 2nd, 3rd and 4th unfolded SFS entries
coal.moment(3, (pg.UnfoldedSFSReward(2), pg.UnfoldedSFSReward(3), pg.UnfoldedSFSReward(4)))

0.005538770711106311

Combining rewards#

Sometimes we may want to combine multiple rewards. To this end, we can use ProductReward or SumReward. Let us for example compute the mean SFS over the first 2 out of 3 demes. That is, we weight the branch lengths by the fraction of lineages in the first two demes.

# 3-deme coalescent with symmetric migration
coal = pg.Coalescent(
    n={'pop_0': 3, 'pop_1': 2, 'pop_2': 0},
    demography=pg.Demography(
        pop_sizes={'pop_0': 3, 'pop_1': 0.5, 'pop_2': 0.1},
        events=[
            pg.SymmetricMigrationRateChanges(
                pops=['pop_0', 'pop_1', 'pop_2'],
                rate=1
            )
        ]
    )
)

sfs = coal.sfs.moment(1, (pg.SumReward([pg.DemeReward('pop_0'), pg.DemeReward('pop_1')]),))

sfs.plot();

<Figure size 440x330 with 0 Axes>
../../_images/b56a855c919f45087933df72719cc61f00963e91932e7868878ed87157c9613e.png

Plotting the marginal spectra for each deme, we indeed see that the mean SFS over the first two demes was obtained above.

pg.Spectra({d: coal.sfs.demes[d].mean for d in coal.demography.pop_names}).plot();

<Figure size 440x330 with 0 Axes>
../../_images/aeeb7ce1e024921a6988436d44fc3d1218e3e63697adecebb47075e7dde33631.png

Note that we have here used moment() of UnfoldedSFSDistribution. This method internally uses a ProductReward to combine the given reward with UnfoldedSFSReward to obtain the SFS for all possible SFS bins. Had we used moment() instead, we would have needed to specify the reward for each SFS bin separately.

demes = pg.SumReward([pg.DemeReward('pop_0'), pg.DemeReward('pop_1')])
sfs_bin = pg.UnfoldedSFSReward(2)

# take product of SFS reward for second SFS bin with sum of rewards for first two demes
assert sfs.data[2] == coal.moment(1, (pg.ProductReward([demes, sfs_bin]),))

Tracing reward accumulation over time#

Since rewards are accumulated over time, we can trace their accumulation over time. This can be useful for debugging purposes or to understand how the reward is accrued over time. To this end, we can use PhaseTypeDistribution’s accumulate() and plot_accumulation(). Note that this is different from the CDF, which instead accumulates probability mass over time.

# accumulation of mean SFS over all demes
coal.sfs.plot_accumulation(k=1);
../../_images/c2242c9bbfa4848da5af37daa4ca8d6c88fe6903945dc08fc78582e89f0ef520.png

Adjusting start and end times of reward accumulation#

By default, rewards are accumulated from time 0 until time of almost sure absorption. We can adjust the start and end times of reward accumulation by specifying the start_time and end_time arguments to Coalescent or PhaseTypeDistribution’s moment().

mean1 = pg.Coalescent(n=10, end_time=2).tree_height.mean
mean2 = pg.Coalescent(n=10).tree_height.moment(1, end_time=2)

assert mean1 == mean2
Contents