Source code for moments.Demographics3D

"""
Three-population demographic models.
"""
import numpy as np
import moments


[docs]def out_of_Africa(params, ns, pop_ids=["YRI", "CEU", "CHB"]): """ The Gutenkunst et al (2009) out-of-Africa that has been reinferred a number of times. :param params: List of parameters, in the order (nuA, TA, nuB, TB, nuEu0, nuEuF, nuAs0, nuAsF, TF, mAfB, mAfEu, mAfAs, mEuAs). :type params: list of floats :param ns: List of population sizes in each population, in order given by `pop_ids`. :type ns: list of ints :param pop_ids: List of population IDs, defaults to ["YRI", "CEU", "CHB"]. :type pop_ids: list of strings, optional """ if pop_ids is not None and len(pop_ids) != 3: raise ValueError("pop_ids must be a list of three population IDs") if len(ns) != 3: raise ValueError("ns must have length 3") ( nuA, TA, nuB, TB, nuEu0, nuEuF, nuAs0, nuAsF, TF, mAfB, mAfEu, mAfAs, mEuAs, ) = params sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1] + ns[2]) fs = moments.Spectrum(sts) # integrate modern human branch with expansion fs.integrate([nuA], TA) # split into African and Eurasian branches and integrate fs = fs.split(0, ns[0], ns[1] + ns[2]) mig_mat = [[0, mAfB], [mAfB, 0]] fs.integrate([nuA, nuB], TB, m=mig_mat) # split Eurasian into CEU and CHB fs = fs.split(1, ns[1], ns[2]) nu_func = lambda t: [ nuA, nuEu0 * np.exp(np.log(nuEuF / nuEu0) * t / TF), nuAs0 * np.exp(np.log(nuAsF / nuAs0) * t / TF), ] mig_mat = [[0, mAfEu, mAfAs], [mAfEu, 0, mEuAs], [mAfAs, mEuAs, 0]] fs.integrate(nu_func, TF, m=mig_mat) fs.pop_ids = pop_ids return fs