k=4
P=[[-1,1],[1,3],[2,4],[5,20],[6,-10],[8,-2],[10,5]]
n=len(P)
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('Given %s samples are $%s$ .'%(n,latex(PP)))
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>'%(exp,str(y[j][0]) ))
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$, $b=%s$<br><br>$A^{T}A = %s$<br><br>'%(latex(A),latex(y),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)
|
Given 7 samples are .
using polynomial of order 4 for fitting
,

Given 7 samples are .
using polynomial of order 4 for fitting
,

|