Python numbers.kept_within函数代码示例

def binprob(n, phi, integer=True):
    """
    Computes all binomial terms for n given the Bernoulli probability phi. 
    Returns the relative frequencies AND the cumulative frequencies (two lists). 
    """

    assert is_posinteger(n),  \
                      "first argument to binProb must be a positive integer!"
    assert 0.0 <= phi <= 1.0, \
                    "Bernoulli probability must be in [0.0, 1.0] in binProb!"

    farray = []
    np1    = n + 1
    for k in range(0, np1):
        x = fbincoeff(n, k, integer)
        farray.append(x)

    parray = []
    carray = []

    y = 0.0
    q = 1.0 - phi
    for k in range(0, np1):
        x  = farray[k] * phi**k * q**(n-k)
        x  = kept_within(0.0, x, 1.0)
        y  = y + x
        y  = kept_within(0.0, y, 1.0)
        parray.append(x)
        carray.append(y)

    return parray, carray

def ipareto_zero(prob, lam, xm=1.0):
    """
    The inverse of the Pareto distribution with the support shifted to [0, inf):
    f = lam * xm**lam / (x+xm)**(lam+1)
    F = 1 - [xm/(x+xm)]**lam
    x in [0, inf)
    lam > 0
    For lam < 1 all moments are infinite
    For lam < 2 all moments are infinite except for the mean
    """

    _assertprob(prob, 'ipareto_zero')
    assert lam > 0.0, "shape parameter lambda in ipareto_zero must be positive!"
    textxm1 = "left support limit parameter xm of unshifted in ipareto_zero"
    textxm2 = "distribution must not be negative in ipareto_zero!"
    assert xm  >= 0.0, textxm1 + textxm2

    q  =  1.0 - prob

    if q == 0.0: return float('inf')

    x  =  xm * (safepow(q, -1.0/lam) - 1.0)

    x  =  kept_within(0.0, x)

    return x

def ipareto(prob, lam, xm=1.0):
    """
    The inverse of the Pareto distribution: 
    f = lam * xm**lam / x**(lam+1) 
    F = 1 - (xm/x)**lam
    x in [xm, inf)
    lam > 0
    For lam < 1 all moments are infinite
    For lam < 2 all moments are infinite except for the mean
    """

    _assertprob(prob, 'ipareto')
    assert lam >  0.0, "shape parameter lambda in ipareto must be positive!"
    assert xm  >= 0.0, \
          "left support limit parameter xm must not be negative in ipareto!"

    q  =  1.0 - prob

    if q == 0.0: return float('inf')
    
    x  =  xm * safepow(q, -1.0/lam)

    x  =  kept_within(xm, x)

    return x

    def rtri_unif_tri(self, a, b, c, d):
        """
        Triangular-uniform-triangular distribution with support on [a, d] and 
        with breakpoints in b and c
                      ------
        pdf:        /        \
                   /           \
            ------               -------                                         
        """

        # Input check -----------------------
        assert a <= b and b <= c and c <= d, \
                                "break points scrambled in rtri_unif_tri!"
        # -----------------------------------


        if d == a: return a


        dcba   =  d + c - b - a
        h      =  2.0 / dcba
        first  =  0.5 * h * (b-a)
        p      =  self.runif01()
        poh    =  0.5 * p * dcba

        if p <= first:
            x  =  sqrt(2.0*(b-a)*poh) + a
        elif first < p <= first + h*(c-b):
            x  =  (c-b)*(poh-0.5*(b-a)) + b
        else:
            x  =  d - sqrt((d-c)*dcba*(1.0-p))

        x  = kept_within(a, x, d)

        return x

    def rtriang(self, left, mode, right):
        """
        Generator of triangularly distributed random numbers on [left, right] 
        with the peak of the pdf at mode. 
        """

        assert left <= mode and mode <= right, \
                                  "mode out of support range in rtriang!"

        p  =  self.runif01()

        span    =  right - left
        spanlo  =  mode  - left
        spanhi  =  right - mode
        #height  =  2.0 / span
        #surf1   =  0.5 * spanlo * height
        #surf1   =  spanlo/float(span)

        #if p <= surf1:
        if p <= spanlo/float(span):
            #x  =  sqrt(2.0*spanlo*p/height)
            x  =  sqrt(spanlo*span*p)
        else:
            #x  =  span - sqrt(2.0*spanhi*(1.0-p)/height)
            x  =  span - sqrt(spanhi*span*(1.0-p))
        x += left

        x  = kept_within(left, x, right)

        return x

def icoxian2(prob, means, probs):
    """
    The Coxian phased distribution, which is based on the exponential.
    probs is a list of probabilities for GOING ON TO THE NEXT PHASE rather 
    than reaching the absorbing state prematurely. The number of means must 
    (of course) be one more than the number of probabilities! 
    
    NB No two means[k] must be equal - if equal means are is desired, use 
    icoxian instead (slower, however). 
    """

    _assertprob(prob, 'icoxian2')
    # Everything else will be checked in dcoxian2 and ccoxian2

   # ------------------------------------
    def _fifi2fid(x):
        x      = kept_within(0.0, x)
        cdf    = ccoxian2(means, probs, x)
        pdf    = dcoxian2(means, probs, x)
        fi     = cdf - prob
        if pdf <= 0.0:
            if fi == 0.0: fi2fid = 1.0
            else:         fi2fid = MAXFLOAT
        else:
            fi2fid = fi/pdf
        return fi, fi2fid
   # ------------------------------------

    x = znewton(_fifi2fid, 0.0, 'icoxian2', tolf=SQRTMACHEPS)

    x = kept_within(0.0, x)
    
    return x

def itriang(prob, left, mode, right):
    """
    The inverse of a triangular distribution between left and right 
    having its peak at x = mode 
    """

    _assertprob(prob, 'itriang')
    # ---
    assert left <= mode and mode <= right, \
                                  "mode out of support range in itriang!"

    span    =  right - left
    spanlo  =  mode  - left
    spanhi  =  right - mode
    #height  =  2.0 / span
    #surf1   =  0.5 * spanlo * height
    #surf1   =  spanlo/float(span)

    #if prob <= surf1:
    if prob <= spanlo/float(span):
        #x  =  sqrt(2.0*spanlo*prob/height)
        x  =  sqrt(spanlo*span*prob)
    else:
        #x  =  span - sqrt(2.0*spanhi*(1.0-prob)/height)
        x  =  span - sqrt(spanhi*span*(1.0-prob))
    x += left

    x  = kept_within(left, x, right)

    return x

def ibinomial(prob, n, phi, normconst=10.0):
    """
    The binomial distribution: p(N=k) = bincoeff * phi**k * (1-phi)**(n-k), 
    n >= 1, k = 0, 1,...., n 
    """

    # Input check -----------
    _assertprob(prob, 'ibinomial')
    assert is_posinteger(n),        "n must be a positive integer in ibinomial!"
    assert 0.0 <= phi and phi <= 1.0, \
                          "success frequency out of support range in ibinomial!"
    assert normconst >= 10.0, \
           "parameter limit for normal approx. in ibinomial must not be < 10.0!"
    # -----------------------

    onemphi = 1.0 - phi

    if phi < 0.5: w = normconst * onemphi / phi
    else:         w = normconst * phi / onemphi

    if n > w:
        k = int(round(inormal(prob, n*phi, sqrt(n*phi*onemphi))))

    else:
        k   = 0
        cdf = binProb(n, phi)[1]
        while True:
            if cdf[k] <= prob:  k = k + 1
            else:               break

    k = kept_within(0, k, n)

    return k

    def rbinomial(self, n, phi):
        """
        The binomial distribution: p(N=k) = bincoeff * phi**k * (1-phi)**(n-k); 
        n >= 1;  k = 0, 1,...., n  where phi is the frequency or "Bernoulli 
        probability".
        
        Algorithm taken from ORNL-RSIC-38, Vol II (1973). 
        """

        assert is_posinteger(n), \
                          "n must be a positive integer in rbinomial!"
        assert 0.0 < phi and phi < 1.0, \
                   "frequency parameter is out of range in rbinomial!"

        normconst = 10.0
        onemphi   = 1.0 - phi
        if phi < 0.5: w = int(round(normconst * onemphi / phi))
        else:         w = int(round(normconst * phi / onemphi))

        if n > w:
            #-------------------------------------------------------
            k = int(round(self.rnormal(n*phi, sqrt(n*phi*onemphi))))

        else:
            #-------------------------------------------------------
            if phi < 0.25:
                k   = -1
                m   =  0
                phi = - safelog(onemphi)
                while m < n:
                    r  = self.rexpo(1.0)
                    j  = 1 + int(r/phi)
                    m += j
                    k += 1
                if m == n:
                    k += 1

            elif phi > 0.75:
                k   =  n + 1
                m   =  0
                phi = - safelog(phi)
                while m < n:
                    r  = self.rexpo(1.0)
                    j  = 1 + int(r/phi)
                    m += j
                    k -= 1
                if m == n:
                    k -= 1

            else: # if 0.25 <= phi and phi <= 0.75:
                k = 0
                m = 0
                while m < n:
                    r  = self.runif01()
                    if r < phi: k += 1
                    m += 1

        k = kept_within(0, k, n)

        return k

def cexppower(loc, scale, alpha, x, lngam1oalpha=False, tolf=FOURMACHEPS, itmax=128):
    """
    The exponential power distribution 
    f  =  (a/s) * exp(-abs([x-l]/s)**a) / [2*gamma(1/a)]
    F  =  1/2 * [1 + sgn(x-l) * Fgamma(1/a, abs([x-l]/s)**a)],   x in R
    s, a > 0
    where Fgamma is the gamma distribution cdf.

    NB It is possible to gain efficiency by providing the value of the 
    natural logarithm of the complete gamma function ln(gamma(1.0/alpha)) 
    as a pre-computed input (may be computed using numlib.specfunc.lngamma) 
    instead of the default 'False'.

    tolf and itmax are the numerical control parameters of cgamma.
    """

    assert scale > 0.0, "scale parameter must be a positive float in cexppower!"
    assert alpha > 0.0, "shape parameter alpha must be a positive float in cexppower!"

    if alpha == 1.0:
        return claplace(loc, scale, x)

    ainv = 1.0 / alpha
    xml = x - loc

    if not lngam1oalpha:
        lng1oa = lngamma(ainv)
    else:
        lng1oa = lngam1oalpha
    cg = cgamma(ainv, 1.0, abs(xml / scale) ** alpha, lng1oa, tolf, itmax)
    cdf = 0.5 * (fsign(xml) * cg + 1.0)

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def ctriang(left, mode, right, x):
    """
    The cdf of the triangular distribution with support 
    on [left, right] and with mode 'mode'. 
    """

    # Input check -----------------------
    assert right > left, "support range must be positive in ctriang!"
    assert left <= mode and mode <= right, "mode must be within support range in ctriang!"
    assert left <= x and x <= right, "variate must be within support range in ctriang!"
    # -----------------------------------

    spant = right - left
    spanl = mode - left
    spanr = right - mode

    if spanr == 0.0:
        cdf = (x - left) ** 2 / float((spant * spanl))

    elif spanl == 0.0:
        cdf = 1.0 - (right - x) ** 2 / float((spant * spanr))

    elif x <= mode:
        cdf = (x - left) ** 2 / float((spant * spanl))

    else:
        cdf = 1.0 - (right - x) ** 2 / float((spant * spanr))

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def cerlang(nshape, phasemean, x):
    """
    The cdf of the Erlang distribution.
    Represents the sum of nshape exponentially distributed random variables, 
    all having "phasemean" as mean
    """

    if nshape == 1:
        cdf = cexpo(phasemean, x)

    else:
        assert is_posinteger(nshape), "shape parameter must be a positive integer in cerlang!"
        assert phasemean >= 0.0, "phase mean must not be negative in cerlang!"
        assert x >= 0.0, "variate must not be negative in cerlang!"
        y = x / float(phasemean)
        cdf = 1.0
        term = 1.0
        cdf = term
        for k in range(1, nshape):
            term = term * y / k
            cdf = cdf + term

        cdf = 1.0 - exp(-y) * cdf

        cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def clevy(scale, x):
    """
    The cdf of the Levy distribution (stable distribution with 
    alpha = 1/2 and beta = 1, aka the Cournot distribution). 
    This is actually the right-skewed Levy!
    f = sqrt(s/2pi) * (1/x)**(3/2) * exp(-s/2x)
    F = erfc(sqrt(s/2x))
    
    s >= 0.0, x >= 0
    """

    assert scale >= 0.0, "scale must not be negative in clevy!"
    assert x >= 0.0, "variate must not be negative in clevy!"

    # The cdf of the Levy can be handled since it is an "incomplete gamma
    # function", but it seems to be more simple than that:

    try:
        cdf = erfc1(sqrt(0.5 * scale / x))
    except (OverflowError, ZeroDivisionError):
        return 0.0

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf

def _stable_sym_tail(alpha, x):
    """
    An asymptotic expression for the tail.
    """

    # calpha = exp(lngamma(alpha)) * sin(PIHALF*alpha) / PI
    calpha = PIINV * exp(lngamma(alpha)) * sin(PIHALF * alpha)

    try:
        cdf = calpha / x ** alpha

    except ZeroDivisionError:
        cdf = log(calpha) - alpha * log(x)
        try:
            cdf = exp(cdf)
        except OverflowError:
            cdf = 0.0

    except OverflowError:
        cdf = log(calpha) - alpha * log(x)
        try:
            cdf = exp(cdf)
        except OverflowError:
            cdf = 0.0

    cdf = 1.0 - cdf

    cdf = kept_within(0.5, cdf, 1.0)
    return cdf

def cextreme_I(type, mu, scale, x):
    """
    Extreme value distribution type I (aka the Gumbel distribution or 
    Gumbel distribution type I):
    F = exp{-exp[-(x-mu)/scale]}       (max variant)
    f = exp[-(x-mu)/scale] * exp{-exp[-(x-mu)/scale]} / scale
    F = 1 - exp{-exp[+(x-mu)/scale]}   (min variant)
    f = exp[+(x-mu)/scale] * exp{-exp[+(x-mu)/scale]} / scale

    type must be 'max' or 'min'
    scale must be > 0.0
    """

    assert scale > 0.0, "scale must be positive in cextreme_I!"

    if type == "max":
        cdf = exp(-exp(-(x - mu) / float(scale)))
    elif type == "min":
        cdf = 1.0 - exp(-exp((x - mu) / float(scale)))
    else:
        raise Error("type must be either 'max' or 'min' in cextreme_I!")

    cdf = kept_within(0.0, cdf, 1.0)

    return cdf