Python binom.pmf方法代码示例

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def test_R(self):
        # test against the values produced by this R code
        # (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Multinom.html)
        # X <- t(as.matrix(expand.grid(0:3, 0:3))); X <- X[, colSums(X) <= 3]
        # X <- rbind(X, 3:3 - colSums(X)); dimnames(X) <- list(letters[1:3], NULL)
        # X
        # apply(X, 2, function(x) dmultinom(x, prob = c(1,2,5)))

        n, p = 3, [1./8, 2./8, 5./8]
        r_vals = {(0, 0, 3): 0.244140625, (1, 0, 2): 0.146484375,
                  (2, 0, 1): 0.029296875, (3, 0, 0): 0.001953125,
                  (0, 1, 2): 0.292968750, (1, 1, 1): 0.117187500,
                  (2, 1, 0): 0.011718750, (0, 2, 1): 0.117187500,
                  (1, 2, 0): 0.023437500, (0, 3, 0): 0.015625000}
        for x in r_vals:
            assert_allclose(multinomial.pmf(x, n, p), r_vals[x], atol=1e-14) 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def test_pmf(self):
        vals0 = multinomial.pmf((5,), 5, (1,))
        assert_allclose(vals0, 1, rtol=1e-8)

        vals1 = multinomial.pmf((3,4), 7, (.3, .7))
        assert_allclose(vals1, .22689449999999994, rtol=1e-8)

        vals2 = multinomial.pmf([[[3,5],[0,8]], [[-1, 9], [1, 1]]], 8,
                (.1, .9))
        assert_allclose(vals2, [[.03306744, .43046721], [0, 0]], rtol=1e-8)

        x = np.empty((0,2), dtype=np.float64)
        vals3 = multinomial.pmf(x, 4, (.3, .7))
        assert_equal(vals3, np.empty([], dtype=np.float64))

        vals4 = multinomial.pmf([1,2], 4, (.3, .7))
        assert_allclose(vals4, 0, rtol=1e-8)
        
        vals5 = multinomial.pmf([3, 3, 0], 6, [2/3.0, 1/3.0, 0])
        assert_allclose(vals5, 0.219478737997, rtol=1e-8) 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def pdf_integral(p1, data):

    # calculate pdens for range of p1
    xj, nj, c, p2, var = data

    dens = pdf(p1, data=data[2:])

    return dens * binom.pmf(xj, nj, p=p1)


# This is recoded from the implementation of xpclr. I do not think it represents
# what is discussed in the paper. Here we take the integral of eq 5 in the paper
# then take the ratio of it to eq 5 without the binomial component.
# additionally they neglect the probability of p1 being 0 or 1, I presume to
# allow the romberg integration to converge ok.
# This calculates the likelihood of a given SNP 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def _mid_p_value(self):
        r"""
        Calculates the mid-p-value of the McNemar test.

        Returns
        -------
        mid_p : float
            The mid-p value.

        Notes
        -----
        The McNemar mid-p test is calculated by subtracting half the point probability of the observed :math:`n_{12}`
        cell of the contingency table from the one-sided :math:`p_{exact}` value using the equation above. The
        resulting p-value is then doubled to obtain the two-sided mid-p-value. Stated more formally, the McNemar
        mid-p test is defined as:

        .. math::

            p_{mid} = 2 \large(\sum^n_{i=b} \binom{n}{i} \frac{1}{2}^i (1 - \frac{1}{2})^{n - i} -
            \frac{1}{2} \binom{n}{b} \frac{1}{2}^b (1 - \frac{1}{2}^{n -b } \large)

        """
        n = self.table[1, 0] + self.table[0, 1]
        mid_p = self.exact_p_value - binom.pmf(self.table[0, 1], n, 0.5)

        return mid_p 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def log_likelihood_binom(E, nd, na):
    """Likelihood function for (nd,na) to be from a binom with p=E (no BG)."""
    return -(np.log(binom.pmf(na, na+nd, E))).sum() 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def test_reduces_binomial(self):
        # test that the multinomial pmf reduces to the binomial pmf in the 2d
        # case
        val1 = multinomial.logpmf((3, 4), 7, (0.3, 0.7))
        val2 = binom.logpmf(3, 7, 0.3)
        assert_allclose(val1, val2, rtol=1e-8)

        val1 = multinomial.pmf((6, 8), 14, (0.1, 0.9))
        val2 = binom.pmf(6, 14, 0.1)
        assert_allclose(val1, val2, rtol=1e-8) 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def test_pmf_broadcasting(self):
        vals0 = multinomial.pmf([1, 2], 3, [[.1, .9], [.2, .8]])
        assert_allclose(vals0, [.243, .384], rtol=1e-8)

        vals1 = multinomial.pmf([1, 2], [3, 4], [.1, .9])
        assert_allclose(vals1, [.243, 0], rtol=1e-8)

        vals2 = multinomial.pmf([[[1, 2], [1, 1]]], 3, [.1, .9])
        assert_allclose(vals2, [[.243, 0]], rtol=1e-8)

        vals3 = multinomial.pmf([1, 2], [[[3], [4]]], [.1, .9])
        assert_allclose(vals3, [[[.243], [0]]], rtol=1e-8)

        vals4 = multinomial.pmf([[1, 2], [1,1]], [[[[3]]]], [.1, .9])
        assert_allclose(vals4, [[[[.243, 0]]]], rtol=1e-8) 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def test_frozen(self):
        # The frozen distribution should agree with the regular one
        np.random.seed(1234)
        n = 12
        pvals = (.1, .2, .3, .4)
        x = [[0,0,0,12],[0,0,1,11],[0,1,1,10],[1,1,1,9],[1,1,2,8]]
        x = np.asarray(x, dtype=np.float64)
        mn_frozen = multinomial(n, pvals)
        assert_allclose(mn_frozen.pmf(x), multinomial.pmf(x, n, pvals))
        assert_allclose(mn_frozen.logpmf(x), multinomial.logpmf(x, n, pvals))
        assert_allclose(mn_frozen.entropy(), multinomial.entropy(n, pvals)) 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def plot_grid_simulation(self, threshold, threshold_type, n_simulations=100, correction=None):
        '''Create a plot of the simulations'''
        if not self.isfit:
            self.fit()
            self.threshold_simulation(threshold=threshold, threshold_type=threshold_type, correction=correction)
        self.run_multiple_simulations(threshold=threshold, threshold_type=threshold_type, n_simulations=n_simulations)

        if self.signal_mask is None:
            f,a = plt.subplots(ncols=3, figsize=(15, 5))
        else:
            f,a = plt.subplots(ncols=4, figsize=(18, 5))
            a[3].hist(self.multiple_tp)
            a[3].set_ylabel('Frequency', fontsize=18)
            a[3].set_xlabel('Percent Signal Recovery', fontsize=18)
            a[3].set_title('Average Signal Recovery', fontsize=18)

        a[0].imshow(self.t_values)
        a[0].set_title('Random Noise', fontsize=18)
        a[0].axes.get_xaxis().set_visible(False)
        a[0].axes.get_yaxis().set_visible(False)
        a[1].imshow(self.thresholded)
        a[1].set_title(f'Threshold: {threshold_type} = {threshold}', fontsize=18)
        a[1].axes.get_xaxis().set_visible(False)
        a[1].axes.get_yaxis().set_visible(False)
        a[2].plot(binom.pmf(np.arange(0, n_simulations, 1), n_simulations, np.mean(self.multiple_fp>0)))
        a[2].axvline(x=np.mean(self.fpr) * n_simulations, color='r', linestyle='dashed', linewidth=2)
        a[2].set_title(f'False Positive Rate = {self.fpr:.2f}', fontsize=18)
        a[2].set_ylabel('Probability', fontsize=18)
        a[2].set_xlabel('False Positive Rate', fontsize=18)
        plt.tight_layout() 

# 需要导入模块: from scipy.stats import binom [as 别名]
# 或者: from scipy.stats.binom import pmf [as 别名]
def calculate_genotype_likelihoods(reference_support, alternate_support, chromosomes, sample_sexes, homozygous_probability, heterozygous_probability):
    total_reads = reference_support + alternate_support
    gl_hom_ref = binom.pmf(n=total_reads, k=reference_support, p=homozygous_probability)
    gl_het = binom.pmf(n=total_reads, k=reference_support, p=heterozygous_probability)
    gl_hom_alt = binom.pmf(n=total_reads, k=reference_support, p=1 - homozygous_probability)

    # Set probabilities for heterozygous and homozygous calls in hemizygous
    # chromosomes to 0 depending on the sex of the sample.
    gl_het[np.array((chromosomes == "chrX") & (sample_sexes == "M"), dtype=bool)] = 0
    gl_het[np.array((chromosomes == "chrY"), dtype=bool)] = 0
    gl_hom_ref[np.array((chromosomes == "chrY") & (sample_sexes == "F"), dtype=bool)] = 0
    gl_hom_alt[np.array((chromosomes == "chrY") & (sample_sexes == "F"), dtype=bool)] = 0

    return np.array((gl_hom_ref, gl_het, gl_hom_alt)).T