我们从Python开源项目中，提取了以下32个代码示例，用于说明如何使用math.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 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 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)
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
