python类gamma()的实例源码

def test_cheb1_scheme(scheme):
    def integrate_exact(k):
        # \int_-1^1 x^k / sqrt(1 - x^2)
        if k == 0:
            return numpy.pi
        if k % 2 == 1:
            return 0.0
        return numpy.sqrt(numpy.pi) * ((-1)**k + 1) \
            * math.gamma(0.5*(k+1)) / math.gamma(0.5*k) \
            / k

    degree = check_degree_1d(
            lambda poly: quadpy.line_segment.integrate(
                    poly, numpy.array([[-1.0], [1.0]]), scheme
                    ),
            integrate_exact,
            scheme.degree + 1
            )
    assert degree >= scheme.degree
    return

def test_cheb2_scheme(scheme):
    def integrate_exact(k):
        # \int_-1^1 x^k * sqrt(1 - x^2)
        if k == 0:
            return 0.5 * numpy.pi
        if k % 2 == 1:
            return 0.0
        return numpy.sqrt(numpy.pi) * ((-1)**k + 1) \
            * math.gamma(0.5*(k+1)) / math.gamma(0.5*k+2) \
            / 4

    degree = check_degree_1d(
            lambda poly: quadpy.line_segment.integrate(
                    poly, numpy.array([[-1.0], [1.0]]), scheme
                    ),
            integrate_exact,
            scheme.degree + 1
            )
    assert degree >= scheme.degree
    return

def incomplete_gamma_lower(a, x, precision=9):
    # Numerical approximation of the lower incomplete gamma function to the given precision
    # https://en.wikipedia.org/wiki/Incomplete_gamma_function
    # http://mathworld.wolfram.com/IncompleteGammaFunction.html

    max_diff = 10 ** -precision

    def summand(k):
        return ((-x) ** k) / factorial(k) / (a + k)

    xs = x ** a
    sum_ = summand(0)
    prev_value = xs * sum_  # for k=0
    sum_ += summand(1)
    cur_value = xs * sum_  # for k=1
    k = 1
    while abs(cur_value - prev_value) > max_diff:
        k += 1
        prev_value = cur_value
        sum_ += summand(k)
        cur_value = xs * sum_

    return cur_value

def get_cuckoos(nest,best,lb,ub,n,dim):

    # perform Levy flights
    tempnest=numpy.zeros((n,dim))
    tempnest=numpy.array(nest)
    beta=3/2;
    sigma=(math.gamma(1+beta)*math.sin(math.pi*beta/2)/(math.gamma((1+beta)/2)*beta*2**((beta-1)/2)))**(1/beta);

    s=numpy.zeros(dim)
    for j in range (0,n):
        s=nest[j,:]
        u=numpy.random.randn(len(s))*sigma
        v=numpy.random.randn(len(s))
        step=u/abs(v)**(1/beta)

        stepsize=0.01*(step*(s-best))

        s=s+stepsize*numpy.random.randn(len(s))


        tempnest[j,:]=numpy.clip(s, lb, ub)

    return tempnest

def _hypervolume_couboid(self, cuboid):
        """Computes the hypervolume of a single fuzzified cuboid."""

        all_dims = [dim for domain in self._core._domains.values() for dim in domain]
        n = len(all_dims)

        # calculating the factor in front of the sum
        weight_product = 1.0
        for (dom, dom_weight) in self._weights._domain_weights.items():
            for (dim, dim_weight) in self._weights._dimension_weights[dom].items():
                weight_product *= dom_weight * sqrt(dim_weight)
        factor = self._mu / (self._c**n * weight_product)

        # outer sum
        outer_sum = 0.0        
        for i in range(0, n+1):
            # inner sum
            inner_sum = 0.0
            subsets = list(itertools.combinations(all_dims, i))
            for subset in subsets:
                # first product
                first_product = 1.0
                for dim in set(all_dims) - set(subset):
                    dom = filter(lambda (x,y): dim in y, self._core._domains.items())[0][0]
                    w_dom = self._weights._domain_weights[dom]
                    w_dim = self._weights._dimension_weights[dom][dim]
                    b = cuboid._p_max[dim] - cuboid._p_min[dim]
                    first_product *= w_dom * sqrt(w_dim) * b * self._c

                # second product
                second_product = 1.0
                reduced_domain_structure = self._reduce_domains(self._core._domains, subset)
                for (dom, dims) in reduced_domain_structure.items():
                    n_domain = len(dims)
                    second_product *= factorial(n_domain) * (pi ** (n_domain/2.0))/(gamma((n_domain/2.0) + 1))

                inner_sum += first_product * second_product

            outer_sum += inner_sum
        return factor * outer_sum

def integrate_monomial_over_enr2(k):
    if numpy.any(k % 2 == 1):
        return 0
    return numpy.prod([math.gamma((kk+1) / 2.0) for kk in k])

def integrate_monomial_over_enr(k):
    if numpy.any(k % 2 == 1):
        return 0
    n = len(k)
    return 2 * math.factorial(sum(k) + n - 1) * numpy.prod([
        math.gamma((kk+1) / 2.0) for kk in k
        ]) / math.gamma((sum(k) + n) / 2)

def plot(scheme, show_axes=True):
    ax = plt.gca()
    plt.axis('equal')

    if not show_axes:
        ax.set_axis_off()

    n = 2
    I0 = 2*math.factorial(n-1)*math.pi**(0.5*n) / math.gamma(0.5*n)

    helpers.plot_disks(
        plt, scheme.points, scheme.weights, I0
        )
    return

def lower_incomplete_gamma2(a,x):
    return gamma(a)-upper_incomplete_gamma2(a,x)

def gammainc(a,x):
    return lower_incomplete_gamma(a,x)/gamma(a)

def gammaincc(a,x):
    return upper_incomplete_gamma(a,x)/gamma(a)

def __init__(self, a, b, K):
        from operator import mul
        self._a = a
        self._b = b
        self._K = K
        theta = 1. / b
        self.gamma_dist = gamma_dist(a, 0, theta)
        # self._alpha = np.array(alpha)
        # self._coef = gamma(np.sum(self._alpha)) / \
        #              reduce(mul, [gamma(a) for a in self._alpha])

def gdir_pdf_integ(self, alpha1, alpha2, alpha3):
        from math import gamma
        from operator import mul
        alpha = np.array([alpha1, alpha2, alpha3])
        coef1 = gamma(np.sum(alpha)) / reduce(mul, [gamma(a) for a in alpha])
        coef2 = reduce(mul, [xx ** (aa - 1)
                             for (xx, aa) in zip(self.x, alpha)])

        coef3 = reduce(mul, [self.gamma_dist.pdf(local) for local in alpha])
        return coef1 * coef2 * coef3

def __init__(self, alpha):
        from math import gamma
        from operator import mul
        self._alpha = np.array(alpha)
        self._coef = gamma(np.sum(self._alpha)) / \
                     reduce(mul, [gamma(a) for a in self._alpha])

def integrand(x, lam, sigma, H):
    C2 = lambda H: np.pi / (H*gamma(2*H)*np.sin(H*np.pi))
    A = (lam*sigma)**2 / C2(H)
    return A * (2*(1-np.cos(x))*np.abs(x)**(-2*H -1)) / (lam**2 - 2*lam*(1-np.cos(x)) + 2*(1-np.cos(x)))

##################################################################################

def integrand(x, lam, sigma, H):
    C2 = lambda H: np.pi / (H*gamma(2*H)*np.sin(H*np.pi))
    A = (lam*sigma)**2 / C2(H)
    return A * (2*(1-np.cos(x))*np.abs(x)**(-2*H -1)) / (lam**2 - 2*lam*(1-np.cos(x)) + 2*(1-np.cos(x)))

##################################################################################

def sigma(beta):
    p=math.gamma(1+beta)* math.sin(math.pi*beta/2)/(math.gamma((1+beta)/2)*beta*(pow(2,(beta-1)/2)))
    return pow(p,1/beta)

def exchange_binary(binary,score):#,alpha,beta,gamma,r):

    #binary in list
    al_binary=binary
    #movement=move(b,alpha,beta,gamma,r)
    movement=math.tanh(score)
    ##al_binary=[case7(b) if random.uniform(0,1) < movement else case8(b) for b in binary]
    if random.uniform(0,1) < movement:
        for i,b in enumerate(binary):
            al_binary[i]=case7(b)
    else:
        for i,b in enumerate(binary):
            al_binary[i]=case8(b)
    return al_binary

def sigma(beta):
    p=math.gamma(1+beta)* math.sin(math.pi*beta/2)/(math.gamma((1+beta)/2)*beta*(pow(2,(beta-1)/2)))
    return pow(p,1/beta)

def exchange_binary(binary,score):#,alpha,beta,gamma,r):

    #binary in list
    al_binary=binary
    #movement=move(b,alpha,beta,gamma,r)
    movement=math.tanh(score)
    ##al_binary=[case7(b) if random.uniform(0,1) < movement else case8(b) for b in binary]
    if random.uniform(0,1) < movement:
        for i,b in enumerate(binary):
            al_binary[i]=case7(b)
    else:
        for i,b in enumerate(binary):
            al_binary[i]=case8(b)
    return al_binary

def test_compute_moments():
    moments = orthopy.line.compute_moments(lambda x: 1, -1, +1, 5)
    assert (
        moments == [2, 0, sympy.Rational(2, 3), 0, sympy.Rational(2, 5)]
        ).all()

    moments = orthopy.line.compute_moments(
            lambda x: 1, -1, +1, 5,
            polynomial_class=orthopy.line.legendre
            )
    assert (moments == [2, 0, 0, 0, 0]).all()

    # Example from Gautschi's "How to and how not to" article
    moments = orthopy.line.compute_moments(
            lambda x: sympy.exp(-x**3/3),
            0, sympy.oo,
            5
            )

    reference = [
        3**sympy.Rational(n-2, 3) * sympy.gamma(sympy.Rational(n+1, 3))
        for n in range(5)
        ]

    assert numpy.all([
        sympy.simplify(m - r) == 0
        for m, r in zip(moments, reference)
        ])
    return

def B(alpha, beta):
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)

def test_incomplete_gamma_sum(self):
        """Test if the sum of lower and upper incomplete gamma functions correctly produce the
        gamma function"""

        rand = Random(372924)
        for _ in range(10):
            a, x = rand.random() * 10, rand.random() * 10
            l = incomplete_gamma_lower(a, x)
            r = incomplete_gamma_upper(a, x)
            g = gamma(a)
            self.assertAlmostEqual(l + r, g, 10)

def test_incomplete_gamma_lower(self):
        """Test the lower incomplete gamma function approximation against pre-computed values."""

        # Pre-computed values from Octave, generated with
        # for a=1:10 for x=1:10 printf("(%d,%d,%.7f),", a, x, gammainc(x,a)*gamma(a)) endfor endfor
        values_table = (
            (1, 1, 0.6321206), (1, 2, 0.8646647), (1, 3, 0.9502129), (1, 4, 0.9816844),
            (1, 5, 0.9932621), (1, 6, 0.9975212), (1, 7, 0.9990881), (1, 8, 0.9996645),
            (1, 9, 0.9998766), (1, 10, 0.9999546), (2, 1, 0.2642411), (2, 2, 0.5939942),
            (2, 3, 0.8008517), (2, 4, 0.9084218), (2, 5, 0.9595723), (2, 6, 0.9826487),
            (2, 7, 0.9927049), (2, 8, 0.9969808), (2, 9, 0.9987659), (2, 10, 0.9995006),
            (3, 1, 0.1606028), (3, 2, 0.6466472), (3, 3, 1.1536198), (3, 4, 1.5237934),
            (3, 5, 1.7506960), (3, 6, 1.8760624), (3, 7, 1.9407277), (3, 8, 1.9724921),
            (3, 9, 1.9875356), (3, 10, 1.9944612), (4, 1, 0.1139289), (4, 2, 0.8572592),
            (4, 3, 2.1166087), (4, 4, 3.3991793), (4, 5, 4.4098445), (4, 6, 5.0927767),
            (4, 7, 5.5094075), (4, 8, 5.7457193), (4, 9, 5.8726411), (4, 10, 5.9379837),
            (5, 1, 0.0878363), (5, 2, 1.2636724), (5, 3, 4.4336821), (5, 4, 8.9079136),
            (5, 5, 13.4281612), (5, 6, 17.1586440), (5, 7, 19.8482014), (5, 8, 21.6088224),
            (5, 9, 22.6808726), (5, 10, 23.2979355), (6, 1, 0.0713022), (6, 2, 1.9876330),
            (6, 3, 10.0701530), (6, 4, 25.7843536), (6, 5, 46.0847214), (6, 6, 66.5184430),
            (6, 7, 83.9150069), (6, 8, 97.0516726), (6, 9, 106.1171375), (6, 10, 111.9496845),
            (7, 1, 0.0599336), (7, 2, 3.2643400), (7, 3, 24.1261454), (7, 4, 79.6852644),
            (7, 5, 171.2279067), (7, 6, 283.4619967), (7, 7, 396.2080398), (7, 8, 494.3705202),
            (7, 9, 571.1177953), (7, 10, 626.2981770), (8, 1, 0.0516560), (8, 2, 5.5274636),
            (8, 3, 59.9986994), (8, 4, 257.7134236), (8, 5, 672.1932373), (8, 6, 1290.3420073),
            (8, 7, 2022.4822690), (8, 8, 2757.0775202), (8, 9, 3407.5592999), (8, 10, 3930.0879411),
            (9, 1, 0.0453682), (9, 2, 9.5738763), (9, 3, 153.3366399), (9, 4, 861.3736786),
            (9, 5, 2745.5353520), (9, 6, 6159.3842425), (9, 7, 10923.0400848),
            (9, 8, 16428.4911932), (9, 9, 21948.0869937), (9, 10, 26900.7105528),
            (10, 1, 0.0404341), (10, 2, 16.8732215), (10, 3, 400.0708927), (10, 4, 2951.0282662),
            (10, 5, 11549.7654353), (10, 6, 30454.3472871), (10, 7, 61509.6342948),
            (10, 8, 102831.3889931), (10, 9, 149721.2962968), (10, 10, 196706.4652125),
        )

        for a, x, inc_gam_value in values_table:
            self.assertAlmostEqual(incomplete_gamma_lower(a, x), inc_gam_value,
                                   places=min(7, int(6 - log10(inc_gam_value))))

def test_incomplete_gamma_upper(self):
        """Test the upper incomplete gamma function approximation against pre-computed values."""

        # Pre-computed values from Octave, generated with:
        #   for a=1:10 for x=1:10
        #       printf("(%d,%d,%.7f),", a, x, gammainc(x,a,'upper')*gamma(a))
        #   endfor endfor
        values_table = (
            (1, 1, 0.3678794), (1, 2, 0.1353353), (1, 3, 0.0497871), (1, 4, 0.0183156),
            (1, 5, 0.0067379), (1, 6, 0.0024788), (1, 7, 0.0009119), (1, 8, 0.0003355),
            (1, 9, 0.0001234), (1, 10, 0.0000454), (2, 1, 0.7357589), (2, 2, 0.4060058),
            (2, 3, 0.1991483), (2, 4, 0.0915782), (2, 5, 0.0404277), (2, 6, 0.0173513),
            (2, 7, 0.0072951), (2, 8, 0.0030192), (2, 9, 0.0012341), (2, 10, 0.0004994),
            (3, 1, 1.8393972), (3, 2, 1.3533528), (3, 3, 0.8463802), (3, 4, 0.4762066),
            (3, 5, 0.2493040), (3, 6, 0.1239376), (3, 7, 0.0592723), (3, 8, 0.0275079),
            (3, 9, 0.0124644), (3, 10, 0.0055388), (4, 1, 5.8860711), (4, 2, 5.1427408),
            (4, 3, 3.8833913), (4, 4, 2.6008207), (4, 5, 1.5901555), (4, 6, 0.9072233),
            (4, 7, 0.4905925), (4, 8, 0.2542807), (4, 9, 0.1273589), (4, 10, 0.0620163),
            (5, 1, 23.9121637), (5, 2, 22.7363276), (5, 3, 19.5663179), (5, 4, 15.0920864),
            (5, 5, 10.5718388), (5, 6, 6.8413560), (5, 7, 4.1517986), (5, 8, 2.3911776),
            (5, 9, 1.3191274), (5, 10, 0.7020645), (6, 1, 119.9286978), (6, 2, 118.0123670),
            (6, 3, 109.9298470), (6, 4, 94.2156464), (6, 5, 73.9152786), (6, 6, 53.4815570),
            (6, 7, 36.0849931), (6, 8, 22.9483274), (6, 9, 13.8828625), (6, 10, 8.0503155),
            (7, 1, 719.9400664), (7, 2, 716.7356600), (7, 3, 695.8738546), (7, 4, 640.3147356),
            (7, 5, 548.7720933), (7, 6, 436.5380033), (7, 7, 323.7919602), (7, 8, 225.6294798),
            (7, 9, 148.8822047), (7, 10, 93.7018230), (8, 1, 5039.9483440), (8, 2, 5034.4725364),
            (8, 3, 4980.0013006), (8, 4, 4782.2865764), (8, 5, 4367.8067627), (8, 6, 3749.6579927),
            (8, 7, 3017.5177310), (8, 8, 2282.9224798), (8, 9, 1632.4407001), (8, 10, 1109.9120589),
            (9, 1, 40319.9546318), (9, 2, 40310.4261237), (9, 3, 40166.6633601),
            (9, 4, 39458.6263214), (9, 5, 37574.4646480), (9, 6, 34160.6157575),
            (9, 7, 29396.9599152), (9, 8, 23891.5088068), (9, 9, 18371.9130063),
            (9, 10, 13419.2894472), (10, 1, 362879.9595659), (10, 2, 362863.1267785),
            (10, 3, 362479.9291073), (10, 4, 359928.9717338), (10, 5, 351330.2345647),
            (10, 6, 332425.6527129), (10, 7, 301370.3657052), (10, 8, 260048.6110069),
            (10, 9, 213158.7037032), (10, 10, 166173.5347875)
        )

        for a, x, inc_gam_value in values_table:
            self.assertAlmostEqual(incomplete_gamma_upper(a, x), inc_gam_value,
                                   places=min(7, int(6 - log10(inc_gam_value))))

def incomplete_gamma_upper(a, x, precision=9):
    # Numerical approximation of the lower incomplete gamma function to the given precision
    # https://en.wikipedia.org/wiki/Incomplete_gamma_function
    # http://mathworld.wolfram.com/IncompleteGammaFunction.html
    return gamma(a) - incomplete_gamma_lower(a, x, precision)

def incomplete_gamma_upper_norm(a, x, precision=9):
    return incomplete_gamma_upper(a, x, precision) / gamma(a)

def epsilon(data, MinPts):
    '''????
    input:  data(mat):????
            MinPts(int):??????????
    output: eps(float):??
    '''
    m, n = np.shape(data)
    xMax = np.max(data, 0)
    xMin = np.min(data, 0)
    eps = ((np.prod(xMax - xMin) * MinPts * math.gamma(0.5 * n + 1)) / (m * math.sqrt(math.pi ** n))) ** (1.0 / n)
    return eps

def train(X, y):
    # This function is used for training our Bayesian model
    # Returns the regression parameters w_opt, and alpha, beta, S_N
    # needed for the predictive distribution

    Phi = X # the measurement matrix of the input variables x (i.e., features)
    t   = y # the vector of observations for the target variable

    (N, M) = np.shape(Phi)

    # Init values for  hyper-parameters alpha, beta
    alpha = 5*10**(-3)
    beta = 5
    max_iter = 100
    k = 0

    PhiT_Phi = np.dot(np.transpose(Phi), Phi)
    s = np.linalg.svd(PhiT_Phi, compute_uv=0) # Just get the vector of singular values s

    ab_old = np.array([alpha, beta])
    ab_new = np.zeros((1,2))
    tolerance = 10**-3
    while( k < max_iter and np.linalg.norm(ab_old-ab_new) > tolerance):
        k += 1
        try:
            S_N = np.linalg.inv(alpha*np.eye(M) + beta*PhiT_Phi)
        except np.linalg.LinAlgError as err:
            print  "******************************************************************************************************"
            print "                           ALERT: LinearAlgebra Error detected!"
            print "      CHECK if your measurement matrix is not leading to a singular alpha*np.eye(M) + beta*PhiT_Phi"
            print "                           GOODBYE and see you later. Exiting ..."
            print  "******************************************************************************************************"
            sys.exit(-1)

        m_N = beta * np.dot(S_N, np.dot(np.transpose(Phi), t))
        gamma = sum(beta*s[i]**2 /(alpha + beta*s[i]**2) for i in range(M))
        #
        # update alpha, beta
        #
        ab_old = np.array([alpha, beta])
        alpha = gamma /np.inner(m_N,m_N)
        one_over_beta = 1/(N-gamma) * sum( (t[n] - np.inner(m_N, Phi[n]))**2 for n in range(N))
        beta = 1/one_over_beta
        ab_new = np.array([alpha, beta])

    S_N = np.linalg.inv(alpha*np.eye(M) + beta*PhiT_Phi)
    m_N = beta * np.dot(S_N, np.dot(np.transpose(Phi), t))
    w_opt = m_N
    return (w_opt, alpha, beta, S_N)

##################################################################################

def train(X, y):
    # This function is used for training our Bayesian model
    # Returns the regression parameters w_opt, and alpha, beta, S_N
    # needed for the predictive distribution

    Phi = X # the measurement matrix of the input variables x (i.e., features)
    t   = y # the vector of observations for the target variable

    (N, M) = np.shape(Phi)

    # Init values for  hyper-parameters alpha, beta
    alpha = 5*10**(-3)
    beta = 5
    max_iter = 100
    k = 0

    PhiT_Phi = np.dot(np.transpose(Phi), Phi)
    s = np.linalg.svd(PhiT_Phi, compute_uv=0) # Just get the vector of singular values s

    ab_old = np.array([alpha, beta])
    ab_new = np.zeros((1,2))
    tolerance = 10**-3
    while( k < max_iter and np.linalg.norm(ab_old-ab_new) > tolerance):
        k += 1
        try:
            S_N = np.linalg.inv(alpha*np.eye(M) + beta*PhiT_Phi)
        except np.linalg.LinAlgError as err:
            print  "******************************************************************************************************"
            print "                           ALERT: LinearAlgebra Error detected!"
            print "      CHECK if your measurement matrix is not leading to a singular alpha*np.eye(M) + beta*PhiT_Phi"
            print "                           GOODBYE and see you later. Exiting ..."
            print  "******************************************************************************************************"
            sys.exit(-1)

        m_N = beta * np.dot(S_N, np.dot(np.transpose(Phi), t))
        gamma = sum(beta*s[i]**2 /(alpha + beta*s[i]**2) for i in range(M))
        #
        # update alpha, beta
        #
        ab_old = np.array([alpha, beta])
        alpha = gamma /np.inner(m_N,m_N)
        one_over_beta = 1/(N-gamma) * sum( (t[n] - np.inner(m_N, Phi[n]))**2 for n in range(N))
        beta = 1/one_over_beta
        ab_new = np.array([alpha, beta])

    S_N = np.linalg.inv(alpha*np.eye(M) + beta*PhiT_Phi)
    m_N = beta * np.dot(S_N, np.dot(np.transpose(Phi), t))
    w_opt = m_N
    return (w_opt, alpha, beta, S_N)

##################################################################################