Draft of the incomplete beta function

betainc.jl

function betainc(x, pin, qin)
    # Are these appropriate?
    eps1 = eps(0.5)
    alneps = log(eps1)
    sml = eps(0.0)
    alnsml = log(sml)

    if x < 0.0 || x > 1.0
        error("x must lie in (0, 1)")
    end

    if pin <= 0.0 || qin <= 0.0
        error("p and q must be greater than 0")
    end

    y = x
    p = pin
    q = qin
    if !(q <= p && x < 0.8 || x < 0.2)
        y = 1.0 - y
        p = qin
        q = pin
    end

    if (p + q) * y / (p + 1.0) >= eps1
        #
        # evaluate the infinite sum first.
        # term will equal y^p/beta(ps,p) * (1.0 - ps)i * y^i / fac(i)
        #
        ps = q - int(q)
        if ps == 0.0
            ps = 1.0
        end
        xb = p * log(y) - lbeta(ps, p) - log(p)
        betai = 0.0

        if xb >= alnsml
            betai = exp(xb)
            term = betai * p
            if ps != 1.0
                # 30
                n = max(alneps / log(y), 4.0e0)
                for i in 1:n
                    term = term * (i - ps) * y / i
                    betai = betai + term / (p + i)
                end
            end
        end

        # Now evaluate the finite sum, maybe
        if q > 1.0
            xb = p * log(y) + q * log(1.0 - y) - lbeta(p, q) - log(q)
            ib = max(xb / alnsml, 0.0e0)
            term = exp(xb - ib * alnsml)
            c = 1.0 / (1.0 - y)
            p1 = q * c / (p + q - 1.0)

            finsum = 0.0
            n = q
            if q == float64(n)
                n = n - 1
            end

            for i in 1:n
                if p1 <= 1.0 && term / eps1 <= finsum
                    break
                end
                term = (q - i + 1) * c * term / (p + q - i)
                if term > 1.0
                    ib = ib - 1
                end
                if term > 1.0
                    term = term * sml
                end
                if ib == 0.0
                    finsum = finsum + term
                end
            end

            betai = betai + finsum
        end

        if y != x || p != pin
            betai = 1.0 - betai
            betai = max(min(betai, 1.0), 0.0)
            return betai
        end
    end

    betai = 0.0
    xb = p * log(max(y, sml)) - log(p) - lbeta(p, q)
    if xb > alnsml && y != 0.0
        betai = exp(xb)
    end
    if y != x || p != pin
        betai = 1.0 - betai
    end
    return betai
end