Computes nodes and weights for quadrature on the beta distribution.
    Default is a=b=1 which is just a uniform distribution

    NOTE: For now I am just following compecon; would be much better to
    find a different way since I don't know what they are doing.

    Parameters
    ----------
    n : scalar : int
        The number of quadrature points

    a : scalar : float, optional(default=1)
        First Beta distribution parameter

    b : scalar : float, optional(default=1)
        Second Beta distribution parameter

    Returns
    -------
    nodes : np.ndarray(dtype=float, ndim=1)
        The quadrature points

    weights : np.ndarray(dtype=float, ndim=1)
        The quadrature weights that correspond to nodes

    Notes
    -----
    Based of original function ``_qnwbeta1`` in CompEcon toolbox by
    Miranda and Fackler

    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.

def _qnwbeta1(n, a=1.0, b=1.0):
    """
    Computes nodes and weights for quadrature on the beta distribution.
    Default is a=b=1 which is just a uniform distribution

    NOTE: For now I am just following compecon; would be much better to
    find a different way since I don't know what they are doing.

    Parameters
    ----------
    n : scalar : int
        The number of quadrature points

    a : scalar : float, optional(default=1)
        First Beta distribution parameter

    b : scalar : float, optional(default=1)
        Second Beta distribution parameter

    Returns
    -------
    nodes : np.ndarray(dtype=float, ndim=1)
        The quadrature points

    weights : np.ndarray(dtype=float, ndim=1)
        The quadrature weights that correspond to nodes

    Notes
    -----
    Based of original function ``_qnwbeta1`` in CompEcon toolbox by
    Miranda and Fackler

    References
    ----------
    Miranda, Mario J, and Paul L Fackler. Applied Computational
    Economics and Finance, MIT Press, 2002.

    """
    # We subtract one and write a + 1 where we actually want a, and a
    # where we want a - 1
    a = a - 1
    b = b - 1

    maxiter = 25

    # Allocate empty space
    nodes = np.zeros(n)
    weights = np.zeros(n)

    # Find "reasonable" starting values.  Why these numbers?
    for i in range(n):
        if i == 0:
            an = a/n
            bn = b/n
            r1 = (1+a) * (2.78/(4+n*n) + .768*an/n)
            r2 = 1 + 1.48*an + .96*bn + .452*an*an + .83*an*bn
            z = 1 - r1/r2
        elif i == 1:
            r1 = (4.1+a) / ((1+a)*(1+0.156*a))
            r2 = 1 + 0.06 * (n-8) * (1+0.12*a)/n
            r3 = 1 + 0.012*b * (1+0.25*abs(a))/n
            z = z - (1-z) * r1 * r2 * r3
        elif i == 2:
            r1 = (1.67+0.28*a)/(1+0.37*a)
            r2 = 1+0.22*(n-8)/n
            r3 = 1+8*b/((6.28+b)*n*n)
            z = z-(nodes[0]-z)*r1*r2*r3
        elif i == n - 2:
            r1 = (1+0.235*b)/(0.766+0.119*b)
            r2 = 1/(1+0.639*(n-4)/(1+0.71*(n-4)))
            r3 = 1/(1+20*a/((7.5+a)*n*n))
            z = z+(z-nodes[-4])*r1*r2*r3
        elif i == n - 1:
            r1 = (1+0.37*b) / (1.67+0.28*b)
            r2 = 1 / (1+0.22*(n-8)/n)
            r3 = 1 / (1+8*a/((6.28+a)*n*n))
            z = z+(z-nodes[-3])*r1*r2*r3
        else:
            z = 3*nodes[i-1] - 3*nodes[i-2] + nodes[i-3]

        ab = a+b

        # Root finding
        its = 0
        z1 = -100
        while abs(z - z1) > 1e-10 and its < maxiter:
            temp = 2 + ab
            p1 = (a-b + temp*z)/2
            p2 = 1

            for j in range(2, n+1):
                p3 = p2
                p2 = p1
                temp = 2*j + ab
                aa = 2*j * (j+ab)*(temp-2)
                bb = (temp-1) * (a*a - b*b + temp*(temp-2) * z)
                c = 2 * (j - 1 + a) * (j - 1 + b) * temp
                p1 = (bb*p2 - c*p3)/aa

            pp = (n*(a-b-temp*z) * p1 + 2*(n+a)*(n+b)*p2)/(temp*(1 - z*z))
            z1 = z
            z = z1 - p1/pp

            if abs(z - z1) < 1e-12:
                break

            its += 1

        if its == maxiter:
            raise ValueError("Max Iteration reached.  Failed to converge")

        nodes[i] = z
        weights[i] = temp/(pp*p2)

    nodes = (1-nodes)/2
    weights = weights * math.exp(gammaln(a+n) + gammaln(b+n)
                                 - gammaln(n+1) - gammaln(n+ab+1))
    weights = weights / (2*math.exp(gammaln(a+1) + gammaln(b+1)
                         - gammaln(ab+2)))

    return nodes, weights