Polynomial curve fitting using Vandermonde matrix - Naguib

import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D from numpy import argmax,argmin from sage.plot.scatter_plot import ScatterPlot global n global k global alpha @interact def f(n=slider(3,20, 1, 7, label="number of samples")): P=[(round(random()*100),round(random()*100)) for _ in range(n)] x=matrix(RDF,n,1) y=matrix(RDF,n,1) PP=matrix(RDF,n,2) for i in range(n): x[i]=P[i][0] y[i]=P[i][1] PP[i,0]=P[i][0] PP[i,1]=P[i][1] html('%s samples are $%s$ created by Sage randomly.'%(n,latex(PP))) @interact def f(k=slider(0,10, 1, 7, label="polynomial order used in fitting")): if k > n-2 : html('Please, choose polynomial order to be considerably less than number of samples to avoid overfitting!') k=n-2 html('using polynomial of order %s for fitting<br><br>'%(k)) ################################################## fitting exp='a_0' for i in range(1,k+1): exp = exp + ' + a_{' + str(i) + '}x^{' + str(i) + '}' html('$y = %s$<br><br>'%(exp)) for j in range(n): exp='a_0' for i in range(1,k+1): exp = exp + ' + ' + str(x[j][0]^i) + 'a_{' + str(i) + '}' html('$%s = %s$<br>'%(str(y[j][0]), exp)) html('<br>') A=matrix(RDF,n,k+1) for j in range(n): for i in range(k+1): A[j,i]=x[j][0]^i html('$Ax=b$<br><br>$A = %s$<br><br>$A^{T}A = %s$<br><br>'%(latex(A),latex(A.transpose()*A))) alpha =(A.transpose()*A).inverse() * A.transpose() * y html('$x = %s$<br><br>'%(latex(alpha))) for i in range(k+1): html('$a_{' + str(i) + '} = ' + str(alpha[i][0])) ################################################## plotting def polynomial(x): ret=0 for i in range(k+1): ret = ret + alpha[i][0]*x^i return ret plotstart=min(x)[0]-5 plotend=max(x)[0]+5 curve=plot(polynomial,(plotstart, plotend)) show(list_plot(P,color='red') + curve)