Python scipy.stats.poisson.pmf() Examples

def compute_enriched_threshold(average_window_readcount):
    # type: (float) -> int
    """
    Computes the minimum number of tags required in window for an island to be enriched.
    """

    current_threshold, survival_function = 0, 1
    for current_threshold in count(start=0, step=1):
        survival_function -= poisson.pmf(current_threshold,
                                         average_window_readcount)
        if survival_function <= WINDOW_P_VALUE:
            break

    island_enriched_threshold = current_threshold + 1

    return island_enriched_threshold 

def compute_window_score(i, poisson_parameter):
    # type: (int, float) -> float

    # No enrichment; poisson param also average
    if i < poisson_parameter:
        return 0

    p_value = poisson.pmf(i, poisson_parameter)

    if p_value > 0:
        window_score = -log(p_value)
    else:
        # log of zero not defined
        window_score = 1000

    return window_score 

def test_update_with_noise_coherent_value_error():
    """Tests the correct error is raised"""
    cutoff = 15
    num_modes = 3
    nbar_vals = np.random.rand(num_modes - 1)
    noise_dists = np.array([poisson.pmf(np.arange(cutoff), nbar) for nbar in nbar_vals])
    hbar = 2
    beta = np.random.rand(num_modes) + 1j * np.random.rand(num_modes)
    means = Means(np.concatenate((beta, beta.conj())), hbar=hbar)
    cov = hbar * np.identity(2 * num_modes) / 2
    cutoff = 10
    probs = probabilities(means, cov, cutoff, hbar=2)
    with pytest.raises(
        ValueError,
        match="The list of probability distributions probs_noise and the tensor of probabilities probs have incompatible dimensions.",
    ):
        update_probabilities_with_noise(noise_dists, probs) 

def test_update_with_noise_coherent(num_modes, parallel, monkeypatch):
    """ Test that adding noise on coherent states gives the same probabilities at some other coherent states"""

    if parallel: # set single-thread use in OpenMP
        monkeypatch.setenv("OMP_NUM_THREADS", "1")

    cutoff = 15
    nbar_vals = np.random.rand(num_modes)
    noise_dists = np.array([poisson.pmf(np.arange(cutoff), nbar) for nbar in nbar_vals])
    hbar = 2
    beta = np.random.rand(num_modes) + 1j * np.random.rand(num_modes)
    means = Means(np.concatenate((beta, beta.conj())), hbar=hbar)
    cov = hbar * np.identity(2 * num_modes) / 2
    cutoff = 10

    probs = probabilities(means, cov, cutoff, parallel=parallel, hbar=2)
    updated_probs = update_probabilities_with_noise(noise_dists, probs)
    beta_expected = np.sqrt(nbar_vals + np.abs(beta) ** 2)
    means_expected = Means(
        np.concatenate((beta_expected, beta_expected.conj())), hbar=hbar
    )
    expected = probabilities(means_expected, cov, cutoff, parallel=parallel, hbar=2)
    assert np.allclose(updated_probs, expected) 

def test_pmf_p1(self):
        poisson_pmf = poisson.pmf(1, 1)
        genpoisson_pmf = sm.distributions.genpoisson_p.pmf(1, 1, 0, 1)
        assert_allclose(poisson_pmf, genpoisson_pmf, rtol=1e-15) 

def level_distribution(self, k, mu):
        _dists = np.array([poisson.pmf(kk, mu) for kk in range(1, k+1)])
        return _dists / np.sum(_dists) 

def single_gap_factor(island_enriched_threshold,
                      poisson_distribution_parameter):
    # type: (int, float) -> float

    poisson_scores = [poisson.pmf(i, poisson_distribution_parameter)
                      for i in range(island_enriched_threshold)]
    return sum(poisson_scores) 

def level_distribution(self, k, mu):
        _dists = np.array([poisson.pmf(kk, mu) for kk in range(1, k + 1)])
        return _dists / np.sum(_dists) 

def test_pmf(self):
        n, p = sm.distributions.zinegbin.convert_params(1, 0.9, 1)
        nb_logpmf = nbinom.pmf(2, n, p)
        tnb_pmf = sm.distributions.zinegbin.pmf(2, 1, 0.9, 2, 0.5)
        assert_allclose(nb_logpmf, tnb_pmf * 2, rtol=1e-7) 

def test_pmf_p2(self):
        n, p = sm.distributions.zinegbin.convert_params(30, 0.1, 2)
        nb_pmf = nbinom.pmf(100, n, p)
        tnb_pmf = sm.distributions.zinegbin.pmf(100, 30, 0.1, 2, 0.01)
        assert_allclose(nb_pmf, tnb_pmf, rtol=1e-5, atol=1e-5) 

def test_pmf(self):
        gp_pmf = sm.distributions.genpoisson_p.pmf(3, 2, 2, 2)
        zigp_pmf = sm.distributions.zigenpoisson.pmf(3, 2, 2, 2, 0.1)
        assert_allclose(gp_pmf, zigp_pmf, rtol=5e-2, atol=5e-2) 

def test_pmf(self):
        poisson_pmf = poisson.pmf(2, 2)
        zipoisson_pmf = sm.distributions.zipoisson.pmf(2, 2, 0.1)
        assert_allclose(poisson_pmf, zipoisson_pmf, rtol=5e-2, atol=5e-2) 

def test_pmf_zero(self):
        poisson_pmf = poisson.pmf(3, 2)
        zipoisson_pmf = sm.distributions.zipoisson.pmf(3, 2, 0)
        assert_allclose(poisson_pmf, zipoisson_pmf, rtol=1e-12) 

def test_pmf_p5(self):
        poisson_pmf = poisson.pmf(10, 2)
        genpoisson_pmf_5 = sm.distributions.genpoisson_p.pmf(10, 2, 1e-25, 5)
        assert_allclose(poisson_pmf, genpoisson_pmf_5, rtol=1e-12) 

def test_pmf_p2(self):
        poisson_pmf = poisson.pmf(2, 2)
        genpoisson_pmf = sm.distributions.genpoisson_p.pmf(2, 2, 0, 2)
        assert_allclose(poisson_pmf, genpoisson_pmf, rtol=1e-15) 

def poisson_probability(n, lam):
    global poisson_cache
    key = n * 10 + lam
    if key not in poisson_cache:
        poisson_cache[key] = poisson.pmf(n, lam)
    return poisson_cache[key] 

def defineOptimization():
    
    model = AbstractModel()

    # Set of orders
    model.S = Set()

    # Budget 
    model.V = Param(within=NonNegativeReals)

    # Purchase cost
    model.c = Param(model.S, within=NonNegativeReals)

    # Holding/refurbishing cost
    model.h = Param(model.S, within=NonNegativeReals)

    # Expected demand
    model.demand = Param(model.S, within=NonNegativeReals)

    # Cost of the lost sale due to the shortage of supply (economic cost)
    model.b = Param(model.S, within=NonNegativeReals)

    # Number of items of each product ordered
    model.s = Var(model.S, within=NonNegativeIntegers, initialize=0)

    # Minimize ordering cost
    def cost_rule(model):

        return sum(model.h[i] * (model.s[i] - model.demand[i]) + 
                   (model.h[i] + model.b[i]) * sum((x - model.s[i]) * poisson.pmf(x,model.demand[i]) * (1/(1+exp(-10*(log(x+1)-log(model.s[i]+1))))) 
                                                   for x in sequence(0, int(model.V / model.c[i]))) 
                                               # we have to sum from 0 and not from s[i] since the indices of sum are determined 
                                               # at compilation time
                                               # 1/(1+exp(...)) term is an approximation of the step function 
                                               # f(x,s[i]) = 1 if x>s[i], f(x,s[i])=0 if x<=s[i]
                   for i in model.S)

    model.cost = Objective(rule=cost_rule)

    # Budget constraint
    def budget_rule(model):
        return summation(model.c, model.s) <= model.V

    model.budget_constraint = Constraint(rule=budget_rule)

    return model 

def determine_thresholds(kernels, ledges, gw_hist, fdr):
    """
    given a 'gw_hist' histogram of observed counts
    for each lambda-chunk for each kernel-type, and
    also given a FDR, calculate q-values for each observed
    count value in each lambda-chunk for each kernel-type.

    Returns
    -------
    threshold_df : dict
      each threshold_df[k] is a Series indexed by la_exp intervals
      (IntervalIndex) and it is all we need to extract "good" pixels from
      each chunk ...
    qvalues : dict
      A dictionary with keys being kernel names and values pandas.DataFrames
      storing q-values: each column corresponds to a lambda-chunk,
      while rows correspond to observed pixels values.


    """
    rcs_hist = {}
    rcs_Poisson = {}
    qvalues = {}
    threshold_df = {}
    for k in kernels:
        # Reverse cumulative histogram for this kernel.
        # First row contains total # of pixels in each lambda-chunk.
        rcs_hist[k] = gw_hist[k].iloc[::-1].cumsum(axis=0).iloc[::-1]

        # Assign a unit Poisson distribution to each lambda-chunk.
        # The expected value is the upper boundary of the lambda-chunk.
        #   poisson.sf = 1 - poisson.cdf, but more precise
        #   poisson.sf(-1,mu) == 1.0, i.e. is equivalent to the
        #   poisson.pmf(gw_hist[k].index, mu)[::-1].cumsum()[::-1]
        rcs_Poisson[k] = pd.DataFrame()
        for mu, column in zip(ledges[1:-1], gw_hist[k].columns):
            renorm_factors = rcs_hist[k].loc[0, column]
            rcs_Poisson[k][column] = renorm_factors * poisson.sf(
                gw_hist[k].index - 1, mu
            )

        # Determine the threshold by checking the value at which 'fdr_diff'
        # first turns positive. Fill NaNs with an "unreachably" high value.
        fdr_diff = fdr * rcs_hist[k] - rcs_Poisson[k]
        very_high_value = len(rcs_hist[k])
        threshold_df[k] = (
            fdr_diff.where(fdr_diff > 0)
            .apply(lambda col: col.first_valid_index())
            .fillna(very_high_value)
            .astype(np.integer)
        )
        # q-values
        # bear in mind some issues with lots of NaNs and Infs after
        # such a brave operation ...
        qvalues[k] = rcs_Poisson[k] / rcs_hist[k]

    return threshold_df, qvalues 

def compute_score_threshold(average_window_readcount,
                            island_enriched_threshold,
                            gap_contribution, boundary_contribution,
                            genome_length_in_bins):
    # type: (float, int, float, float, float) -> float
    """
    What does island_expectations do?
    """

    required_p_value = poisson.pmf(island_enriched_threshold,
                                   average_window_readcount)

    prob = boundary_contribution * required_p_value

    score = -log(required_p_value)

    current_scaled_score = int(round(score / BIN_SIZE))

    island_expectations_d = {}  # type: Dict[int, float]
    island_expectations_d[current_scaled_score] = prob * genome_length_in_bins
    island_expectations_d[
        0] = boundary_contribution * genome_length_in_bins / gap_contribution

    current_max_scaled_score = current_scaled_score

    interval = int(1 / BIN_SIZE)
    partial_cumu = 0.0
    logging.info("Finding the score required to consider an island enriched.")
    while (partial_cumu > E_VALUE_THRESHOLD or partial_cumu < 1e-100):

        current_scaled_score += interval
        current_max_scaled_score = current_scaled_score - interval
        # logging.debug(island_expectations_d)

        if current_scaled_score > current_max_scaled_score:

            # logging.debug(island_expectations_d)
            island_expectations_d = add_to_island_expectations_dict(
                average_window_readcount, current_max_scaled_score,
                island_enriched_threshold, island_expectations_d,
                gap_contribution)
            partial_cumu = 0.0001
            current_max_scaled_score += 1000

            if max(island_expectations_d) > interval:
                partial_cumu = sum(
                    [val
                     for idx, val in island_expectations_d.items()
                     if idx > current_max_scaled_score - interval])
            else:
                partial_cumu = sum(island_expectations_d.values())

    logging.debug("Computing cumulative distribution.")
    score_threshold = generate_cumulative_dist(island_expectations_d,
                                               current_max_scaled_score + 1)
    logging.info("Enriched score threshold for islands: " + str(
        score_threshold))
    return score_threshold