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)

def mc2mnc(mc_):
    '''convert central to non-central moments, uses recursive formula
    optionally adjusts first moment to return mean

    '''
    n = len(mc_)
    mean = mc_[0]
    mc = [1] + list(mc_)    # add zero moment = 1
    mc[1] = 0  # define central mean as zero for formula
    mnc = [1, mean] # zero and first raw moments
    for nn,m in enumerate(mc[2:]):
        n=nn+2
        mnc.append(0)
        for k in range(n+1):
            mnc[n] += scipy.comb(n,k,exact=1) * mc[k] * mean**(n-k)

    return mnc[1:]

def mnc2mc(mnc_, wmean = True):
    '''convert non-central to central moments, uses recursive formula
    optionally adjusts first moment to return mean

    '''
    n = len(mnc_)
    mean = mnc_[0]
    mnc = [1] + list(mnc_)    # add zero moment = 1
    mu = [] #np.zeros(n+1)
    for n,m in enumerate(mnc):
        mu.append(0)
        #[scipy.comb(n-1,k,exact=1) for k in range(n)]
        for k in range(n+1):
            mu[n] += (-1)**(n-k) * scipy.comb(n,k,exact=1) * mnc[k] * mean**(n-k)
    if wmean:
        mu[1] = mean
    return mu[1:]

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

    References
    ----------
    Kenneth Lange: Numerical Analysis for Statisticians, page 40
    (http://books.google.ca/books?id=gm7kwttyRT0C&pg=PA40&lpg=PA40&dq=convert+cumulants+to+moments&source=web&ots=qyIaY6oaWH&sig=cShTDWl-YrWAzV7NlcMTRQV6y0A&hl=en&sa=X&oi=book_result&resnum=1&ct=result)


    '''
    mc = [1,0.0] #kappa_[0]]  #insert 0-moment and mean
    kappa = [1] + list(kappa_)
    for nn,m in enumerate(kappa[2:]):
        n = nn+2
        mc.append(0)
        for k in range(n-1):
            mc[n] += scipy.comb(n-1,k,exact=1) * kappa[n-k]*mc[k]

    mc[1] = kappa_[0] # insert mean as first moments by convention
    return mc[1:]

def diallelic_approximation_d(N_small, g, m0, m1):
    """
    This is experimental.
    The numerical integration should be replaced
    by a call to the confluent hypergeometric function hyp1f1.
    See also
    http://functions.wolfram.com/HypergeometricFunctions/
    Hypergeometric1F1/03/01/04/01/ .
    Also
    www.cs.unc.edu/Research/Image/MIDAG/p01/biostat/Digital_1.pdf .
    Also a gsl implementation gsl_sf_hyperg_1F1_int in hyperg_1F1.c
    specifically hyperg_1F1_ab_posint for positive integers a and b.
    Also
    http://mathworld.wolfram.com/
    ConfluentHypergeometricFunctionoftheFirstKind.html
    """
    hist = np.zeros(N_small + 1)
    for n0 in range(1, N_small):
        n1 = N_small - n0
        prefix = scipy.comb(n0+n1, n0) * special.beta(n0, n1)
        hist[n0] += m0 * prefix * diallelic_d_helper(n0, n1, g)
        hist[n0] += m1 * prefix * diallelic_d_helper(n1, n0, -g)
    return hist[1:-1] / np.sum(hist[1:-1])

def Adams_Moulton(k):
    r""" 
        Construct the k-step, Adams-Moulton method.
        The methods are implicit and have order k+1.
        They have the form:

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

        They are generated using equation (1.9) and the equation in 
        Exercise 3 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]_
    """
    from scipy import comb
    alpha=np.zeros(k+1)
    beta=np.zeros(k+1)
    alpha[k]=1.
    alpha[k-1]=-1.
    gamma=np.zeros(k+1)
    gamma[0]=1.
    beta[k]=1.
    betaj=np.zeros(k+1)
    for j in range(1,k+1):
        gamma[j]= -sum(gamma[:j]/np.arange(j+1,1,-1))
        for i in range(0,j+1):
            betaj[k-i]=(-1.)**i*comb(j,i)*gamma[j]
        beta=beta+betaj
    name=str(k)+'-step Adams-Moulton method'
    return LinearMultistepMethod(alpha,beta,name=name)