Python MatrixUtil.get_stationary_distribution方法代码示例

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_two_allele_distribution(N_big, N_small, f0, f1, f_subsample):
    """
    Assumes small genic selection.
    Assumes small mutation.
    The mutational bias does not affect the distribution.
    @param N_big: total number of alleles in the population
    @param N_small: number of alleles sampled from the population
    @param f0: fitness of allele 0
    @param f1: fitness of allele 1
    @param f_subsample: subsampling function
    @return: distribution over all non-fixed population states
    """
    # construct a transition matrix
    nstates = N_big + 1
    P = np.zeros((nstates, nstates))
    for i in range(nstates):
        p0, p1 = wrightfisher.genic_diallelic(f0, f1, i, N_big - i)
        if i == 0:
            P[i, 1] = 1.0
        elif i == N_big:
            P[i, N_big - 1] = 1.0
        else:
            for j in range(nstates):
                logp = StatsUtil.binomial_log_pmf(j, N_big, p0)
                P[i, j] = math.exp(logp)
    # find the stationary distribution
    v = MatrixUtil.get_stationary_distribution(P)
    MatrixUtil.assert_distribution(v)
    if not np.allclose(v, np.dot(v, P)):
        raise ValueError('expected a left eigenvector with eigenvalue 1')
    # return the stationary distribution conditional on dimorphism
    print v
    distn = f_subsample(v, N_small)
    return distn[1:-1] / np.sum(distn[1:-1])

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
 def test_exact_fixation_probability(self):
     s = 0.03
     N_diploid = 3
     N_haploid = N_diploid * 2
     # Define a transition matrix with boundary conditions.
     # This approach is nice because the computationally intensive part
     # is just finding the equilibrium distribution of a process,
     # and you also get the conditional equilibrium distribution
     # of the dimorphic states.
     P = np.exp(wfengine.create_genic_diallelic(N_diploid, s))
     P[0, 0] = 0
     P[0, 1] = 1
     P[N_haploid, N_haploid] = 0
     P[N_haploid, 1] = 1
     v = MatrixUtil.get_stationary_distribution(P)
     fpa = v[-1] / (v[0] + v[-1])
     # Define a transition matrix with boundary conditions.
     # This approach is nice because it gives fixation probabilities
     # conditional on each possible initial allele frequency.
     P = np.exp(wfengine.create_genic_diallelic(N_diploid, s))
     A = P - np.eye(N_haploid + 1)
     b = np.zeros(N_haploid + 1)
     A[0, 0] = 1
     A[N_haploid, N_haploid] = 1
     b[0] = 0
     b[N_haploid] = 1
     x = linalg.solve(A, b)
     fpb = x[1]
     # Compare the transition probabilities.
     #print fpa
     #print get_pfix_approx(N_diploid, s)
     self.assertTrue(np.allclose(fpa, fpb))

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_large_population_approximation(n, k, gammas, M):
    """
    @param n: sample size
    @param k: number of alleles e.g. 4 for A, C, G, T
    @param gammas: scaled selections near zero, positive is better
    @param M: mutation rate matrix
    @return: xxx
    """
    # Approximate the fixation probabilities.
    F = get_scaled_fixation_probabilities(gammas)
    # Compute the rate matrix as the hadamard product
    # of mutation rates and fixation probabilities,
    # and adjust the diagonal.
    Q = M*F
    Q -= np.diag(np.sum(Q, 1))
    # This is kind of a hack,
    # I should just get the stationary distribution directly from Q
    # without the expm.
    v = MatrixUtil.get_stationary_distribution(linalg.expm(Q))
    # Get sample allele frequencies associated with the
    # transitions between the fixed states of alleles 0 and 1.
    m0 = v[0] * M[0,1]
    m1 = v[1] * M[1,0]
    g = gammas[1] - gammas[0]
    return diallelic_approximation(n, g, m0, m1)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
 def test_expected_generations_until_absorption_no_selection(self):
     s = 0
     N_diploid = 100
     N_haploid = N_diploid * 2
     # Define a transition matrix with boundary conditions.
     # This approach is nice because the computationally intensive part
     # is just finding the equilibrium distribution of a process,
     # and you also get the conditional equilibrium distribution
     # of the dimorphic states.
     P = np.exp(wfengine.create_genic_diallelic(N_diploid, s))
     P[0, 0] = 0
     P[0, 1] = 1
     P[N_haploid, N_haploid] = 0
     P[N_haploid, 1] = 1
     v = MatrixUtil.get_stationary_distribution(P)
     fpa = v[-1] / (v[0] + v[-1])
     # Use the Kimura approximation to get the expected
     # number of generations until absorption for a neutral allele.
     p = 1 / float(N_haploid)
     #p_fixation = v[-1] / (v[0] + v[-1])
     #p_loss = v[0] / (v[0] + v[-1])
     # expected number of generations until fixation excluding loss
     t1 = -(1/p)*2*N_haploid*(1-p)*math.log(1-p)
     # expected number of generations until loss excluding fixation
     p_fixation = get_pfix_approx(N_diploid, s)
     p_loss = 1 - p_fixation
     t0 = -2*N_haploid*(p/(1-p))*math.log(p)
     t = p_loss * t0 + p_fixation * t1

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_plot_array(N_diploid, Nr, theta_values, Ns_values):
    """
    Compute expected hitting times.
    Theta is 4*N*mu, and the units of time are 4*N*mu generations.
    @param N_diploid: diploid population size
    @param Nr: recombination rate
    @param theta_values: mutation rates
    @param Ns_values: selection values
    @return: arr[i][j] gives time for Ns_values[i] and theta_values[j]
    """
    # set up the state space
    k = 4
    M = multinomstate.get_sorted_states(2 * N_diploid, k)
    T = multinomstate.get_inverse_map(M)
    nstates = M.shape[0]
    lmcs = wfengine.get_lmcs(M)
    # precompute rate matrices
    R_rate = wfcompens.create_recomb(M, T)
    M_rate = wfcompens.create_mutation(M, T)
    # precompute a recombination probability matrix
    R_prob = linalg.expm(Nr * R_rate / float((2 * N_diploid) ** 2))
    #
    arr = []
    for theta in theta_values:
        # Compute the expected number of mutation events per generation.
        mu = theta / 2
        # Precompute the mutation matrix
        # and the product of mutation and recombination.
        M_prob = linalg.expm(mu * M_rate / float(2 * 2 * N_diploid))
        MR_prob = np.dot(M_prob, R_prob)
        #
        row = []
        for Ns in Ns_values:
            s = Ns / float(N_diploid)
            lps = wfcompens.create_selection(s, M)
            S_prob = np.exp(wfengine.create_genic(lmcs, lps, M))
            P = np.dot(MR_prob, S_prob)
            # compute the stationary distribution
            v = MatrixUtil.get_stationary_distribution(P)
            # compute the transition matrix limit at time infinity
            # P_inf = np.outer(np.ones_like(v), v)
            # compute the fundamental matrix Z
            # Z = linalg.inv(np.eye(nstates) - (P - P_inf)) - P_inf
            #
            # Use broadcasting instead of constructing P_inf.
            Z = linalg.inv(np.eye(nstates) - (P - v)) - v
            # compute the hitting time from state AB to state ab.
            i = 0
            j = 3
            hitting_time_generations = (Z[j, j] - Z[i, j]) / v[j]
            hitting_time = hitting_time_generations * theta
            row.append(hitting_time)
        arr.append(row)
    return arr

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def moran_distn_helper(M, T, R_mut):
    """
    @param M: index to states
    @param T: states to index
    @param R_mut: scaled mutation rate matrix
    @return: stationary distribution of the process
    """
    N = np.sum(M[0])
    R_drift = 0.5 * get_moran_drift(M, T) / float(N)
    R = R_mut + R_drift
    P = scipy.linalg.expm(R)
    v = MatrixUtil.get_stationary_distribution(P)
    return v

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def wright_distn_helper(M, T, R_mut):
    """
    @param M: index to states
    @param T: states to index
    @param R_mut: scaled mutation rate matrix
    @return: stationary distribution of the process
    """
    #FIXME: remove dependence on T
    lmcs = wrightcore.get_lmcs(M)
    lps = wrightcore.create_selection_neutral(M)
    log_drift = wrightcore.create_neutral_drift(lmcs, lps, M)
    P_drift = np.exp(log_drift)
    P_mut = scipy.linalg.expm(R_mut)
    P = np.dot(P_mut, P_drift)
    v = MatrixUtil.get_stationary_distribution(P)
    return v

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
 def __call__(self, X):
     """
     @param X: six params defining mutation and selection
     @return: negative log likelihood
     """
     # define the hardcoded number of alleles
     k = 4
     # unpack the params
     params = X.tolist()
     theta, ka, kb, g0, g1, g2 = params
     if any(x < 0 for x in (theta, ka, kb)):
         return float('inf')
     mutation, fitnesses = kaizeng.params_to_mutation_fitness(
             self.N, params)
     # get the transition matrix
     P = kaizeng.get_transition_matrix(self.N, k, mutation, fitnesses)
     v = MatrixUtil.get_stationary_distribution(P)
     return -StatsUtil.multinomial_log_pmf(v, self.observed_counts)

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_two_allele_distribution(N_diploid, f0, f1):
    """
    Assumes small genic selection.
    Assumes small mutation.
    The mutational bias does not affect the distribution.
    @param N_diploid: diploid population size
    @param f0: fitness of allele 0
    @param f1: fitness of allele 1
    @return: distribution over all non-fixed population states
    """
    #"""
    # get the transition matrix without mutation
    s = 1 - f1 / f0
    P = np.exp(wfengine.create_genic_diallelic(N_diploid, s))
    # add mutation
    P[0, 0] = 0
    P[0, 1] = 1
    P[2*N_diploid, 2*N_diploid] = 0
    P[2*N_diploid, 2*N_diploid-1] = 1
    #"""
    """
    # construct something like a transition matrix
    N_haploid = N_diploid * 2
    nstates = N_haploid + 1
    P = np.zeros((nstates, nstates))
    for i in range(nstates):
        p0, p1 = wrightfisher.genic_diallelic(f0, f1, i, N_haploid-i)
        if i == 0:
            P[i, 1] = 1.0
        elif i == N_haploid:
            P[i, N_haploid-1] = 1.0
        else:
            for j in range(nstates):
                logp = StatsUtil.binomial_log_pmf(j, N_haploid, p0)
                P[i, j] = math.exp(logp)
    """
    # find the stationary distribution
    v = MatrixUtil.get_stationary_distribution(P)
    if not np.allclose(v, np.dot(v, P)):
        raise ValueError
    # return the stationary distribution conditional on dimorphism
    return v[1:-1] / np.sum(v[1:-1])

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def do_full_simplex_then_collapse(mutrate, popsize):
    #mutrate = 0.01
    #mutrate = 0.2
    #mutrate = 10
    #mutrate = 100
    #mutrate = 1
    N = popsize
    k = 4
    M = np.array(list(multinomstate.gen_states(N, k)), dtype=int)
    T = multinomstate.get_inverse_map(M)
    # Create the joint site pair mutation rate matrix.
    R = mutrate * wrightcore.create_mutation(M, T)
    # Create the joint site pair drift transition matrix.
    lmcs = wrightcore.get_lmcs(M)
    lps = wrightcore.create_selection_neutral(M)
    log_drift = wrightcore.create_neutral_drift(lmcs, lps, M)
    # Define the drift and mutation transition matrices.
    P_drift = np.exp(log_drift)
    P_mut = scipy.linalg.expm(R)
    # Define the composite per-generation transition matrix.
    P = np.dot(P_mut, P_drift)
    # Solve a system of equations to find the stationary distribution.
    v = MatrixUtil.get_stationary_distribution(P)
    for state, value in zip(M, v):
        print state, value
    # collapse the two middle states
    nstates_collapsed = multinomstate.get_nstates(N, k-1)
    M_collapsed = np.array(list(multinomstate.gen_states(N, k-1)), dtype=int)
    T_collapsed = multinomstate.get_inverse_map(M_collapsed)
    v_collapsed = np.zeros(nstates_collapsed)
    for i, bigstate in enumerate(M):
        AB, Ab, aB, ab = bigstate.tolist()
        Ab_aB = Ab + aB
        j = T_collapsed[AB, Ab_aB, ab]
        v_collapsed[j] += v[i]
    for state, value in zip(M_collapsed, v_collapsed):
        print state, value
    # draw an equilateral triangle
    #drawtri(M_collapsed, T_collapsed, v_collapsed)
    #test_mesh()
    return M_collapsed, T_collapsed, v_collapsed

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def do_collapsed_simplex(scaled_mut, N):
    """
    @param N: population size
    """
    k = 3
    M = np.array(list(multinomstate.gen_states(N, k)), dtype=int)
    T = multinomstate.get_inverse_map(M)
    # Create the joint site pair mutation rate matrix.
    # This is scaled so that there are about popsize mutations per generation.
    R_mut_raw = wrightcore.create_mutation_collapsed(M, T)
    R_mut = (scaled_mut / float(N)) * R_mut_raw
    # Create the joint site pair drift transition matrix.
    lmcs = wrightcore.get_lmcs(M)
    lps = wrightcore.create_selection_neutral(M)
    #log_drift = wrightcore.create_neutral_drift(lmcs, lps, M)
    # Define the drift and mutation transition matrices.
    #P_drift = np.exp(log_drift)
    #P_mut = scipy.linalg.expm(R)
    # Define the composite per-generation transition matrix.
    #P = np.dot(P_mut, P_drift)
    # Solve a system of equations to find the stationary distribution.
    #v = MatrixUtil.get_stationary_distribution(P)
    # Try a new thing.
    # The raw drift matrix is scaled so that there are about N*N
    # replacements per generation.
    generation_rate = 1.0
    R_drift_raw = wrightcore.create_moran_drift_rate_k3(M, T)
    R_drift = (generation_rate / float(N)) * R_drift_raw
    #FIXME: you should get the stationary distn directly from the rate matrix
    P = scipy.linalg.expm(R_mut + R_drift)
    v = MatrixUtil.get_stationary_distribution(P)
    """
    for state, value in zip(M, v):
        print state, value
    """
    # draw an equilateral triangle
    #drawtri(M, T, v)
    return M, T, v

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_sample_distn(N, n, mutation_rate, fitness_ratio, h):
    """
    @param N: haploid pop size
    @param n: allele sample size
    @param mutation_rate: mutation rate
    @param fitness_ratio: fitness ratio
    @param h: dominance parameter 0.5 when additive
    @return: a distribution over n+1 mono-/di-morphic sample states
    """
    s = 1.0 - fitness_ratio
    P = np.exp(wfengine.create_diallelic_recessive(N // 2, s, h))
    MatrixUtil.assert_transition_matrix(P)
    # allow mutation out of the fixed states
    P[0, 0] = 1.0 - mutation_rate
    P[0, 1] = mutation_rate
    P[-1, -1] = 1.0 - mutation_rate
    P[-1, -2] = mutation_rate
    MatrixUtil.assert_transition_matrix(P)
    # get the population stationary distribution
    v_large = MatrixUtil.get_stationary_distribution(P)
    # get the allele distribution
    v_small = StatsUtil.subsample_pmf_without_replacement(v_large, n)
    return v_small

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_response_content(fs):
    np.set_printoptions(linewidth=200)
    out = StringIO()
    nsamples = 1
    arr = []
    #
    nsites = 50000
    N = 15*2
    k = 4
    params = (0.002, 1, 1, 0, 0, 0)
    #params = (0.008, 1, 1, 0.5, 1, 1.5)
    mutation, fitnesses = kaizeng.params_to_mutation_fitness(N, params)
    #
    tm = time.time()
    P = kaizeng.get_transition_matrix(N, k, mutation, fitnesses)
    print 'time to construct transition matrix:', time.time() - tm
    #
    tm = time.time()
    v = MatrixUtil.get_stationary_distribution(P)
    print 'time to get stationary distribution:', time.time() - tm
    #
    tm = time.time()
    counts = np.random.multinomial(nsites, v)
    print 'time to sample multinomial counts:', time.time() - tm
    #
    tm = time.time()
    logp = StatsUtil.multinomial_log_pmf(v, counts)
    print 'time to get multinomial log pmf:', time.time() - tm
    #
    for i in range(nsamples):
        counts = np.random.multinomial(nsites, v)
        X0 = np.array(params)
        g = G(N, counts)
        Xopt = optimize.fmin(g, X0)
        arr.append(Xopt)
    print >> out, np.array(arr)
    return out.getvalue()

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
 def test_tricky_distribution(self):
     N_diploid = 5
     k = 4
     mutation, fitness = get_test_mutation_fitness()
     v_tricky = get_stationary_distribution_tricky(
             N_diploid, k, mutation, fitness)
     P = get_transition_matrix(
             N_diploid, k, mutation, fitness)
     v_plain = MatrixUtil.get_stationary_distribution(P)
     print 'initial blocks of the stationary distribution:'
     print v_tricky[:4]
     print v_plain[:4]
     print
     print 'normalized initial blocks of the stationary distribution:'
     print v_tricky[:4] / np.sum(v_tricky[:4])
     print v_plain[:4] / np.sum(v_plain[:4])
     print
     print 'next block of the stationary distribution:'
     print v_tricky[4:4+N_diploid*2-1]
     print v_plain[4:4+N_diploid*2-1]
     print
     print 'normalized next block of the stationary distribution:'
     print v_tricky[4:4+N_diploid*2-1] / np.sum(v_tricky[4:4+N_diploid*2-1])
     print v_plain[4:4+N_diploid*2-1] / np.sum(v_plain[4:4+N_diploid*2-1])
     print
     #print 'subsequent blocks of the stationary distribution:'
     #print v_tricky[4:]
     #print v_plain[4:]
     #print
     #print 'normalized subsequent blocks of the stationary distribution:'
     #print v_tricky[4:] / np.sum(v_tricky[4:])
     #print v_plain[4:] / np.sum(v_plain[4:])
     #print
     #print 'ratio of normalizers:'
     #print np.sum(v_tricky[:4]) / np.sum(v_plain[:4])
     #print
     self.assertTrue(np.allclose(v_tricky, v_plain))

# 需要导入模块: import MatrixUtil [as 别名]
# 或者: from MatrixUtil import get_stationary_distribution [as 别名]
def get_response_content(fs):
    N_diploid = fs.N_diploid
    N = N_diploid * 2
    k = 2
    gamma = fs.gamma
    # define the fitnesses and the selection value
    f0 = 1.0
    f1 = 1.0 - gamma / N
    s = 1 - f1 / f0
    if f1 <= 0:
        raise ValueError('the extreme selection caused a non-positive fitness')
    # get a wright fisher transition matrix
    P = np.exp(wfengine.create_genic_diallelic(N_diploid, s))
    """
    # condition on no fixation
    for i in range(N):
        P[i] /= 1 - P[i, N]
    # remove the fixed state from the transition matrix
    P = P[:N, :N]
    """
    # add mutations
    P[0, 0] = 0
    P[0, 1] = 1
    P[N, N] = 0
    P[N, 1] = 1
    # compute the stationary distribution
    v = MatrixUtil.get_stationary_distribution(P)
    # get the distribution over dimorphic states
    h = v[1:-1]
    h /= np.sum(h)
    # look at continuous approximations
    w = np.zeros(N+1)
    for i in range(1, N):
        x = i / float(N)
        #x0 = i / float(N)
        #x1 = (i + 1) / float(N)
        #value = sojourn_definite(x0, x1, gamma)
        value = sojourn_kernel(x, gamma)
        w[i] = value
    w = w[1:-1]
    w /= np.sum(w)
    # get the array for the R plot
    arr = [h.tolist(), w.tolist()]
    # define the r script
    out = StringIO()
    print >> out, 'title.string <- "allele 1 vs allele 2"'
    print >> out, 'mdat <-', RUtil.matrix_to_R_string(arr)
    print >> out, mk_call_str(
            'barplot',
            'mdat',
            'legend.text=' + mk_call_str(
                'c',
                '"exact discrete distribution"',
                '"continuous approximation"',
                #'"two-allele large N limit"',
                #'"two-allele"',
                #'"four-allele without mutational bias"',
                #'"four-allele with mutational bias (kappa_{1,2}=2)"',
                #'"four-allele with mutational bias, large N limit"',
                ),
            'args.legend = list(x="topright", bty="n")',
            'names.arg = 1:%s' % (N-1),
            main='title.string',
            xlab='"frequency of allele 1"',
            ylab='"frequency"',
            col=mk_call_str(
                'c',
                #'"red"',
                #'"white"',
                '"black"',
                #'"gray"',
                '"red"',
                ),
            beside='TRUE',
            border='NA',
            )
    #print >> out, 'box()'
    script = out.getvalue().rstrip()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter_no_table(
            script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data