def eval_func(x, basis_func, coef):
    n = len(coef)
    assert len(basis_func) == n
    p = 0.0
    for j in range(n):
        p += coef[j] * basis_func[j](x)
    return p

def best_fit(basis_func, xi, yi):
    m = len(basis_func); N = len(xi); assert len(yi) == N
    vdm = zeros((m,m), float)
    rhs = zeros((m,), float)
    for j in range(m):
        for k in range(m):
            for i in range(N):
                vdm[j,k] += basis_func[j](xi[i]) * basis_func[k](xi[i])
            for i in range(N):
                rhs[j] += yi[i] * basis_func[j](xi[i])
    return linalg.solve(vdm,rhs)

def least_squares(m, xi, yi):
    monomial_basis = [lambda x, k=j: x**k for j in range(m)]
    coef = best_fit(monomial_basis, xi, yi)
    x_plot = linspace(min(xi), max(xi), 100)
    y_plot = [eval_func(x, monomial_basis, coef) for x in x_plot]
    return x_plot, y_plot, coef

def lsplot(x, y, order):
    data = [[],[]]
    data[0], data[1], coef = least_squares(order, x, y)
    print "coefficients are", coef
    plot(x, y, 'o')
    plot(data[0], data[1], '-')
    show()

lsplot(t, Damped_Newton(x), 5)