Python math.gamma() Examples

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

Example 1
Project: bruno   Author: christabella   File: utils_stats.py    MIT License 6 votes vote down vote up
def mvt_pdf(X, phi, K, nu):
    """
    Multivariate student-t density:
    output:
        the density of the given element
    input:
        x = parameter (d dimensional numpy array or scalar)
        mu = mean (d dimensional numpy array or scalar)
        K = scale matrix (dxd numpy array)
        nu = degrees of freedom
    """
    d = X.shape[-1]
    num = math.gamma((d + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * ((X - phi).dot(np.linalg.inv(K)).dot(np.transpose(X - phi))), -(d + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi, d / 2.) * pow(np.linalg.det(K), 0.5)
    return num / denom 
Example 2
Project: spotpy   Author: thouska   File: likelihoods.py    MIT License 6 votes vote down vote up
def AR_1_Coeff(data):
        """
        The autocovariance coefficient called as rho, for an AR(1) model can be calculated as shown here:

        .. math::

            \\rho(1) = \\frac{\\gamma(1)}{\\gamma(0)}

        For further information look for example in "Zeitreihenanalyse", pages 17, by Matti Schneider, Sebastian Mentemeier,
        SS 2010.

        :param data: numerical list
        :type data: list
        :return: autocorrelation coefficient
        :rtype: float
        """
        return TimeSeries.acf(data, 1) / TimeSeries.acf(data, 0) 
Example 3
Project: quadpy   Author: nschloe   File: test_line_segment.py    MIT License 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: scheme.integrate(poly, numpy.array([[-1.0], [1.0]])),
        integrate_exact,
        scheme.degree + 1,
    )
    assert degree >= scheme.degree
    return 
Example 4
Project: quadpy   Author: nschloe   File: test_tools.py    MIT License 6 votes vote down vote up
def test_integrate():
    moments = quadpy.tools.integrate(lambda x: [x ** k for k in range(5)], -1, +1)
    assert (moments == [2, 0, sympy.S(2) / 3, 0, sympy.S(2) / 5]).all()

    moments = quadpy.tools.integrate(
        lambda x: orthopy.line_segment.tree_legendre(x, 4, "monic", symbolic=True),
        -1,
        +1,
    )
    assert (moments == [2, 0, 0, 0, 0]).all()

    # Example from Gautschi's "How to and how not to" article
    moments = quadpy.tools.integrate(
        lambda x: [x ** k * sympy.exp(-(x ** 3) / 3) for k in range(5)], 0, sympy.oo
    )
    S = numpy.vectorize(sympy.S)
    gamma = numpy.vectorize(sympy.gamma)
    n = numpy.arange(5)
    reference = 3 ** (S(n - 2) / 3) * gamma(S(n + 1) / 3)
    assert numpy.all([sympy.simplify(m - r) == 0 for m, r in zip(moments, reference)]) 
Example 5
Project: lenstronomy   Author: sibirrer   File: jeans_equation.py    MIT License 6 votes vote down vote up
def power_law_anisotropy(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        """
        equation (19) in Suyu+ 2010
        :param r:
        :return:
        """
        # first term
        theta_E = kwargs_profile['theta_E']
        gamma = kwargs_profile['gamma']
        r_ani = kwargs_anisotropy['r_ani']
        a = 0.551 * kwargs_light['r_eff']
        rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
        prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
        prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
        hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
        hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
        fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
        sigma2_dim_less = prefac1 * prefac2 * fac
        return sigma2_dim_less * (self.cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000)**2 
Example 6
Project: lenstronomy   Author: sibirrer   File: analytic_kinematics.py    MIT License 6 votes vote down vote up
def vel_disp(self, gamma, theta_E, r_eff, r_ani, rendering_number=1000):
        """
        computes the averaged LOS velocity dispersion in the slit (convolved)

        :param gamma: power-law slope of the mass profile (isothermal = 2)
        :param theta_E: Einstein radius of the lens (in arcseconds)
        :param r_eff: half light radius of the Hernquist profile (or as an approximation of any other profile to be described as a Hernquist profile
        :param r_ani: anisotropy radius
        :param kwargs_aperture: keyword arguments describing the aperture of the collected spectral
        :param rendering_number: number of spectral renderings drawn from the light distribution that go through the
            slit of the observations

        :return: LOS integrated velocity dispersion in units [km/s]
        """
        sigma_s2_sum = 0
        rho0_r0_gamma = self._rho0_r0_gamma(theta_E, gamma)
        for i in range(0, rendering_number):
            sigma_s2_draw = self.vel_disp_one(gamma, rho0_r0_gamma, r_eff, r_ani)
            sigma_s2_sum += sigma_s2_draw
        sigma_s2_average = sigma_s2_sum / rendering_number
        return np.sqrt(sigma_s2_average) 
Example 7
Project: lenstronomy   Author: sibirrer   File: analytic_kinematics.py    MIT License 6 votes vote down vote up
def vel_disp_one(self, gamma, rho0_r0_gamma, r_eff, r_ani):
        """
        computes one realisation of the velocity dispersion realized in the slit

        :param gamma: power-law slope of the mass profile (isothermal = 2)
        :param rho0_r0_gamma: combination of Einstein radius and power-law slope as equation (14) in Suyu+ 2010
        :param r_eff: half light radius of the Hernquist profile (or as an approximation of any other profile to be described as a Hernquist profile
        :param r_ani: anisotropy radius
        :param kwargs_aperture: keyword arguments describing the aperture of the collected spectral
        :param FWHM: full width at half maximum of the seeing conditions, described as a Gaussian
        :return: projected velocity dispersion of a single drawn position in the potential [km/s]
        """
        a = 0.551 * r_eff
        while True:
            r = self.P_r(a)  # draw r
            R, x, y = self.R_r(r)  # draw projected R
            x_, y_ = self._psf.displace_psf(x, y)
            bool = self.aperture.aperture_select(x_, y_)
            if bool is True:
                break
        sigma_s2 = self.sigma_s2(r, R, r_ani, a, gamma, rho0_r0_gamma)
        return sigma_s2 
Example 8
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 6 votes vote down vote up
def __init__(self, params=None, core_model='moffat', aureole_model='power'):
        """
        Parameters
        ----------
        params : a dictionary containing keywords of PSF parameter
        core_model : model of PSF core (moffat)
        aureole_model : model of aureole ("power" or "multi-power")	
        
        """
        self.core_model = core_model
        self.aureole_model = aureole_model
        
        self.params = params
        
        # Build attribute for parameters from dictionary keys 
        for key, val in params.items():
            exec('self.' + key + ' = val')
            
        if hasattr(self, 'fwhm'):
            self.gamma = fwhm_to_gamma(self.fwhm, self.beta)
            self.params['gamma'] = self.gamma
            
        if hasattr(self, 'gamma'):
            self.fwhm  = gamma_to_fwhm(self.gamma, self.beta)
            self.params['fwhm'] = self.fwhm 
Example 9
Project: bruno   Author: IraKorshunova   File: utils_stats.py    MIT License 6 votes vote down vote up
def mvt_pdf(X, phi, K, nu):
    """
    Multivariate student-t density:
    output:
        the density of the given element
    input:
        x = parameter (d dimensional numpy array or scalar)
        mu = mean (d dimensional numpy array or scalar)
        K = scale matrix (dxd numpy array)
        nu = degrees of freedom
    """
    d = X.shape[-1]
    num = math.gamma((d + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * ((X - phi).dot(np.linalg.inv(K)).dot(np.transpose(X - phi))), -(d + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi, d / 2.) * pow(np.linalg.det(K), 0.5)
    return num / denom 
Example 10
Project: tripp   Author: mjamesruggiero   File: hypothesis.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def B(alpha, beta):
    """a normalizing constant so that the total probability is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 11
Project: bruno   Author: christabella   File: utils_stats.py    MIT License 5 votes vote down vote up
def student_pdf_1d(X, phi, K, nu):
    """
    Univariate student-t density:
    output:
        the density of the given element
    input:
        x = parameter scalar
        mu = mean
        var = scale matrix
        nu = degrees of freedom
    """
    num = math.gamma((1. + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * (1. / K * (X - phi) ** 2), -(1. + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * K, 0.5)
    return num / denom 
Example 12
Project: bruno   Author: christabella   File: test_latent_dist.py    MIT License 5 votes vote down vote up
def student_pdf_1d(X, mu, var, nu):
    if nu > 50:
        return gauss_pdf_1d(X, mu, var)
    num = math.gamma((1. + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * (1. / var * (X - mu) ** 2), -(1. + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * var, 0.5)
    return num / denom 
Example 13
Project: GPEXP   Author: goroda   File: utilities.py    GNU General Public License v2.0 5 votes vote down vote up
def distribFuncNDCircle(points, radius):
    nDim = points.shape[1]
    area = np.pi**(nDim/2.0)/math.gamma(nDim/2.0 + 1)*radius**nDim
    prob = 1.0/area
    out = np.zeros((len(points)))
    out[points[:,0]**2.0 + points[:,1]**2.0 < radius**2.0] = prob

    return out 
Example 14
Project: wind-turbine-mdo   Author: lewisli   File: WindDistribution.py    MIT License 5 votes vote down vote up
def ComputeScaleFunction(Mean, Shape):
	return Mean/math.gamma(1 + 1.0/Shape) 
Example 15
Project: ironpython2   Author: IronLanguages   File: test_math.py    Apache License 2.0 5 votes vote down vote up
def test_gamma(self):
        self.assertAlmostEqual(math.gamma(0.5), math.sqrt(math.pi), places=15)
        for i in xrange(1, 20):
            self.assertEqual(math.factorial(i-1), math.gamma(i))
        self.assertEqual(math.gamma(float('inf')), float('inf'))
        self.assertRaises(ValueError, math.gamma, float('-inf'))
        self.assertTrue(math.isnan(math.gamma(float('nan'))))
        for i in xrange(0, -1001, -1):
            self.assertRaises(ValueError, math.gamma, i) 
Example 16
Project: ironpython2   Author: IronLanguages   File: test_math.py    Apache License 2.0 5 votes vote down vote up
def test_gamma_lgamma(self):
        tolerance = 13
        for x in itertools.count(0.001, 0.001):
            if x > 5.0:
                break
            self.assertAlmostEqual(math.lgamma(x), math.log(math.gamma(x)), places=tolerance)
            self.assertAlmostEqual(math.lgamma(x*x), math.log(math.gamma(x*x)), places=tolerance)
            self.assertAlmostEqual(math.lgamma(2.0**x), math.log(math.gamma(2.0**x)), places=tolerance)

            # Test negative values too, but not integers
            if x % 1.0 != 0.0:
                self.assertAlmostEqual(math.lgamma(-x), math.log(abs(math.gamma(-x))), places=tolerance)
                self.assertAlmostEqual(math.lgamma(-x*x), math.log(abs(math.gamma(-x*x))), places=tolerance)
                self.assertAlmostEqual(math.lgamma(-2.0**x), math.log(abs(math.gamma(-2.0**x))), places=tolerance) 
Example 17
Project: medium   Author: yusueliu   File: plot_dirichlet.py    MIT License 5 votes vote down vote up
def __init__(self, alpha):
        
        self._alpha = np.array(alpha)
        self._coef = gamma(np.sum(self._alpha)) / reduce(mul, [gamma(a) for a in self._alpha]) 
Example 18
Project: ParticleFlowBayesRule   Author: xinshi-chen   File: distributions.py    MIT License 5 votes vote down vote up
def get_gamma(X, bandwidth):
    with torch.no_grad():
        x_norm = torch.sum(X ** 2, dim=1, keepdim=True)
        x_t = torch.transpose(X, 0, 1)
        x_norm_t = x_norm.view(1, -1)
        t = x_norm + x_norm_t - 2.0 * torch.matmul(X, x_t)
        dist2 = F.relu(Variable(t)).detach().data

        d = dist2.cpu().numpy()
        d = d[np.isfinite(d)]
        d = d[d > 0]
        median_dist2 = float(np.median(d))
        gamma = 0.5 / median_dist2 / bandwidth
        return gamma 
Example 19
Project: ParticleFlowBayesRule   Author: xinshi-chen   File: distributions.py    MIT License 5 votes vote down vote up
def get_kernel_mat(x, landmarks, gamma):
    d = pairwise_distances(x, landmarks)
    k = torch.exp(d * -gamma)
    k = k.view(x.shape[0], -1)
    return k 
Example 20
Project: ParticleFlowBayesRule   Author: xinshi-chen   File: distributions.py    MIT License 5 votes vote down vote up
def MMD(x, y, bandwidth=1.0):
    y = y.detach()
    gamma = get_gamma(y.detach(), bandwidth)
    kxx = get_kernel_mat(x, x, gamma)
    idx = torch.arange(0, x.shape[0], out=torch.LongTensor())
    kxx = kxx * (1 - torch.eye(x.shape[0]).to(DEVICE))
    kxx = torch.sum(kxx) / x.shape[0] / (x.shape[0] - 1)

    kyy = get_kernel_mat(y, y, gamma)
    idx = torch.arange(0, y.shape[0], out=torch.LongTensor())
    kyy[idx, idx] = 0.0
    kyy = torch.sum(kyy) / y.shape[0] / (y.shape[0] - 1)
    kxy = torch.sum(get_kernel_mat(y, x, gamma)) / x.shape[0] / y.shape[0]
    mmd = kxx + kyy - 2 * kxy
    return mmd 
Example 21
Project: ParticleFlowBayesRule   Author: xinshi-chen   File: distributions.py    MIT License 5 votes vote down vote up
def log_pdf(self, log_x):
        return (self.alpha - 1) * log_x - self.beta * torch.exp(log_x) + self.alpha * np.log(self.beta) - np.log(math.gamma(self.alpha)) 
Example 22
Project: ParticleFlowBayesRule   Author: xinshi-chen   File: distributions.py    MIT License 5 votes vote down vote up
def get_samples(self, num_samples):
        rand_vars = np.random.gamma(self.alpha, 1.0 / self.beta, size=[num_samples, 1])  # numpy uses scale parameter
        return torch.tensor(np.log(rand_vars), 
                            dtype=t_float).to(DEVICE) 
Example 23
Project: pynam   Author: hbp-unibi   File: entropy.py    GNU General Public License v3.0 5 votes vote down vote up
def ncrr(x, y):
    """
    Implementation of the binomial coefficient (ncr) function for real arguments
    x and y (respectively corresponding to n and r).
    """
    return math.gamma(x+1.0) / (math.gamma(y+1.0) * math.gamma(x-y+1.0)) 
Example 24
Project: data-science-from-scratch   Author: kevntao   File: hypothesis_and_inference.py    The Unlicense 5 votes vote down vote up
def B(alpha, beta):
    """a normalizing constant so that the total probability is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 25
Project: mne-features   Author: mne-tools   File: univariate.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _hurst_exp_helper(x, n_splits=20):
    """Helper function for :func:`compute_hurst_exp`.

    Compute the Hurst exponent from a univariate time series. The Hurst
    exponent is defined as the slope of the least-squares regression line
    going through a cloud of `n_splits` points. Each point is obtained by
    considering sub-series of `x` of `n_splits` different lenghts.

    Parameters
    ----------
    x : ndarray, shape (n_times,)

    Returns
    -------
    output : ndarray, shape (n_splits,)
    """
    n_times = x.shape[0]
    _splits = np.floor(np.logspace(start=4, stop=np.log2(n_times / 2),
                                   num=n_splits, base=2.))
    splits = np.unique(_splits).astype(int)
    reg = np.zeros((splits.size,))
    for j, n in enumerate(splits):
        a = x.copy()
        d = int(floor(n_times / n))
        a = np.lib.stride_tricks.as_strided(a, shape=(d, n),
                                            strides=(n * a.strides[-1],
                                                     a.strides[-1]))
        _rs = _hurst_exp_compute_rs(a)
        _rs = _rs[~np.isnan(_rs)]
        reg[j] = np.log(np.mean(_rs))
        s = sum([sqrt((n - i) / i) for i in range(1, n)]) * ((n - 0.5) / n)
        if n <= 340:
            corr = (gamma((n - 1) / 2.) / (sqrt(np.pi) * gamma(n / 2.))) * s
        else:
            corr = ((n - 0.5) / n) * (1. / sqrt(np.pi * n / 2.)) * s
        reg[j] -= (np.log(corr) - np.log(n) / 2)
    return _slope_lstsq(np.log(splits), reg) 
Example 26
Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    GNU General Public License v2.0 5 votes vote down vote up
def lower_incomplete_gamma2(a,x):
    return gamma(a)-upper_incomplete_gamma2(a,x) 
Example 27
Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    GNU General Public License v2.0 5 votes vote down vote up
def gammainc(a,x):
    return lower_incomplete_gamma(a,x)/gamma(a) 
Example 28
Project: sp800_22_tests   Author: dj-on-github   File: gamma_functions.py    GNU General Public License v2.0 5 votes vote down vote up
def gammaincc(a,x):
    return upper_incomplete_gamma(a,x)/gamma(a) 
Example 29
Project: python3_ios   Author: holzschu   File: bayes_update.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def beta_pdf(x, a, b):
    return (x**(a-1) * (1-x)**(b-1) * math.gamma(a + b)
            / (math.gamma(a) * math.gamma(b))) 
Example 30
Project: BNMTF   Author: ThomasBrouwer   File: test_bnmtf_vb_optimised.py    Apache License 2.0 5 votes vote down vote up
def test_update_exp_F():
    for k in range(0,K):
        BNMTF = bnmtf_vb_optimised(R,M,K,L,priors)
        BNMTF.initialise(init_S,init_FG)
        BNMTF.tauF = 4*numpy.ones((I,K))  # muF = [[0.5]], tauF = [[4.]]
        BNMTF.update_exp_F(k) #-mu*sqrt(tau) = -0.5*2 = -1. lambda(1) = 0.241971 / (1-0.1587) = 0.2876155949126352. gamma = 0.37033832534958433
        for i in range(0,I):        
            assert abs(BNMTF.expF[i,k] - (0.5 + 1./2. * 0.2876155949126352)) < 0.00001
            assert abs(BNMTF.varF[i,k] - 1./4.*(1.-0.37033832534958433)) < 0.00001 
Example 31
Project: BNMTF   Author: ThomasBrouwer   File: test_bnmtf_vb_optimised.py    Apache License 2.0 5 votes vote down vote up
def test_update_exp_S():
    for k,l in itertools.product(xrange(0,K),xrange(0,L)):
        BNMTF = bnmtf_vb_optimised(R,M,K,L,priors)
        BNMTF.initialise(init_S,init_FG) 
        BNMTF.tauS = 4*numpy.ones((K,L)) # muS = [[1./3.]], tauS = [[4.]]
        BNMTF.update_exp_S(k,l) #-mu*sqrt(tau) = -2./3., lambda(..) = 0.319448 / (1-0.2525) = 0.4273551839464883, gamma = 
        assert abs(BNMTF.expS[k,l] - (1./3. + 1./2. * 0.4273551839464883)) < 0.00001
        assert abs(BNMTF.varS[k,l] - 1./4.*(1. - 0.4675359092102624)) < 0.00001 
Example 32
Project: BNMTF   Author: ThomasBrouwer   File: test_bnmtf_vb_optimised.py    Apache License 2.0 5 votes vote down vote up
def test_update_exp_G():
    for l in range(0,L):
        BNMTF = bnmtf_vb_optimised(R,M,K,L,priors)
        BNMTF.initialise(init_S,init_FG) 
        BNMTF.tauG = 4*numpy.ones((J,L)) # muG = [[1./4.]], tauG = [[4.]]
        BNMTF.update_exp_G(l) #-mu*sqrt(tau) = -0.5., lambda(..) = 0.352065 / (1-0.3085) = 0.5091323210412148, gamma = 0.5137818808494219
        for j in range(0,J):        
            assert abs(BNMTF.expG[j,l] - (1./4. + 1./2. * 0.5091323210412148)) < 0.0001
            assert abs(BNMTF.varG[j,l] - 1./4.*(1. - 0.5137818808494219)) < 0.0001 
Example 33
Project: BNMTF   Author: ThomasBrouwer   File: test_bnmf_vb_optimised.py    Apache License 2.0 5 votes vote down vote up
def test_update_exp_U():
    for k in range(0,K):
        BNMF = bnmf_vb_optimised(R,M,K,priors)
        BNMF.initialise()
        BNMF.tauU = 4*numpy.ones((I,K)) # muU = [[0.5]], tauU = [[4.]]
        BNMF.update_exp_U(k) #-mu*sqrt(tau) = -0.5*2 = -1. lambda(1) = 0.241971 / (1-0.1587) = 0.2876155949126352. gamma = 0.37033832534958433
        for i in range(0,I):        
            assert abs(BNMF.expU[i,k] - (0.5 + 1./2. * 0.2876155949126352)) < 0.00001
            assert abs(BNMF.varU[i,k] - 1./4.*(1.-0.37033832534958433)) < 0.00001 
Example 34
Project: paradoc   Author: betaveros   File: discrete.py    MIT License 5 votes vote down vote up
def factorial(n: Union[int, float]) -> Union[int, float]:
    try:
        import sympy.functions.combinatorial.factorials as fs
    except ModuleNotFoundError:
        if isinstance(n, int):
            p = 1
            for i in range(1, n + 1): p *= i
            return p
        else:
            return math.gamma(n + 1)

    return int(fs.factorial(n)) 
Example 35
Project: geometer   Author: jan-mue   File: curve.py    MIT License 5 votes vote down vote up
def _alpha(n):
        return math.pi**(n/2) / math.gamma(n/2 + 1) 
Example 36
Project: PyOpenDial   Author: KAIST-AILab   File: math_utils.py    MIT License 5 votes vote down vote up
def log_gamma(x):
        """
        Returns the log-gamma value using Lanczos approximation formula

        :param x: the point
        :return: the log-gamma value for the point
        """
        return math.lgamma(x) 
Example 37
Project: PyOpenDial   Author: KAIST-AILab   File: math_utils.py    MIT License 5 votes vote down vote up
def gamma(x):
        """
        Returns the value of the gamma function: Gamma(x) = integral( t^(x-1) e^(-t),
        t = 0 .. infinity)

        :param x: the point
        :return: the gamma value for the point
        """
        return math.gamma(x) 
Example 38
Project: PyOpenDial   Author: KAIST-AILab   File: math_utils.py    MIT License 5 votes vote down vote up
def get_volume(radius, dim):
        """
        Returns the volume of an N-ball of a certain radius.

        :param radius: the radius
        :param dim: the number of dimensions to consider
        :return: the resulting volume
        """
        numerator = math.pow(math.pi, dim / 2.)
        denominator = MathUtils.gamma((dim / 2.) + 1)
        radius2 = math.pow(radius, dim)
        return radius2 * numerator / denominator 
Example 39
Project: quadpy   Author: nschloe   File: _helpers.py    MIT License 5 votes vote down vote up
def plot(self, show_axes=True):
        import matplotlib.pyplot as plt

        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)

        plot_disks(plt, self.points, self.weights, I0)
        return 
Example 40
Project: quadpy   Author: nschloe   File: _helpers.py    MIT License 5 votes vote down vote up
def plot(self, show_axes=True):
        import matplotlib.pyplot as plt

        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)

        plot_disks(plt, self.points, self.weights, I0)
        return 
Example 41
Project: quadpy   Author: nschloe   File: weights_from_least_squares.py    MIT License 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)
    )


# def simplex_monomials(degree):
#     exponents = numpy.concatenate(
#         [quadpy.helpers.partition(d, 2) for d in range(degree + 1)]
#     )
#
#     exact_vals = numpy.array(
#         [integrate_monomial_over_unit_simplex(k) for k in exponents]
#     )
#
#     def fun(x):
#         k = exponents.T
#         # <https://stackoverflow.com/a/46689653/353337>
#         s = x.shape[1:] + k.shape[1:]
#         return (
#             (x.reshape(x.shape[0], -1, 1) ** k.reshape(k.shape[0], 1, -1))
#             .prod(0)
#             .reshape(s)
#         )
#
#     return fun, exact_vals 
Example 42
Project: quadpy   Author: nschloe   File: improve_precision.py    MIT License 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 43
Project: quadpy   Author: nschloe   File: helpers.py    MIT License 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 44
Project: quadpy   Author: nschloe   File: helpers.py    MIT License 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 45
Project: PyFlow   Author: wonderworks-software   File: MathLib.py    Apache License 2.0 5 votes vote down vote up
def gamma(x=('FloatPin', 0.0), Result=(REF, ('BoolPin', False))):
        '''Return the Gamma function at `x`.'''
        try:
            Result(True)
            return math.gamma(x)
        except:
            Result(False)
            return -1 
Example 46
Project: Machine-Learning-with-Python   Author: devAmoghS   File: hypothesis_inference.py    MIT License 5 votes vote down vote up
def B(alpha, beta):
    """a normalizing constant so that the total probability is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta) 
Example 47
Project: m-stp   Author: MukeunKim   File: core.py    MIT License 5 votes vote down vote up
def gamma(x):
    if isinstance(x, LineValue): lx = x.get_value()
    else: lx = x
    if lx == NAN: return LineValue(NAN)
    return LineValue(math.gamma(lx)) 
Example 48
Project: lenstronomy   Author: sibirrer   File: jeans_equation.py    MIT License 5 votes vote down vote up
def _rho0_r0_gamma(self, theta_E, gamma):
        # equation (14) in Suyu+ 2010
        return -1 * math.gamma(gamma/2.)/(np.sqrt(np.pi)*math.gamma((gamma-3)/2.)) * theta_E**gamma/self.cosmo.arcsec2phys_lens(theta_E) * self.cosmo.epsilon_crit * const.M_sun/const.Mpc**3  # units kg/m^3 
Example 49
Project: lenstronomy   Author: sibirrer   File: analytic_kinematics.py    MIT License 5 votes vote down vote up
def _rho0_r0_gamma(self, theta_E, gamma):
        # equation (14) in Suyu+ 2010
        return -1 * math.gamma(gamma/2) / (np.sqrt(np.pi)*math.gamma((gamma-3)/2.)) * theta_E ** gamma / \
               self._cosmo.arcsec2phys_lens(theta_E) * self._cosmo.epsilon_crit * const.M_sun / const.Mpc ** 3 
Example 50
Project: lenstronomy   Author: sibirrer   File: analytic_kinematics.py    MIT License 5 votes vote down vote up
def sigma_r2(self, r, a, gamma, rho0_r0_gamma, r_ani):
        """
        equation (19) in Suyu+ 2010
        """
        # first term
        prefac1 = 4*np.pi * const.G * a**(-gamma) * rho0_r0_gamma / (3-gamma)
        prefac2 = r * (r + a)**3/(r**2 + r_ani**2)
        hyp1 = vel_util.hyp_2F1(a=2+gamma, b=gamma, c=3+gamma, z=1./(1+r/a))
        hyp2 = vel_util.hyp_2F1(a=3, b=gamma, c=1+gamma, z=-a/r)
        fac = r_ani**2/a**2 * hyp1 / ((2+gamma) * (r/a + 1)**(2+gamma)) + hyp2 / (gamma*(r/a)**gamma)
        return prefac1 * prefac2 * fac * (self._cosmo.arcsec2phys_lens(1.) * const.Mpc / 1000) ** 2 
Example 51
Project: cms   Author: broadinstitute   File: stats.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def gammainc_halfint(s, x) :
    """ Lower incomplete gamma function = 
            integral from 0 to x of t ** (s-1) exp(-t) dt divided by gamma(s),
        i.e., the fraction of gamma that you get if you integrate only until
            x instead of all the way to infinity.
        Implemented here only if s is a positive multiple of 0.5.
    """
    # scipy equivalent: scipy.special.gammainc(s,x)
    
    if s <= 0 :
        raise ValueError('%s is not positive' % s)
    if x < 0 :
        raise ValueError('%s < 0' % x)
    if s * 2 != int(s * 2) :
        raise NotImplementedError('%s is not a multiple of 0.5' % s)
    
    # Handle integers analytically
    if s == int(s) :
        term = 1
        total = 1
        for k in range(1, int(s)) :
            term *= x /  k
            total += term
        return 1 - exp(-x) * total

    # Otherwise s is integer + 0.5. Decrease to 0.5 using recursion formula:
    result = 0.0
    while s > 1 :
        result -= x ** (s - 1) * exp(-x) / gamma(s)
        s = s - 1
    # Then use gammainc(0.5, x) = erf(sqrt(x))
    result += erf(sqrt(x))
    return result 
Example 52
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def make_grid(self, image_size, pixel_scale=2.5):
        """ Build grid for drawing """
        self.image_size = image_size
        self.yy, self.xx = np.mgrid[:image_size, :image_size]
        self.cen = ((image_size-1)/2., (image_size-1)/2.)
        self.pixel_scale = pixel_scale
        
        for key, val in self.params.items():
            if (key == 'gamma') | ('theta' in key):
                val = val / pixel_scale
                exec('self.' + key + '_pix' + ' = val') 
Example 53
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def generate_core(self, folding_threshold=1.e-10):
        """ Generate Galsim PSF of core. """
        self.fwhm = self.gamma * 2. * np.sqrt(2**(1./self.beta)-1)
        gsparams = galsim.GSParams(folding_threshold=folding_threshold)
        psf_core = galsim.Moffat(beta=self.beta, fwhm=self.fwhm,
                                flux=1., gsparams=gsparams) # in arcsec
        self.psf_core = psf_core
        return psf_core 
Example 54
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def draw_core2D_in_real(self, star_pos, Flux):
        """ 2D drawing function of the core in real space given positions and flux (of core) of target stars """
        from astropy.modeling.models import Moffat2D
        gamma, alpha = self.gamma_pix, self.beta
        Amps = np.array([moffat2d_Flux2Amp(gamma, alpha, Flux=flux)
                       for flux in Flux])
        f_core_2d_s = np.array([Moffat2D(amplitude=amp, x_0=x0, y_0=y0,
                                                gamma=gamma, alpha=alpha)
                                for ((x0,y0), amp) in zip(star_pos, Amps)])
            
        return f_core_2d_s 
Example 55
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def gamma_to_fwhm(gamma, beta):
    """ in arcsec """
    return gamma / fwhm_to_gamma(1, beta) 
Example 56
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def moffat1d_normed(x, gamma, alpha):
    """ Moffat for 1d array, flux normalized = 1 """
    from astropy.modeling.models import Moffat1D
    Mof_mod_1d = Moffat1D(amplitude=1, x_0=0, gamma=gamma, alpha=alpha)
    norm_mof = quad(Mof_mod_1d, 0, np.inf)[0] 
    y = Mof_mod_1d(x) / norm_mof
    return y 
Example 57
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def C_mof2Dto1D(r_core, beta):
    """ gamma in pixel """
    return 1. / (beta-1) * 2 * np.sqrt(np.pi) * r_core * math.gamma(beta) / math.gamma(beta - 1. / 2) 
Example 58
Project: mrf   Author: AstroJacobLi   File: modeling.py    MIT License 5 votes vote down vote up
def C_mof1Dto2D(r_core, beta):
    """ gamma in pixel """
    return 1. / C_mof2Dto1D(r_core, beta) 
Example 59
Project: recipes   Author: cwlseu   File: possion.py    GNU General Public License v3.0 5 votes vote down vote up
def bernonlli(p, n=1, k=1):
	return C(n, k)*(p**k)*((1-p)**(n-k))

# gamma函数是阶乘在实数上的延拓,gamma(x+1) = x*gamma(x) 
Example 60
Project: recipes   Author: cwlseu   File: possion.py    GNU General Public License v3.0 5 votes vote down vote up
def gamma(x):
	if x <= 0:
		raise ValueError("Gamma Parameter Error!")
	if x == 0.5:
		return math.pi**0.5
	else:
		return (x-1)*gamma(x-1) 
Example 61
Project: Turing   Author: TuringApp   File: stats.py    MIT License 5 votes vote down vote up
def binomial(n, k):
    """Calculates the binomial coefficient."""
    return math.gamma(n + 1) / (math.gamma(k + 1) * math.gamma(n - k + 1)) 
Example 62
Project: Turing   Author: TuringApp   File: stats.py    MIT License 5 votes vote down vote up
def gamma(x):
    return math.gamma(x) 
Example 63
Project: Turing   Author: TuringApp   File: stats.py    MIT License 5 votes vote down vote up
def beta(a, b):
    return (gamma(a) * gamma(b)) / gamma(a + b) 
Example 64
Project: opytimizer   Author: gugarosa   File: distribution.py    GNU General Public License v3.0 5 votes vote down vote up
def generate_levy_distribution(beta=0.1, size=1):
    """Generates a n-dimensional array based on a Lévy distribution.

    References:
        X.-S. Yang and S. Deb. Computers & Operations Research.
        Multiobjective Cuckoo Search for Design Optimization (2013).

    Args:
        beta (float): Skewness parameter.
        size (int): Size of array.

    Returns:
        A Lévy distribution n-dimensional array.

    """

    # Calculates the equation's numerator
    num = gamma(1+beta) * sin(pi*beta/2)

    # Calculates the equation's denominator
    den = gamma((1+beta)/2) * beta * (2 ** ((beta-1)/2))

    # Calculates the sigma for further distribution generation
    sigma = (num/den) ** (1/beta)

    # Calculates the 'u' distribution
    u = r.generate_gaussian_random_number(size=size) * sigma

    # Calculates the 'v' distribution
    v = r.generate_gaussian_random_number(size=size)

    # Finally, we can calculate the Lévy distribution
    step = u / np.fabs(v) ** (1 / beta)

    return step 
Example 65
Project: Seriously   Author: Mego   File: SeriouslyCommands.py    MIT License 5 votes vote down vote up
def _factorial(n):
    try:
        return rmath.factorial(n)
    except ValueError:
        return rmath.gamma(n+1) 
Example 66
Project: Functional-Python-Programming-Second-Edition   Author: PacktPublishing   File: chi_sq.py    MIT License 5 votes vote down vote up
def fact(k: int) -> int:
    """Simple factorial.

    >>> fact(1)
    1
    >>> fact(2)
    2
    >>> fact(3)
    6
    >>> fact(4)
    24
    """
    if k < 2:
        return 1
    return reduce(operator.mul, range(2,k+1))

# The implementation uses ``reduce( operator.mul, ... )`` to compute
# the product of a sequence of integer values.
# We've included the ``@lru_cache`` because this is used often,
# and the small domain of possible values leads to some benefit
# from the cache.

# We could also use ``math.factorial()``. In order to make use of
# the cache, we would need to do something like this.
#
# ..  parsed-literal::
#
#     fact = lru_cache(128)(math.factorial)
#
# This would create a similarly cached factorial function.

# Incomplete Gamma
# -----------------

# The incomplete (lower) gamma function is this:
#
# ..  math::
#
#     \gamma(s,z) = \sum_{k=0}^{\infty} \dfrac {(-1)^k} {k!} \; \dfrac {z^{s+k}} {s+k}
# 
Example 67
Project: Functional-Python-Programming-Second-Edition   Author: PacktPublishing   File: chi_sq.py    MIT License 5 votes vote down vote up
def Gamma_Half(k: float) -> float:
    """Gamma(k) with special case for k = n+1/2; k-1/2=n.

    >>> import math
    >>> round(Gamma_Half(2),1)
    1.0
    >>> round(Gamma_Half(3),1)
    2.0
    >>> round(Gamma_Half(4),1)
    6.0
    >>> round(Gamma_Half(5),1)
    24.0

    >>> round(Gamma_Half(.5), 7)
    1.7724539
    >>> round(math.sqrt(math.pi), 7)
    1.7724539
    >>> round(Gamma_Half(1.5), 7)
    0.8862269
    >>> round(math.sqrt(math.pi)/2, 7)
    0.8862269
    """
    ε = 1E-6
    if abs(k-int(k)-.5) < ε:
        n = int(k-.5)
        return fact(2*n)/(4**n*fact(n))*math.sqrt(math.pi)
    else:
        return float(Gamma2(k))

# If the value is an :math:`n+\dfrac{1}{2} \pm \epsilon`, we'll use the special
# close-form expression. If the value is not close to :math:`n+\dfrac{1}{2}`,
# we'll use a more general approximation.

# The math.gamma() Version
# ---------------------------

# Here's a test case for the ``math.gamma()`` function. 
Example 68
Project: bruno   Author: IraKorshunova   File: utils_stats.py    MIT License 5 votes vote down vote up
def student_pdf_1d(X, phi, K, nu):
    """
    Univariate student-t density:
    output:
        the density of the given element
    input:
        x = parameter scalar
        mu = mean
        var = scale matrix
        nu = degrees of freedom
    """
    num = math.gamma((1. + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * (1. / K * (X - phi) ** 2), -(1. + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * K, 0.5)
    return num / denom 
Example 69
Project: bruno   Author: IraKorshunova   File: test_latent_dist.py    MIT License 5 votes vote down vote up
def student_pdf_1d(X, mu, var, nu):
    if nu > 50:
        return gauss_pdf_1d(X, mu, var)
    num = math.gamma((1. + nu) / 2.) * pow(
        1. + (1. / (nu - 2)) * (1. / var * (X - mu) ** 2), -(1. + nu) / 2.)
    denom = math.gamma(nu / 2.) * pow((nu - 2) * math.pi * var, 0.5)
    return num / denom 
Example 70
Project: chinook   Author: rpday   File: electron_configs.py    MIT License 5 votes vote down vote up
def Slater_exec(Z_ind,orb):

    '''

    Define an executable Slater-type orbital wavefunction which takes only
    the radius as an input argument. In this way, the usser specifies Z 
    and the orbital label string, and this generates a lambda function
    
    :math:`R(r) = (\\frac{2Z_{e}}{n_{e}})^{n_{e}} \\sqrt{\\frac{2Z_{e}}{n_{e}\\Gamma(2n_{e}+1)}} r^{n_{e}-1} e^{-\\frac{Z_{e} r}{n_{e}}}`    
 
    *args*:

        - **Z_ind**: int, atomic number
        
        - **orb**: string, 'nlxx' as per standard orbital definitions used
        throughout the library
    
    *return*:

        - executable function of position (in units of Angstrom)
    
    ***
    '''
    ne = ndic[int(orb[0])]
    xi = Z_eff(Z_ind,orb)/ne
    tmp = (2*xi)**ne*np.sqrt(2*xi/gamma(float(2*ne)+1))
    
    def lambda_gen():

        return lambda r: tmp*(r)**(ne-1)*np.exp(-xi*r)
    return lambda_gen() 
Example 71
Project: mmse-port   Author: rajivpoddar   File: mmse.py    BSD 2-Clause "Simplified" License 4 votes vote down vote up
def MMSESTSA(signal, fs, IS=0.25, W=1024, NoiseMargin=3, saved_params=None):
    SP = 0.4
    wnd = np.hamming(W)

    y = segment(signal, W, SP, wnd)
    Y = np.fft.fft(y, axis=0)
    YPhase = np.angle(Y[0:int(np.fix(len(Y)/2))+1,:])
    Y = np.abs(Y[0:int(np.fix(len(Y)/2))+1,:])
    numberOfFrames = Y.shape[1]

    NoiseLength = 9
    NoiseCounter = 0
    alpha = 0.99

    NIS = int(np.fix(((IS * fs - W) / (SP * W) + 1)))
    N = np.mean(Y[:,0:NIS].T).T
    LambdaD = np.mean((Y[:,0:NIS].T) ** 2).T

    if saved_params != None:
        NIS = 0
        N = saved_params['N']
        LambdaD = saved_params['LambdaD']
        NoiseCounter = saved_params['NoiseCounter']

    G = np.ones(N.shape)
    Gamma = G

    Gamma1p5 = math.gamma(1.5)
    X = np.zeros(Y.shape)

    for i in range(numberOfFrames):
        Y_i = Y[:,i]

        if i < NIS:
            SpeechFlag = 0
            NoiseCounter = 100
        else:
            SpeechFlag, NoiseCounter = vad(Y_i, N, NoiseCounter, NoiseMargin)

        if SpeechFlag == 0:
            N = (NoiseLength * N + Y_i) / (NoiseLength + 1)
            LambdaD = (NoiseLength * LambdaD + (Y_i ** 2)) / (1 + NoiseLength)

        gammaNew = (Y_i ** 2) / LambdaD
        xi = alpha * (G ** 2) * Gamma + (1 - alpha) * np.maximum(gammaNew - 1, 0)

        Gamma = gammaNew
        nu = Gamma * xi / (1 + xi)

        # log MMSE algo
        #G = (xi/(1 + xi)) * np.exp(0.5 * expn(1, nu))

        # MMSE STSA algo
        G = (Gamma1p5 * np.sqrt(nu)) / Gamma * np.exp(-1 * nu / 2) * ((1 + nu) * bessel(0, nu / 2) + nu * bessel(1, nu / 2))
        Indx = np.isnan(G) | np.isinf(G)
        G[Indx] = xi[Indx] / (1 + xi[Indx])

        X[:,i] = G * Y_i

    output = OverlapAdd2(X, YPhase, W, SP * W)
    return output, {'N': N, 'LambdaD': LambdaD, 'NoiseCounter': NoiseCounter} 
Example 72
Project: data_algebra   Author: WinVector   File: SQLite.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def prepare_connection(self, conn):
        # https://docs.python.org/3/library/sqlite3.html#sqlite3.Connection.create_function
        conn.create_function("is_bad", 1, _check_scalar_bad)
        # math fns
        conn.create_function("acos", 1, math.acos)
        conn.create_function("acosh", 1, math.acosh)
        conn.create_function("asin", 1, math.asin)
        conn.create_function("asinh", 1, math.asinh)
        conn.create_function("atan", 1, math.atan)
        conn.create_function("atanh", 1, math.atanh)
        conn.create_function("ceil", 1, math.ceil)
        conn.create_function("cos", 1, math.cos)
        conn.create_function("cosh", 1, math.cosh)
        conn.create_function("degrees", 1, math.degrees)
        conn.create_function("erf", 1, math.erf)
        conn.create_function("erfc", 1, math.erfc)
        conn.create_function("exp", 1, math.exp)
        conn.create_function("expm1", 1, math.expm1)
        conn.create_function("fabs", 1, math.fabs)
        conn.create_function("factorial", 1, math.factorial)
        conn.create_function("floor", 1, math.floor)
        conn.create_function("frexp", 1, math.frexp)
        conn.create_function("gamma", 1, math.gamma)
        conn.create_function("isfinite", 1, math.isfinite)
        conn.create_function("isinf", 1, math.isinf)
        conn.create_function("isnan", 1, math.isnan)
        conn.create_function("lgamma", 1, math.lgamma)
        conn.create_function("log", 1, math.log)
        conn.create_function("log10", 1, math.log10)
        conn.create_function("log1p", 1, math.log1p)
        conn.create_function("log2", 1, math.log2)
        conn.create_function("modf", 1, math.modf)
        conn.create_function("radians", 1, math.radians)
        conn.create_function("sin", 1, math.sin)
        conn.create_function("sinh", 1, math.sinh)
        conn.create_function("sqrt", 1, math.sqrt)
        conn.create_function("tan", 1, math.tan)
        conn.create_function("tanh", 1, math.tanh)
        conn.create_function("trunc", 1, math.trunc)
        conn.create_function("atan2", 2, math.atan2)
        conn.create_function("copysign", 2, math.copysign)
        conn.create_function("fmod", 2, math.fmod)
        conn.create_function("gcd", 2, math.gcd)
        conn.create_function("hypot", 2, math.hypot)
        conn.create_function("isclose", 2, math.isclose)
        conn.create_function("ldexp", 2, math.ldexp)
        conn.create_function("pow", 2, math.pow) 
Example 73
Project: data_algebra   Author: WinVector   File: SQLite.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def prepare_connection(self, conn):
        # https://docs.python.org/3/library/sqlite3.html#sqlite3.Connection.create_function
        conn.create_function("is_bad", 1, _check_scalar_bad)
        # math fns
        conn.create_function("acos", 1, math.acos)
        conn.create_function("acosh", 1, math.acosh)
        conn.create_function("asin", 1, math.asin)
        conn.create_function("asinh", 1, math.asinh)
        conn.create_function("atan", 1, math.atan)
        conn.create_function("atanh", 1, math.atanh)
        conn.create_function("ceil", 1, math.ceil)
        conn.create_function("cos", 1, math.cos)
        conn.create_function("cosh", 1, math.cosh)
        conn.create_function("degrees", 1, math.degrees)
        conn.create_function("erf", 1, math.erf)
        conn.create_function("erfc", 1, math.erfc)
        conn.create_function("exp", 1, math.exp)
        conn.create_function("expm1", 1, math.expm1)
        conn.create_function("fabs", 1, math.fabs)
        conn.create_function("factorial", 1, math.factorial)
        conn.create_function("floor", 1, math.floor)
        conn.create_function("frexp", 1, math.frexp)
        conn.create_function("gamma", 1, math.gamma)
        conn.create_function("isfinite", 1, math.isfinite)
        conn.create_function("isinf", 1, math.isinf)
        conn.create_function("isnan", 1, math.isnan)
        conn.create_function("lgamma", 1, math.lgamma)
        conn.create_function("log", 1, math.log)
        conn.create_function("log10", 1, math.log10)
        conn.create_function("log1p", 1, math.log1p)
        conn.create_function("log2", 1, math.log2)
        conn.create_function("modf", 1, math.modf)
        conn.create_function("radians", 1, math.radians)
        conn.create_function("sin", 1, math.sin)
        conn.create_function("sinh", 1, math.sinh)
        conn.create_function("sqrt", 1, math.sqrt)
        conn.create_function("tan", 1, math.tan)
        conn.create_function("tanh", 1, math.tanh)
        conn.create_function("trunc", 1, math.trunc)
        conn.create_function("atan2", 2, math.atan2)
        conn.create_function("copysign", 2, math.copysign)
        conn.create_function("fmod", 2, math.fmod)
        conn.create_function("gcd", 2, math.gcd)
        conn.create_function("hypot", 2, math.hypot)
        conn.create_function("isclose", 2, math.isclose)
        conn.create_function("ldexp", 2, math.ldexp)
        conn.create_function("pow", 2, math.pow) 
Example 74
Project: BNMTF   Author: ThomasBrouwer   File: test_bnmf_vb_optimised.py    Apache License 2.0 4 votes vote down vote up
def test_elbo():
    I,J,K = 5,3,2
    R = numpy.ones((I,J))
    M = numpy.ones((I,J))
    M[0,0], M[2,2], M[3,1] = 0, 0, 0 # size Omega = 12
    
    lambdaU = 2*numpy.ones((I,K))
    lambdaV = 3*numpy.ones((J,K))
    alpha, beta = 3, 1
    priors = { 'alpha':alpha, 'beta':beta, 'lambdaU':lambdaU, 'lambdaV':lambdaV }
    
    expU = 5*numpy.ones((I,K))
    expV = 6*numpy.ones((J,K))
    varU = 11*numpy.ones((I,K))
    varV = 12*numpy.ones((J,K))
    exptau = 8.
    explogtau = 9.
    
    muU = 14*numpy.ones((I,K))
    muV = 15*numpy.ones((J,K))
    tauU = numpy.ones((I,K))/100.
    tauV = numpy.ones((J,K))/101.
    alpha_s = 20.
    beta_s = 21.
    
    # expU * expV = [[60]]
    # (R - expU*expV)^2 = 12*59^2 = 41772
    # Var[U*V] = 12*K*((11+5^2)*(12+6^2)-5^2*6^2) = 12*2*828 = 19872
    
    # -muU*sqrt(tauU) = -14*math.sqrt(100) = -1.4
    # -muV*sqrt(tauV) = -15*math.sqrt(101) = -1.4925557853149838
    # cdf(-1.4) = 0.080756659233771066
    # cdf(-1.4925557853149838) = 0.067776752211548219
    
    ELBO = 12./2.*(explogtau - math.log(2*math.pi)) - 8./2.*(41772+19872) \
         + 5*2*(math.log(2.) - 2.*5.) + 3*2*(math.log(3.) - 3.*6.) \
         + 3.*numpy.log(1.) - numpy.log(math.gamma(3.)) + 2.*9. - 1.*8. \
         - 20.*numpy.log(21.) + numpy.log(math.gamma(20.)) - 19.*9. + 21.*8. \
         - 0.5*5*2*math.log(1./100.) + 0.5*5*2*math.log(2*math.pi) + 5*2*math.log(1.-0.080756659233771066) \
         + 0.5*5*2*1./100.*(11.+81.) \
         - 0.5*3*2*math.log(1./101.) + 0.5*3*2*math.log(2*math.pi) + 3*2*math.log(1.-0.067776752211548219) \
         + 0.5*3*2*1./101.*(12.+81.)
         
    BNMF = bnmf_vb_optimised(R,M,K,priors)
    BNMF.expU = expU
    BNMF.expV = expV
    BNMF.varU = varU
    BNMF.varV = varV
    BNMF.exptau = exptau
    BNMF.explogtau = explogtau
    BNMF.muU = muU
    BNMF.muV = muV
    BNMF.tauU = tauU
    BNMF.tauV = tauV
    BNMF.alpha_s = alpha_s
    BNMF.beta_s = beta_s
    assert BNMF.elbo() == ELBO 
Example 75
Project: quadpy   Author: nschloe   File: improve_precision.py    MIT License 4 votes vote down vote up
def stroud_1967_5():
    from quadpy.helpers import rd

    def f(x):
        degree = 5
        n = 7

        lmbda, xi, mu, gamma = x
        eta = 0
        A = 1 / 9
        B = 1 / 72
        C = B

        # data = [
        #     (B, rd(n, [(+lmbda, 1), (+xi, n - 1)])),
        #     (B, rd(n, [(-lmbda, 1), (-xi, n - 1)])),
        #     (C, rd(n, [(+mu, 2), (+gamma, n - 2)])),
        #     (C, rd(n, [(-mu, 2), (-gamma, n - 2)])),
        #     (2 * A, numpy.full((1, n), eta)),
        # ]
        # points, weights = untangle(data)
        # weights *= numpy.sqrt(numpy.pi) ** n

        data = [
            (B, rd(n, [(+lmbda, 1), (+xi, n - 1)])),
            (B, rd(n, [(-lmbda, 1), (-xi, n - 1)])),
            (C, rd(n, [(+mu, 2), (+gamma, n - 2)])),
            (C, rd(n, [(-mu, 2), (-gamma, n - 2)])),
            (2 * A, numpy.full((1, n), eta)),
        ]

        points, weights = untangle(data)
        weights *= numpy.sqrt(numpy.pi) ** n

        A = numpy.concatenate(orthopy.enr2.tree(points.T, degree, symbolic=False))

        out = numpy.dot(A, weights)
        out[0] -= numpy.sqrt(numpy.sqrt(numpy.pi)) ** n

        norm_v = numpy.sqrt(numpy.vdot(out, out))
        return norm_v

    x0 = [
        2.009_505_637_083_749e00,
        2.774_548_295_173_737e-01,
        -1.062_215_595_206_724e00,
        6.698_352_123_613_097e-01,
        # 2.009_505_6,
        # 0.277_454_83,
        # -1.062_215_60,
        # 0.669_835_21,
    ]

    out = minimize(f, x0, method="Powell", tol=1.0e-20, options={"maxiter": 20000})
    print(out.status, out.nfev, out.message, "Function value", out.fun)
    assert out.success
    print()
    for x in out.x:
        print(f"{x:.15e}")
    return 
Example 76
Project: quadpy   Author: nschloe   File: test_tools.py    MIT License 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 = quadpy.tools.coefficients_from_gauss(points, weights)
    # alpha, beta = quadpy.tools.chebyshev(moments)

    errors_alpha, errors_beta = quadpy.tools.check_coefficients(moments, alpha, beta)

    assert numpy.max(errors_alpha) > 1.0e-2
    assert numpy.max(errors_beta) > 1.0e-2 
Example 77
Project: Functional-Python-Programming-Second-Edition   Author: PacktPublishing   File: chi_sq.py    MIT License 4 votes vote down vote up
def cdf(x: float, k: int) -> float:
    """χ² cumulative distribution function.

    :param x: χ² value -- generally sum (obs[i]-exp[i])**2/exp[i]
        for parallel sequences of observed and expected values.
    :param k: degrees of freedom >= 1; generally len(data)-1

    From http://en.wikipedia.org/wiki/Chi-squared_distribution

    >>> round(cdf(0.004, 1), 2)
    0.95
    >>> round(cdf(10.83, 1), 3)
    0.001
    >>> round(cdf(3.94, 10), 2)
    0.95
    >>> round(cdf(29.59, 10), 3)
    0.001
    >>> expected=[0.95, 0.90, 0.80, 0.70, 0.50, 0.30, 0.20, 0.10, 0.05, 0.01, 0.001]
    >>> chi2= [0.004, 0.02, 0.06, 0.15, 0.46, 1.07, 1.64, 2.71, 3.84, 6.64, 10.83]
    >>> act= [round(x,3) for x in map(cdf, chi2, [1]*len(chi2))]
    >>> act
    [0.95, 0.888, 0.806, 0.699, 0.498, 0.301, 0.2, 0.1, 0.05, 0.01, 0.001]

    From http://www.itl.nist.gov/div898/handbook/prc/section4/prc45.htm

    >>> round(cdf(19.18, 6), 5)
    0.00387
    >>> round(cdf(12.5916, 6), 2)
    0.05

    From http://www.itl.nist.gov/div898/handbook/prc/section4/prc46.htm
    >>> round(cdf(12.131, 4), 5)
    0.01639
    >>> round(cdf(9.488, 4), 2)
    0.05

    """
    return 1-gamma(k/2, x/2)/Gamma_Half(k/2)

# The calcuation is 1 minus the ratio the partial
# gamma to the full gamma.

# Unit Test Cases
# ===============

# We'll use the doctest comments in each function defined above.
# Additionally, we'll use some strings with doctest test cases.