fitting curve for given points

587 days ago by mathmodel2017

k=4 P=[[-1,1],[1,3],[2,4],[5,20],[6,-10],[8,-2],[10,5],[8,6]] 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 8 samples are  .
using polynomial of order 4 for fitting














,





Given 8 samples are  .
using polynomial of order 4 for fitting














,





df <- data.frame(x=c(-1,1,2,5,6,8,10,8), y=c(1,3,4,20,-10,-2,5,6)) df$x2 <- df$x ** 2 df$x3 <- df$x ** 3 df$x4 <- df$x ** 4 mylm <- lm(y~x4+x3+x2+x, data=df) mylm 
       
[1]   1   1   4  25  36  64 100  64


Call:
lm(formula = y ~ x4 + x3 + x2 + x, data = df)

Coefficients:
(Intercept)           x4           x3           x2            x  
    0.61775      0.02331     -0.36490      1.22148      1.49498  
[1]   1   1   4  25  36  64 100  64


Call:
lm(formula = y ~ x4 + x3 + x2 + x, data = df)

Coefficients:
(Intercept)           x4           x3           x2            x  
    0.61775      0.02331     -0.36490      1.22148      1.49498  
a <- seq(-5,12) b <- 0.61775+1.49498*a+1.22148*a^2-0.36490*a^3+0.02331*a^4 plot(a,b)