一元多项式回归

设拟合多项式为φm(x)

残差平方和

增广矩阵法

由于需要构造系数为1，将进行大量浮点运算，况且矩阵的元素肯能很大，会出现大吃小的可能，会引入大量误差

逆矩阵法

依旧计算量相当大，浮点运算误差影响更大

R代码

三次多项式拟合

x=matrix(c(5.0, 15.0, 55.0, 225.0 , 15.0, 55.0, 225.0, 979.0 , 55.0, 225.0, 979.0, 4425.0, 225.0, 979.0, 4425.0, 20515.0 ), nr=4,nc=4) b=matrix(c(15.0, 56.0, 231.0 ,1007.0 ),nr=4,nc=1) c=solve(x)%*%b a=matrix(c(1.5,-1,0.5),nr=3,nc=1) > c [,1] [1,] 2. [2,] -2. [3,] 1. [4,] -0. 

针对对称矩阵有LDLT分解法（LT为L的转置）

java代码

public class CurveFit { 
    public static void main(String[] args) { 
    double[] x=new double[] { 
   1,2,3,4,5}; double[] y=new double[] { 
   1,1.5,3,4.5,5}; show(curveFit(x,y,3)); } public static double [] curveFit(double[] x,double[] y,int pow) { 
    double[][] matrix=createMatrix(x,y,pow); LDLTDCMP(matrix); double [] xy=iteration(y,x,pow+1); LDLTBKSB(matrix,xy); return xy; } public static double[][] createMatrix(double[] x,double[] y,int pow) { 
    double [] xx=iteration(x,x,2*pow); double[][] matrix=new double[pow+1][pow+1]; matrix[0][0]=x.length; for(int i=0;i<pow+1;i++) { 
    for(int k=0;k<pow+1;k++) { 
    if(i+k>0) { 
    matrix[i][k]=xx[i+k-1]; } } } return matrix; } public static double[] iteration(double[] start,double [] value,int count) { 
    double[] itr = null; double[] ret = new double [count]; for(int i=0;i<count;i++) { 
    if(i==0) { 
    itr=start.clone(); }else { 
    itr=mul(itr,value); } ret[i]=sum(itr); } return ret; } public static double[] mul(double[] a,double[] b) { 
    double[] c=a.clone(); if(a.length==b.length) { 
    for(int i=0;i<a.length;i++) { 
    c[i]*=b[i]; } } return c; } public static double sum(double[] a) { 
    double ret=0; for(int i=0;i<a.length;i++) { 
    ret+=a[i]; } return ret; } public static void LDLTDCMP ( double [][]a) { 
    int k=0,m=0,i=0,n=a.length; for (k = 0; k < n; k++) { 
    for (m = 0; m < k; m++) { 
    a[k][k] = a[k][k] - a[m][k] * a[k][m]; } if (a[k][k] == 0) { 
    return; } for (i = k + 1; i < n; i++) { 
    for (m = 0; m < k; m++) a[k][i] = a[k][i] - a[m][i] * a[k][m]; a[i][k] = a[k][i] / a[k][k]; } } } public static void LDLTBKSB ( double[][]a, double []b) { 
    int i; int k; int n=a.length; for (i = 0; i < n; i++){ 
    for (k = 0; k < i; k++) b[i] = b[i] - a[i][k] * b[k]; } for (i = n - 1; i >= 0; i--){ 
    b[i] = b[i] / a[i][i]; for (k = i + 1; k < n; k++) b[i] = b[i] - a[k][i] * b[k]; } } public static void show(double[][] array) { 
    for(int i=0;i<array.length;i++) { 
    double[] itme=array[i]; for(double out:itme) { 
    System.out.print(out+"\t"); } System.out.println(); }System.out.println(); }; public static void show(double[] array) { 
    for(double out:array) { 
    System.out.print(out+"\t"); } System.out.println(); System.out.println(); } } 

2.0047 -2.3396 1.00235 -0.66693

R语言多项式回归

x=c(1,2,3,4,5) y = c(1,1.5,3,4.5,5) xy=data.frame(x=x,y=y) #model <- lm(xy$y ~ poly(x,3)) model <- lm(xy$y ~ x + I(x^2) + I(x^3)) 

回归模型

> model Call: lm(formula = xy$y ~ x + I(x^2) + I(x^3)) Coefficients: (Intercept) x I(x^2) I(x^3) 2.5000 -2.8333 1.5000 -0.1667 

模型摘要

> summary(model) Call: lm(formula = xy$y ~ x + I(x^2) + I(x^3)) Residuals: 1 2 3 4 5 3.052e-16 -1.221e-15 1.831e-15 -1.221e-15 3.052e-16 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.500e+00 1.256e-14 1.990e+14 3.20e-15 * x -2.833e+00 1.642e-14 -1.726e+14 3.69e-15 * I(x^2) 1.500e+00 6.094e-15 2.461e+14 2.59e-15 * I(x^3) -1.667e-01 6.729e-16 -2.477e+14 2.57e-15 * --- Signif. codes: 0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.554e-15 on 1 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 6.39e+29 on 3 and 1 DF, p-value: 9.196e-16 

p值9.196e-16接近于0说明完美拟合。

拟合图像

x=c(1,2,3,4,5) y = c(1,1.5,3,4.5,5) f1=function(x){ -0.3+1.1*x } f3=function(x){ 2.5-2.83*x+1.5*x^2-0.167*x^3 } plot(x,y) curve(f1,from=0,to=6,col="blue",add=T) curve(f3,from=0,to=6,col="red",add=T) text(4,3.5, "y = -0.3+1.1*x",col="blue") text(3,1.5, "y = 2.5-2.83x+1.5x^2-0.167X^3",col="red") 

过拟合问题