def set_constants( self ):
        '''
        We assume the user gives us an the parameters
        gamma and kappa such that the covariance is
        ( - gamma * Laplacian  +  kappa^2  )^{-2}.
        Then we modify these parameters such that
        
        '''
        
        self.kappa = math.sqrt( self.alpha / self.gamma )
        assert np.isreal( self.kappa )

        # We factor out gamma, so we scale kappa
        # accordingly. Later we compensate
       
        
        # Here we compensate - we now have covariance
        # [ gamma * (-Delta + kappa^2 / gamma ) ]^2
        self.sig2 = ( 
            math.gamma( self.nu )      / 
            math.gamma( 2 )            /
            (4*math.pi)**(self.dim/2.) /
            self.alpha**( self.nu )    /
            self.gamma**( self.dim/2.)
            )
        
        self.sig  = math.sqrt( self.sig2 )
        self.factor = self.sig2 * 2**(1-self.nu) / math.gamma( self.nu )
        self.ran = ( 0.0, 1.3 * self.sig2 )

    def calculate(v, x):

        # Approximation of the boys function for small x
        if x <= 25:
            i = 0
            ans = 0
            while 1 > 0:
                seq = (gamma(v + 0.5) / gamma(v + i + 1.5)) * x**i
                if seq < 1e-10:
                    break
                ans += seq
                i += 1
            ans *= (1/2) * exp(-x)
            return ans

        # Approximation of the boys function for large x
        elif x > 25:
            i = 0
            ans = 0
            while 1 > 0:
                seq = (gamma(v + 0.5) / gamma(v - i + 1.5)) * x**(-i)
                if seq < 1e-10:
                    break
                ans += seq
                i += 1
            ans *= (1/2) * exp(-x)
            ans = (gamma(v + 0.5) / (2*x**(v + 0.5))) - ans
            return ans

 def __init__(self, alpha):
     '''Creates Dirichlet distribution with parameter `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 dirichlet_pdf(alpha):
    k = len(alpha)
    gamma_sum = gamma(sum(alpha))
    product_gamma = reduce(multiply, [gamma(a) for a in alpha])
    beta = product_gamma / gamma_sum

    return lambda x: reduce(multiply, [x[i] ** (alpha[i] - 1) for i in xrange(k)])

    def uniform_normalstd_multiple_conds_with_shared_sigma(self, sigma_0, mu_0, n_subjs, seed, use_metropolis):
        """test estimation of Normal distribution std with uniform prior
            sigma_0 - the value of the std noe
            mu_0 - the value of the mu node
            use_metropolis - should it use metropolis to evaluate the sampled mean
                instead of the UniformPriorNormalstd
        """

        np.random.seed(seed)
        n_conds = len(mu_0)

        nodes, x_values = self.create_nodes_for_PriorNormalstd(n_subjs, sigma_0, mu_0, prior=pm.Uniform)
        sigma = nodes['sigma']
        mm = pm.MCMC(nodes)

        if use_metropolis:
            mm.sample(20000,5000)
        else:
            mm.use_step_method(kabuki.steps.UniformPriorNormalstd, sigma)
            mm.sample(10000)



        #calc the new distrbution
        alpha = (n_subjs*n_conds - 1) / 2.
        beta = 0
        for i_cond in range(n_conds):
            cur_x_values = x_values[i_cond*n_subjs:(i_cond+1)*n_subjs]
            beta  += sum([(x - mu_0[i_cond])**2 for x in cur_x_values]) / 2.
        true_mean = math.gamma(alpha-0.5)/math.gamma(alpha)*np.sqrt(beta)
        anal_var = beta / (alpha - 1) - true_mean**2
        true_std = np.sqrt(anal_var)

        self.assert_results(sigma, sigma_0, true_mean, true_std)
        return mm

def invgammapdf(x, alpha, beta):
	alpha = float(alpha)
	beta = float(beta)
	if not np.isscalar(x):
		return (beta**alpha / math.gamma(alpha))*np.array([(xi**(-alpha - 1))*math.exp(-beta/xi) for xi in x])
	else:
		return (beta**alpha / math.gamma(alpha))*(x**(-alpha - 1))*math.exp(-beta/x)

def variance(r0=None,L0=None,atmosphere=None):
    if atmosphere is not None:
        r0 = atmosphere.r0
        L0 = atmosphere.L0
    L0r0ratio= (L0/r0)**(5./3)
    return (24*math.gamma(6./5)/5.)**(5./6)* \
        (math.gamma(11./6)*math.gamma(5./6)/(2.*math.pi**(8./3)))*L0r0ratio

    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 get_curve(self):
        # If is already computed just return it
        if self.curve is not None:
            return self.curve

        length = 20     # maximum length in sec - see later if this need to be calculated differently
        k = self.k
        theta = self.theta
        theta_up = self.theta_up

        # This is our time vector (just the length of the gamma atom):
        t = np.linspace(0, length, length*self.fs, endpoint=False)

        # np.vectorize is not really vectorized, it's just nicer way to loop
        gamma_function_up = np.vectorize(lambda tt: 1/(math.gamma(k)*theta_up**k)*tt**(k-1)*math.exp(-tt/theta_up))
        gamma_function_down = np.vectorize(lambda tt: 1/(math.gamma(k)*theta**k)*tt**(k-1)*math.exp(-tt/theta))
        gamma_atom_up = gamma_function_up(t)
        gamma_atom_down = gamma_function_down(t)

        # stick them together : )
        gamma_atom_up = gamma_atom_up[:np.argmax(gamma_atom_up)] / np.max(gamma_atom_up)
        gamma_atom_down = gamma_atom_down[np.argmax(gamma_atom_down):] / np.max(gamma_atom_down)
        gamma_atom = np.concatenate((gamma_atom_up, gamma_atom_down))

        gamma_atom /= linalg.norm(gamma_atom)  # this preserves array and eliminates for

        return gamma_atom

def G_3(d, i_1, i_2, i_3):
    p = G_2(d, i_1, i_2) * 1./(3*(d+1) + 2*(i_1 + i_2 + i_3))
    p *= math.gamma(d + i_3 + 1)
    p /= math.gamma(i_3 + 1)
    p *= math.gamma(d + 2 + i_1 + i_2 + i_3)
    p /= math.gamma(2*d + 2 + i_1 + i_2 + i_3)
    return p

def G_2(d, i_1, i_2):
    p = G_1(d, i_1) * 1./(2*(d + 1 + i_1 + i_2))
    p *= math.gamma(d + i_2 + 1)
    p /= math.gamma(i_2 + 1)
    p *= math.gamma(i_1 + i_2 + (d+3)/2.)
    p /= math.gamma(i_1 + i_2 + (d+1)*(3./2.))
    return p

def crp_lh(theta, partition):
    n = sum([len(s) for s in partition])
    lh = 0
    lh += safety_log( math.gamma(theta)*theta**len(partition) / math.gamma(theta + n) )
    for subset in partition:
        lh += safety_log(math.gamma(len(subset)))
    return lh

def incomplete_gamma2( dA, dX ):

	if ( dA < 0 ) or ( dX < 0 ):
		return None
	if not dX:
		return 0
	xam = -dX + dA * math.log( dX )
	if ( xam > 700 ) or ( dA > 170 ):
		return 1
	if dX <= ( dA + 1 ):
		r = s = 1.0 / dA
		for k in range( 1, 61 ):
			r *= float(dX) / ( dA + k )
			s += r
			if abs( r / s ) < 1e-15:
				break
		ga = math.gamma( dA )
		gin = math.exp( xam ) * s
		return ( gin / ga )

	t0 = 0
	for k in range( 60, 0, -1 ):
		t0 = float(k - dA) / ( 1 + ( float(k) / ( dX + t0 ) ) )
	gim = math.exp( xam ) / ( dX + t0 )
	ga = math.gamma( dA )
	return ( 1 - ( gim / ga ) )

    def _calculate_exponential_params(self, moment_1=2, moment_2=4):
        """ Calculate Exponential DSD parameters.

        Calculate Exponential DSD parameters using method of moments. The choice of moments
        is given in the parameters. Uses method from [1]

        Parameters:
        moment_1: float
            First moment to use.
        moment_2: float
            Second moment to use.

        References:
        ------
        [1] Zhang, et. al., 2008, Diagnosing the Intercept Parameter for Exponential Raindrop Size
            Distribution Based on Video Disdrometer Observations: Model Development. J. Appl.
            Meteor. Climatol.,
            https://doi.org/10.1175/2008JAMC1876.1
        """

        m1 = self._calc_mth_moment(moment_1)
        m2 = self._calc_mth_moment(moment_2)

        num = m1 * gamma(moment_2 + 1)
        den = m2 * gamma(moment_1 + 1)

        Lambda = np.power(np.divide(num, den), (1 / (moment_2 - moment_1)))
        N0 = m1 * np.power(Lambda, moment_1 + 1) / gamma(moment_1 + 1)

        return Lambda, N0

def c2(psi):
    r"""Second Stumpff function.

    For positive arguments:

    .. math::

        c_2(\psi) = \frac{1 - \cos{\sqrt{\psi}}}{\psi}

    """
    eps = 1.0
    if psi > eps:
        res = (1 - np.cos(np.sqrt(psi))) / psi
    elif psi < -eps:
        res = (np.cosh(np.sqrt(-psi)) - 1) / (-psi)
    else:
        res = 1.0 / 2.0
        delta = (-psi) / gamma(2 + 2 + 1)
        k = 1
        while res + delta != res:
            res = res + delta
            k += 1
            delta = (-psi) ** k / gamma(2 * k + 2 + 1)

    return res

def c3(psi):
    r"""Third Stumpff function.

    For positive arguments:

    .. math::

        c_3(\psi) = \frac{\sqrt{\psi} - \sin{\sqrt{\psi}}}{\sqrt{\psi^3}}

    """
    eps = 1.0
    if psi > eps:
        res = (np.sqrt(psi) - np.sin(np.sqrt(psi))) / (psi * np.sqrt(psi))
    elif psi < -eps:
        res = (np.sinh(np.sqrt(-psi)) - np.sqrt(-psi)) / (-psi * np.sqrt(-psi))
    else:
        res = 1.0 / 6.0
        delta = (-psi) / gamma(2 + 3 + 1)
        k = 1
        while res + delta != res:
            res = res + delta
            k += 1
            delta = (-psi) ** k / gamma(2 * k + 3 + 1)

    return res

  def test_triangularelement( self ):
    MAXORDER = 7
    roottrans = transform.roottrans( 'test', (0,0) )
    elem = element.Element( element.SimplexReference(2), roottrans )
    F = lambda a,b: gamma(1+a)*gamma(1+b)/gamma(3+a+b)

    self._test( MAXORDER, elem, F )

  def test_tetrahedralelement( self ):
    MAXORDER = 8
    roottrans = transform.roottrans( 'test', (0,0,0) )
    elem = element.Element( element.SimplexReference(3), roottrans )
    F = lambda a,b,c: gamma(1+a)*gamma(1+b)*gamma(1+c)/gamma(4+a+b+c)

    self._test ( MAXORDER, elem, F )

 def test_gamma(self):
     import math
     assert raises(ValueError, math.gamma, 0.0)
     assert math.gamma(5.0) == 24.0
     assert math.gamma(6.0) == 120.0
     assert raises(ValueError, math.gamma, -1)
     assert math.gamma(0.5) == math.pi ** 0.5

def compute_gamma(a=0.5, h=2.0, A=math.sqrt(2), resolution=500,
                  range=[0,7]):
    """Return plot and mean/st.dev. value of the gamma density."""
    print 'range:', type(range), range
    gah = math.gamma(a + 1./h)
    mean = A*gah/math.gamma(a)
    stdev = A/math.gamma(a)*math.sqrt(
        math.gamma(a + 2./h)*math.gamma(a) - gah**2)
    x = linspace(0, range[1]*stdev, resolution+1)
    y = gamma_density(x, a, h, A)
    plt.figure()  # needed to avoid adding curves in plot
    plt.plot(x, y)
    plt.title('a=%g, h=%g, A=%g' % (a, h, A))
    if not os.path.isdir('static'):
        os.mkdir('static')
    else:
        # Remove old plot files
        for filename in glob.glob(os.path.join('static', '*.png')):
            os.remove(filename)
    # Use time since Jan 1, 1970 in filename in order make
    # a unique filename that the browser has not chached
    t = str(time.time())
    plotfile1 = os.path.join('static', 'density_%s.png' % t)
    plotfile2 = os.path.join('static', 'cumulative_%s.png' % t)
    plt.savefig(plotfile1)
    y = gamma_cumulative(x, a, h, A)
    plt.figure()
    plt.plot(x, y)
    plt.grid(True)
    plt.savefig(plotfile2)
    return plotfile1, plotfile2, mean, stdev