Computes the least squares fitting for a polynomial of order k, for a data-set in a file with columns for x and y. The condition k< n should be satisfied.

poly_lsq.py

"""

Computes the least squares fitting for a polynomial of order k.
The condition k< n should be satisfied.

"""

from scipy import zeros, dot, array, loadtxt
from scipy import shape, mean
from numpy.linalg import norm, solve

data = loadtxt('data.txt')


x = data[:,0]
y = data[:,1]

n=shape(x)[0]
k = 4

A = zeros((k + 1,k + 1))
b = zeros((k + 1))

for i in range(0,k+1):
    b[i] = dot(x**i,y)
    for j in range(0,k+1):
        r = i + j
        A[i,j] = sum(x**r)

sol = solve(A,b)
y_mean = mean(y)
f = zeros((n))
for i in range(0,k+1):
    f = f + sol[i]*x**i

print "Coefficients, from 0 to n: ", sol
print "Coefficient of determination (R2): ", 1 - sum((f -  y)**2)/sum((y - y_mean)**2)