1d quadrature weights and nodes for Gamma distributed random variable

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

    a : scalar : float, optional(default=1.0)
        Shape parameter of the gamma distribution parameter. Must be positive

    b : scalar : float, optional(default=1.0)
        Scale parameter of the gamma distribution parameter. Must be positive

    tol : scalar : float, optional(default=3e-14)
        Tolerance parameter for newton iterations for each node

    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 ``qnwgamma1`` in CompEcon toolbox by
    Miranda and Fackler

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

def _qnwgamma1(n, a=1.0, b=1.0, tol=3e-14):
    """
    1d quadrature weights and nodes for Gamma distributed random variable

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

    a : scalar : float, optional(default=1.0)
        Shape parameter of the gamma distribution parameter. Must be positive

    b : scalar : float, optional(default=1.0)
        Scale parameter of the gamma distribution parameter. Must be positive

    tol : scalar : float, optional(default=3e-14)
        Tolerance parameter for newton iterations for each node

    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 ``qnwgamma1`` in CompEcon toolbox by
    Miranda and Fackler

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

    """
    a -= 1

    maxit = 25

    factor = -math.exp(gammaln(a+n) - gammaln(n) - gammaln(a+1))
    nodes = np.zeros(n)
    weights = np.zeros(n)

    # Create nodes
    for i in range(n):
        # Reasonable starting values
        if i == 0:
            z = (1+a) * (3+0.92*a) / (1 + 2.4*n + 1.8*a)
        elif i == 1:
            z = z + (15 + 6.25*a) / (1 + 0.9*a + 2.5*n)
        else:
            j = i-1
            z = z + ((1 + 2.55*j) / (1.9*j) + 1.26*j*a / (1 + 3.5*j)) * \
                (z - nodes[j-1]) / (1 + 0.3*a)

        # root finding iterations
        its = 0
        z1 = -10000
        while abs(z - z1) > tol and its < maxit:
            p1 = 1.0
            p2 = 0.0
            for j in range(1, n+1):
                # Recurrance relation for Laguerre polynomials
                p3 = p2
                p2 = p1
                p1 = ((2*j - 1 + a - z)*p2 - (j - 1 + a)*p3) / j

            pp = (n*p1 - (n+a)*p2) / z
            z1 = z
            z = z1 - p1/pp
            its += 1

        if its == maxit:
            raise ValueError('Failure to converge')

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

    return nodes*b, weights