Polynomial curve fitting using Vandermonde matrix - Naguib

1915 days ago by bigdata2016

Following the example in this week lecture, this is a polynomial fitting using Vandermonde matrix

 

I'm not sure why the interface is not working here.

it works in http://sage.skku.edu

 

i recommend you copy the code over there to run it properly.

 

If you have any suggestions or ideas to make the interface work here, please kindly feel free to modify the code

 

^_^

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) 
       
number of samples 

Click to the left again to hide and once more to show the dynamic interactive window