Christoffel-Darboux recurrence relation for orthogonal polynomials using numpy

orth_poly.py

# this is a translation of orthpoly.ado from Stata 11. Their license likely applies.

def orthpoly(X, deg, weights=None):
    """
    Christoffel-Darboux recurrence relation for orthogonal polynomials.
    """
    nobs = len(X)
    orth_poly = np.ones((nobs, deg+1))
    for i in range(1, deg+1):
        t = X*orth_poly[:,i-1]**2
        b = np.average(t, weights=weights)
        t = X*t
        a = np.average(t, weights=weights)
        if i > 1:
            k = i - 2
            t = X*orth_poly[:,i-1]*orth_poly[:,i-2]
            c = np.average(t, weights=weights)
        else:
            k = 0
            c = 0
        a = 1/np.sqrt(a - b**2 - c**2)
        b *= -a
        c *= a
        orth_poly[:,i] = (a*X + b)*orth_poly[:,i-1] - c * orth_poly[:,k]
    return orth_poly[:,1:]

A one-liner free of any license issues, but not numerically sound given the Vandermonde expansion.

poly = lambda x, degree : np.linalg.qr(np.vander(x, degree + 1)[:, ::-1])[0][:, 1:]