vendredi 1 juillet 2016

Where does "ainv" come from?

Python's random module has the following method in it:

def gammavariate(self, alpha, beta):
    """Gamma distribution.  Not the gamma function!

    Conditions on the parameters are alpha > 0 and beta > 0.

    The probability distribution function is:

                x ** (alpha - 1) * math.exp(-x / beta)
      pdf(x) =  --------------------------------------
                  math.gamma(alpha) * beta ** alpha

    """

    # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2

    # Warning: a few older sources define the gamma distribution in terms
    # of alpha > -1.0
    if alpha <= 0.0 or beta <= 0.0:
        raise ValueError('gammavariate: alpha and beta must be > 0.0')

    random = self.random
    if alpha > 1.0:

        # Uses R.C.H. Cheng, "The generation of Gamma
        # variables with non-integral shape parameters",
        # Applied Statistics, (1977), 26, No. 1, p71-74

        ainv = _sqrt(2.0 * alpha - 1.0)
        bbb = alpha - LOG4
        ccc = alpha + ainv

        while 1:
            u1 = random()
            if not 1e-7 < u1 < .9999999:
                continue
            u2 = 1.0 - random()
            v = _log(u1/(1.0-u1))/ainv
            x = alpha*_exp(v)
            z = u1*u1*u2
            r = bbb+ccc*v-x
            if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
                return x * beta

    elif alpha == 1.0:
        # expovariate(1)
        u = random()
        while u <= 1e-7:
            u = random()
        return -_log(u) * beta

    else:   # alpha is between 0 and 1 (exclusive)

        # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle

        while 1:
            u = random()
            b = (_e + alpha)/_e
            p = b*u
            if p <= 1.0:
                x = p ** (1.0/alpha)
            else:
                x = -_log((b-p)/alpha)
            u1 = random()
            if p > 1.0:
                if u1 <= x ** (alpha - 1.0):
                    break
            elif u1 <= _exp(-x):
                break
        return x * beta

While porting the code over to Java for use, my IDE has been complaining that there is a typo (misspelling) for the ainv variable name. Where does this name come from, and what is the justification behind it? If anyone can find reference material for Applied Statistics, (1977), 26, No. 1, p71-74, being able to read those pages may be helpful.




Aucun commentaire:

Enregistrer un commentaire