Python math.gamma() Examples

The following are code examples for showing how to use math.gamma(). They are extracted from open source Python projects. You can vote up the examples you like or vote down the ones you don't like. You can also save this page to your account.

Example 1
Project: quadpy   Author: nschloe   File: test_line_segment.py    (license) View Source Project 6 votes vote down vote up
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 
Example 2
Project: quadpy   Author: nschloe   File: test_line_segment.py    (license) View Source Project 6 votes vote down vote up
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 
Example 3
Project: py-prng   Author: czechnology   File: math_tools.py    (license) View Source Project 6 votes vote down vote up
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 
Example 4
Project: EvoloPy-NN   Author: 7ossam81   File: CS.py    (license) View Source Project 6 votes vote down vote up
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 
Example 5
Project: ConceptualSpaces   Author: lbechberger   File: concept.py    (license) View Source Project 5 votes vote down vote up
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 
Example 6
Project: quadpy   Author: nschloe   File: helpers.py    (license) View Source Project 5 votes vote down vote up
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]) 
Example 7
Project: quadpy   Author: nschloe   File: helpers.py    (license) View Source Project 5 votes vote down vote up
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) 
Example 8
Project: quadpy   Author: nschloe   File: tools.py    (license) View Source Project 5 votes vote down vote up
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 
Example 9
Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    (license) View Source Project 5 votes vote down vote up
def lower_incomplete_gamma2(a,x):
    return gamma(a)-upper_incomplete_gamma2(a,x) 
Example 10
Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    (license) View Source Project 5 votes vote down vote up
def gammainc(a,x):
    return lower_incomplete_gamma(a,x)/gamma(a) 
Example 11
Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    (license) View Source Project 5 votes vote down vote up
def gammaincc(a,x):
    return upper_incomplete_gamma(a,x)/gamma(a) 
Example 12
Project: PyBGMM   Author: junlulocky   File: gDirichlet.py    (license) View Source Project 5 votes vote down vote up
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]) 
Example 13
Project: PyBGMM   Author: junlulocky   File: gDirichlet.py    (license) View Source Project 5 votes vote down vote up
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 
Example 14
Project: PyBGMM   Author: junlulocky   File: Dirichlet.py    (license) View Source Project 5 votes vote down vote up
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]) 
Example 15
Project: TTP-Dev   Author: samemu   File: blr_v12.py    (license) View Source Project 5 votes vote down vote up
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)))

################################################################################## 
Example 16
Project: TTP-Dev   Author: samemu   File: blr_v11.py    (license) View Source Project 5 votes vote down vote up
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)))

################################################################################## 
Example 17
Project: binary_swarm_intelligence   Author: Sanbongawa   File: binary_optimization.py    (license) View Source Project 5 votes vote down vote up
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) 
Example 18
Project: binary_swarm_intelligence   Author: Sanbongawa   File: binary_optimization.py    (license) View Source Project 5 votes vote down vote up
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 
Example 19
Project: binary_swarm_intelligence   Author: Sanbongawa   File: binary_optimization_multi.py    (license) View Source Project 5 votes vote down vote up
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) 
Example 20
Project: binary_swarm_intelligence   Author: Sanbongawa   File: binary_optimization_multi.py    (license) View Source Project 5 votes vote down vote up
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 
Example 21
Project: orthopy   Author: nschloe   File: test_line_orthopy.py    (license) View Source Project 5 votes vote down vote up
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 
Example 22
Project: Data-Science   Author: titu1994   File: Probability.py    (license) View Source Project 5 votes vote down vote up
def B(alpha, beta):
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 23
Project: py-prng   Author: czechnology   File: test_math_tools.py    (license) View Source Project 5 votes vote down vote up
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) 
Example 24
Project: py-prng   Author: czechnology   File: test_math_tools.py    (license) View Source Project 5 votes vote down vote up
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)))) 
Example 25
Project: py-prng   Author: czechnology   File: test_math_tools.py    (license) View Source Project 5 votes vote down vote up
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)))) 
Example 26
Project: py-prng   Author: czechnology   File: math_tools.py    (license) View Source Project 5 votes vote down vote up
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) 
Example 27
Project: py-prng   Author: czechnology   File: math_tools.py    (license) View Source Project 5 votes vote down vote up
def incomplete_gamma_upper_norm(a, x, precision=9):
    return incomplete_gamma_upper(a, x, precision) / gamma(a) 
Example 28
Project: Python-Machine-Learning-Algorithm   Author: zhaozhiyong19890102   File: dbscan.py    (license) View Source Project 5 votes vote down vote up
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 
Example 29
Project: TTP-Dev   Author: samemu   File: blr_v12.py    (license) View Source Project 4 votes vote down vote up
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)

################################################################################## 
Example 30
Project: TTP-Dev   Author: samemu   File: blr_v11.py    (license) View Source Project 4 votes vote down vote up
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)

################################################################################## 
Example 31
Project: binary_swarm_intelligence   Author: Sanbongawa   File: binary_optimization.py    (license) View Source Project 4 votes vote down vote up
def BFFA(Eval_Func,n=20,m_i=25,minf=0,dim=None,prog=False,gamma=1.0,beta=0.20,alpha=0.25):
    """
    input:{ Eval_Func: Evaluate_Function, type is class
            n: Number of population, default=20
            m_i: Number of max iteration, default=300
            minf: minimazation flag, default=0, 0=maximization, 1=minimazation
            dim: Number of feature, default=None
            prog: Do you want to use a progress bar?, default=False
            }

    output:{Best value: type float 0.967
            Best position: type list(int) [1,0,0,1,.....]
            Nunber of 1s in best position: type int [0,1,1,0,1] ? 3
            }
    """
    estimate=Eval_Func().evaluate
    if dim==None:
        dim=Eval_Func().check_dimentions(dim)
    #flag=dr
    global_best=float("-inf") if minf == 0 else float("inf")
    pb=float("-inf") if minf == 0 else float("inf")

    global_position=tuple([0]*dim)
    gen=tuple([0]*dim)
    #gamma=1.0
    #beta=0.20
    #alpha=0.25
    gens_dict = {tuple([0]*dim):float("-inf") if minf == 0 else float("inf")}
    #gens_dict[global_position]=0.001
    gens=random_search(n,dim)
    #vs = [[random.choice([0,1]) for i in range(length)] for i in range(N)]
    for gen in gens:
        if tuple(gen) in gens_dict:
            score = gens_dict[tuple(gen)]
        else:
            score=estimate(gen)
            gens_dict[tuple(gen)]=score
        if score > global_best:
            global_best=score
            global_position=dc(gen)
    if prog:
        miter=tqdm(range(m_i))
    else:
        miter=range(m_i)
    for it in miter:
        for i,x in enumerate(gens):
            for j,y in enumerate(gens):
                if gens_dict[tuple(y)] < gens_dict[tuple(x)]:
                    gens[j]=exchange_binary(y,gens_dict[tuple(y)])
                gen = gens[j]
                if tuple(gen) in gens_dict:
                    score = gens_dict[tuple(gen)]
                else:
                    score=estimate(gens[j])
                    gens_dict[tuple(gen)]=score
                if score > global_best if minf==0 else score < global_best:
                    global_best=score
                    global_position=dc(gen)
    return global_best,global_position,global_position.count(1) 
Example 32
Project: orthopy   Author: nschloe   File: test_line_orthopy.py    (license) View Source Project 4 votes vote down vote up
def test_gautschi_how_to_and_how_not_to():
    '''Test Gautschi's famous example from

    W. Gautschi,
    How and how not to check Gaussian quadrature formulae,
    BIT Numerical Mathematics,
    June 1983, Volume 23, Issue 2, pp 209–216,
    <https://doi.org/10.1007/BF02218441>.
    '''
    points = numpy.array([
        1.457697817613696e-02,
        8.102669876765460e-02,
        2.081434595902250e-01,
        3.944841255669402e-01,
        6.315647839882239e-01,
        9.076033998613676e-01,
        1.210676808760832,
        1.530983977242980,
        1.861844587312434,
        2.199712165681546,
        2.543839804028289,
        2.896173043105410,
        3.262066731177372,
        3.653371887506584,
        4.102376773975577,
        ])
    weights = numpy.array([
        3.805398607861561e-2,
        9.622028412880550e-2,
        1.572176160500219e-1,
        2.091895332583340e-1,
        2.377990401332924e-1,
        2.271382574940649e-1,
        1.732845807252921e-1,
        9.869554247686019e-2,
        3.893631493517167e-2,
        9.812496327697071e-3,
        1.439191418328875e-3,
        1.088910025516801e-4,
        3.546866719463253e-6,
        3.590718819809800e-8,
        5.112611678291437e-11,
        ])

    # weight function exp(-t**3/3)
    n = len(points)
    moments = numpy.array([
        3.0**((k-2)/3.0) * math.gamma((k+1) / 3.0)
        for k in range(2*n)
        ])

    alpha, beta = orthopy.line.coefficients_from_gauss(points, weights)
    # alpha, beta = orthopy.line.chebyshev(moments)

    errors_alpha, errors_beta = \
        orthopy.line.check_coefficients(moments, alpha, beta)

    assert numpy.max(errors_alpha) > 1.0e-2
    assert numpy.max(errors_beta) > 1.0e-2
    return