“Custom” multivariate distributions

Show code cell source Hide code cell source

import demes
import demesdraw
import fwdpy11
import numpy as np

The previous two sections cover cases where the methods for generating deviates from a multivariate distribution are straightforward and agreed upon.

In order to simulate multivariate distributions of effect sizes based on fwdpy11.Sregion types, we follow a fairly intuitive approach described in [XKS00]. Briefly, the multivariate Gaussian kernel is used to produce deviates. Then, the quantiles from the cummulative distribution of each marginal Gaussian are used to generate a deviate from the desired output distribution of interest.

For a simulation with n populations we need:

• A list of n fwdpy11.Sregion objects

• An array of n means for the multivariate Gaussian

• An n-by-n covariance matrix for the multivariate Gaussian

We will use the same demographic model as in the previous vignettes. The details can be viewed by expanding the “click to show” button below”.

Show code cell source Hide code cell source

yaml = """
description: Island model forever
time_units: generations
demes:
  - name: A
    epochs:
     - start_size: 100
  - name: B
    epochs:
     - start_size: 100
migrations:
  - demes: [A, B]
    rate: 0.10
"""
g = demes.loads(yaml)
model = fwdpy11.discrete_demography.from_demes(g, burnin=1)
initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
pdict = {
    "nregions": [],
    "recregions": [],
    "sregions": None, # Will get filled in below
    "rates": (0, 5e-3, None),
    "demography": model,
    "simlen": model.metadata["total_simulation_length"],
    "gvalue": fwdpy11.Multiplicative(ndemes=2, scaling=2),
}
rng = fwdpy11.GSLrng(123512)

The following generates exponentially distributed effect sizes in each deme with a high correlation across demes:

mvdes = fwdpy11.mvDES(
    [fwdpy11.ExpS(0, 1, 1, -0.5)] * 2,
    np.zeros(2),
    np.matrix([1, 0.9, 0.9, 1]).reshape((2, 2)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.0481598  -0.01990405]
[-1.31708542 -1.10171055]
[-0.0336445 -0.0317481]
[-0.14208652 -0.10416488]

We can mix and match our distributions. Here, the distribution of effect sizes in deme 0 is exponential and the distribution in deme 1 is gamma. The two distributions have means with opposite signs and the magnitudes of the marginal deviates negatively covary:

mvdes = fwdpy11.mvDES(
    [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.GammaS(0, 1, 1, mean=0.1, shape_parameter=1)],
    np.zeros(2),
    np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.33083002  0.10031833]
[-0.42548846  0.05563312]
[-0.13448292  0.25725126]
[-0.05196185  0.15740337]
[-0.13115791  0.08624639]
[-0.10758304  0.14649511]
[-0.24521457  0.16010134]
[-0.06690352  0.19844306]
[-0.2526577   0.09564946]
[-0.29293677  0.05486382]
[-0.13782504  0.091609  ]
[-0.37426508  0.08506987]
[-0.26032863  0.10840349]
[-0.50714586  0.08579448]
[-0.17467581  0.08643293]
[-0.20649866  0.09632986]
[-0.48077505  0.13611904]
[-0.00917911  0.44995212]
[-0.16189138  0.13520003]
[-1.20496755  0.01927578]
[-1.65766443  0.00423084]
[-0.12988024  0.12932187]
[-0.0454018   0.23482635]
[-0.19691462  0.07874511]
[-0.00514447  0.41722543]

The type fwdpy11.ConstantS has intuitive behavior:

mvdes = fwdpy11.mvDES(
    [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.ConstantS(0, 1, 1, -0.1)],
    np.zeros(2),
    np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
)
pdict["sregions"] = [mvdes]
params = fwdpy11.ModelParams(**pdict)
pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
rng = fwdpy11.GSLrng(1010)
fwdpy11.evolvets(rng, pop, params, 10)
for i in pop.tables.mutations:
    print(pop.mutations[i.key].esizes)

[-0.04863043 -0.1       ]
[-0.60824008 -0.1       ]
[-0.00696691 -0.1       ]
[-0.08694534 -0.1       ]
[-0.01399284 -0.1       ]
[-0.15168546 -0.1       ]
[-0.01237223 -0.1       ]
[-0.45890092 -0.1       ]
[-0.45196062 -0.1       ]
[-0.16847123 -0.1       ]
[-0.29808442 -0.1       ]