Python scipy.comb函数代码示例

    def pdf(self, x, k, n, p):
        '''distribution of success runs of length k or more

        Parameters
        ----------
        x : float
            count of runs of length n
        k : int
            length of runs
        n : int
            total number of observations or trials
        p : float
            probability of success in each Bernoulli trial

        Returns
        -------
        pdf : float
            probability that x runs of length of k are observed

        Notes
        -----
        not yet vectorized

        References
        ----------
        Muselli 1996, theorem 3
        '''

        q = 1-p
        m = np.arange(x, (n+1)//(k+1)+1)[:,None]
        terms = (-1)**(m-x) * comb(m, x) * p**(m*k) * q**(m-1) \
                * (comb(n - m*k, m - 1) + q * comb(n - m*k, m))
        return terms.sum(0)

def mv_hypergeometric(x,m):
    """
    x : number of draws for each category.
    m : size of each category.
    """
    x = np.asarray(x)
    m = np.asarray(m)
    return log(comb(m,x).prod()/comb(m.sum(), x.sum()))

 def runs_prob_odd(self, r):
     n0, n1 = self.n0, self.n1
     k = (r+1)//2
     tmp0 = comb(n0-1, k-1)
     tmp1 = comb(n1-1, k-2)
     tmp3 = comb(n0-1, k-2)
     tmp4 = comb(n1-1, k-1)
     return (tmp0 * tmp1 + tmp3 * tmp4)  / self.comball

def hypergProb(k, N, m, n):
    """ Wikipedia: There is a shipment of N objects in which m are defective. The hypergeometric distribution describes the probability that in a sample of n distinctive objects drawn from the shipment exactly k objects are defective. """
    #return float(choose(m, k) * choose(N-m, n-k)) / choose(N, n)
    hp = float(scipy.comb(m, k) * scipy.comb(N-m, n-k)) / scipy.comb(N, n)
    if scipy.isnan(hp):
        stderr.write("error: not possible to calculate hyperg probability in util.py for k=%d, N=%d, m=%d, n=%d\n" %(k, N, m,n))
        stdout.write("error: not possible to calculate hyperg probability in util.py for k=%d, N=%d, m=%d, n=%d\n" %(k, N, m,n))
    return hp

def dumb_factor(goal, n):
  # Assumes factors are not equal to each other and bounded above by
  # upper_bound.

  comb0 = comb(n, 2)
  for i in range(0, n):
    comb1 = comb(n - i, 2)
    for j in range(i + 1, n):
      if goal == round(comb0 - comb1 + (j - i - 1)):
        return i, j

def hypergeometric(x, n, m, N):
    """
    x : number of successes drawn
    n : number of draws
    m : number of successes in total
    N : successes + failures in total.
    """
    if x < max(0, n-(N-m)):
        return 0.
    elif x > min(n, m):
        return 0.
    else:
        return comb(N-m, x) * comb(m, n-x) / comb(N,n)

    def _calculate_orders(self):
        k = self.k
        n = self.n
        m = self.m
        dim = self.dim
        
        # Calculate the length of each order
        self.order_idx       = np.zeros(n+2, dtype=int) 
        self.order_length    = np.zeros(n+1, dtype=int)
        self.row_counter     = 0

        for ordi in xrange(n+1):    
            self.order_length[ordi] = (sp.comb(n, ordi+1, exact=1) * 
                                        ((m-1)**(ordi+1)))
            self.order_idx[ordi] = self.row_counter
            self.row_counter += self.order_length[ordi]

        self.order_idx[n+1] = dim+1

        # Calculate nnz for A
        # not needed for lil sparse format
        x = (m*np.ones(n))**np.arange(n-1,-1,-1)
        x = x[:k]
        y = self.order_length[:k]
        self.Annz = np.sum(x*y.T)

def triple(g, P, T, r, family):
	if g.intersection_matrix[g.schubert_list.index(P)][g.schubert_list.index(T)] == 0:
		return 0
	else:
		delta = 1
		quad, lin = num_equations(g,P, T)
		subspace_eqns = g.m + r - 1
		#if g.type == 'D' and quad == 0 and r == g.k:
		#	delta = int(family + h(g,P,T))%2
		#if g.type == 'B':
		#	if quad > 0:
		#		quad -=1
		#		if r > g.k:
		#			subspace_eqns +=1
		if g.OG:
			if r > g.k:
				subspace_eqns += 1
				if quad > 0:					
					quad -= 1
			if r == g.k and g.type == 'D':
				if quad == 0:
					delta = int(family + h(g,P,T))%2
				if quad > 0:
					subspace_eqns +=1
					quad -=1
		#subspace_eqns += (1-delta)
		triple_list = []
		for j in range(int(g.N - quad - lin - subspace_eqns)):
			triple_list.append((-1)**j * 2**(quad - j) * sp.comb(quad,j))
		return delta*sum(triple_list)

def backward_difference_formula(k):
    r""" 
        Construct the k-step backward differentiation method.
        The methods are implicit and have order k.
        They have the form:

        `\sum_{j=0}^{k} \alpha_j y_{n+k-j+1} = h \beta_j f(y_{n+1})`

        They are generated using equation (1.22') from Hairer & Wanner III.1,
        along with the binomial expansion.

        .. note::
            Accuracy is lost when evaluating the order conditions
            for methods with many steps.  This could be avoided by using 
            SAGE rationals instead of NumPy doubles for the coefficient
            representation.

        **References**:
            #.[hairer1993]_ pp. 364-365
    """
    from scipy import comb
    alpha=np.zeros(k+1)
    beta=np.zeros(k+1)
    beta[k]=1.
    gamma=np.zeros(k+1)
    gamma[0]=1.
    alphaj=np.zeros(k+1)
    for j in range(1,k+1):
        gamma[j]= 1./j
        for i in range(0,j+1):
            alphaj[k-i]=(-1.)**i*comb(j,i)*gamma[j]
        alpha=alpha+alphaj
    name=str(k)+'-step BDF method'
    return LinearMultistepMethod(alpha,beta,name=name)

def triple(g, P, T, r, family):
	#if not g.index_set_leq(T,P):
	if g.intersection_matrix[g.schubert_list.index(P)][g.schubert_list.index(T)] == 0:
		return 0
	else:
		delta = 1
		quad, lin = num_equations(g,P, T)
		subspace_eqns = g.m + r - 1
		if g.OG and r > g.k:
			subspace_eqns += 1
			if quad > 0:
				quad -= 1
		if g.type == 'D' and r == g.k:
			subspace_eqns = g.n+1
			#if quad == 0 or single_ruling(g,P,T):
			if quad == 0:
				#subspace_eqns -= int((family + h(g,P,T)))%2
				delta = int((family + h(g,P,T)))%2
				subspace_eqns = g.n
			if quad > 0:
				quad -= 1
		triple_list = []
		for j in range(int(g.N - quad - lin - subspace_eqns)):
			triple_list.append((-1)**j * 2**(quad - j) * sp.comb(quad,j))
		return delta*sum(triple_list)

	def triple(self, P, T, r, family):
		#if not self.index_set_leq(T,P):
		if self.intersection_matrix[self.schubert_list.index(P)][self.schubert_list.index(T)] == 0:
			return 0
		else:
			delta = 1
			quad, lin = self.outsourced_num_equations(P, T)
			subspace_eqns = self.m + r - 1
			if self.OG and r > self.k:
				subspace_eqns += 1
				if quad > 0:
					quad -= 1
			if self.type == 'D' and r == self.k:
				subspace_eqns = self.n+1
				#if quad == 0 or self.single_ruling(P,T):
				if quad == 0:
					#subspace_eqns -= int((family + self.h(P,T)))%2
					delta = int((family + self.h(P,T)))%2
					subspace_eqns = self.n
				if quad > 0:
					quad -= 1
			triple_list = []
			for j in range(int(self.N - quad - lin - subspace_eqns)):
				triple_list.append((-1)**j * 2**(quad - j) * sp.comb(quad,j))
			return delta*sum(triple_list)

def similarity_matrix(db, sim_func, check_func):
    matrix = numpy.zeros((len(db["genomes"]), len(db["genomes"])))

    sys.stderr.write("populating matrix\n")
    num_work_to_do = int(scipy.comb(len(db["genomes"]), 2))
    work_done = 0

    for nameA, nameB in itertools.combinations(db["genomes"], 2):
        nameA_num = db["genomes"].index(nameA)
        nameB_num = db["genomes"].index(nameB)
        
        sizeA = db["cds_counts"][nameA_num]
        sizeB = db["cds_counts"][nameB_num]

        hits_count = check_func((db["hit_matrix"][nameA_num][nameB_num],
                                 db["hit_matrix"][nameB_num][nameA_num]))

        matrix[nameA_num][nameB_num] = sim_func(sizeA, sizeB, hits_count)
        matrix[nameB_num][nameA_num] = sim_func(sizeB, sizeA, hits_count)
        
        work_done += 1
        
        if work_done % 100 == 0:
            sys.stderr.write("\r  %s/%s completed" % (work_done, num_work_to_do))

    sys.stderr.write("\r  %s/%s completed\n" % (work_done, num_work_to_do))

    return matrix

def beta_binomial(k, n, a, b):
    """The pmf/pdf of the Beta-binomial distribution.

    Computation based on beta function.

    See: http://en.wikipedia.org/wiki/Beta-binomial_distribution
    and http://mathworld.wolfram.com/BetaBinomialDistribution.html
    
    k = a vector of non-negative integers <= n
    n = an integer
    a = an array of non-negative real numbers
    b = an array of non-negative real numbers
    """
    return (comb(n, k) * beta(k+a, n-k+b) / beta(a,b)).prod(0)

def mnc2cum(mnc_):
    '''convert non-central moments to cumulants
    recursive formula produces as many cumulants as moments

    http://en.wikipedia.org/wiki/Cumulant#Cumulants_and_moments
    '''
    mnc = [1] + list(mnc_)
    kappa = [1]
    for nn,m in enumerate(mnc[1:]):
        n = nn+1
        kappa.append(m)
        for k in range(1,n):
            kappa[n] -= scipy.comb(n-1,k-1,exact=1) * kappa[k]*mnc[n-k]

    return kappa[1:]

def Test():
    # Double check stirlings approximation implementation
    test_vals = [(30, 5),
                 (18, 7),
                 (91, 32)]
    for n, k in test_vals:
        print np.log(scipy.comb(n, k))
        print ln_stirling_binomial(n, k)
    
    m = np.matrix([[50, 39],
                   [66, 5]])
    print CalcPValues(m)
    
    prob_m = np.matrix([[0.25, 0.25],
                        [0.25, 0.25]])
    print CalcPValues(m, prob_m)