拟合曲线实现

使偏差平方和最小的原则称为最小二乘原则，按照该原则拟合曲线称为最小二乘法。首先我了解了线性回归相关数学原理。线性回归:线性回归假设数据集中特征与结果存在着线性关系；等式：y=mx+cy为结果，X为特征，m为系数，c为误差在数学中m为梯度c为截距。这个等式为我们假设的，我们需要找到*c使得mx+c得到的结果与真实的y误差最小，这里使用平方差来衡量估计值与真实值得误差(如果只用差值就可能会存在负数)；用于计算真实值与预测值的误差的函数称为：平方损失函数(squardlossfunction)；这里用L表示损失函数，所以有： 「Lfl=(yn-(mxr+c)):整个数据集上的平均损失为：L=77汇2G/(不”))N仁我们要求得最匹配的m与0使得L最小；数学表达式可以表示为：1丁,/argmin—汇/(三；仁加))mCnn-1最小二乘法用于求目标函数的最优值，它通过最小化误差的平方和寻找匹配项。这里将用最小二乘法用于求得线性回归的最优解。数据集有1…N个数据组成，每个数据由、构成，x表示特征，y为结果；这里将线性回归模型定义为：y(x;m.c)=mx+c平均损失函数定义有：］二, //、NM= +2mx(c-yn)+c2+加二七)n-i要求得L的最小，其关于c与m的偏导数定为0,所以求偏导数，得出后让导数等于0,并对c与m求解便能得到最小的L此时的0与m便是最匹配该模型的；关于c的偏导数，把L的等式中不包含c的项去掉得：]翼翦(c=2乂c+2-)整理式子把不和外移得到：c2整理式子把不和外移得到：c2+包含下标n的往累加对C求偏导数得:||=2c+2T冬+遮乂）关于m的偏导数：求关于m的偏导数，因此把L等式中不包含项去掉得:1K-2ynmxn+2cmx-2ynmxn+2cmxn）整理式子把不包含下标n的往累加和外移得到:1n1nNF N二对m求偏导数得：dm nJ4八'令关于C的偏导数等于0,求解：2c+2吟辱）,2c+2吟辱）,（少卜。1—m—

N从上求解得到的值可以看出，上面式子中存在两个平均值，因此该等式也可以改写成:令:歹4（»=色三）得c=y—mx令关于m的偏导数等于0,求解：关于m的偏导数依赖于c,又因为已经求得了关于c偏导数的解，因此把求关于c偏导数的解代数关于m的偏导数式子得：KT 1N 9K,=2m一之（x：）h 之二（c-乂）dm NJ、nJ'"1N 2一.=2%五g（x：）+犷gx”G-M—%）=2"4冬*）+好信七）-mx—（^xn使用元=Nxj化商表达式:也=2"4冬*）+好信七）-mx—（^xn使用元=Nxj化商表达式:也=2mdmM-强办）合并含有合并含有m的项化简:叱=2mdm叱=2mdm1二、77工（“：）一元无求解:求解:2m62m6与（6可卜2河-2-②"J-VX1N一汽（x：）-启nJ、Q-VX1N一汽（x：）-启nJ、QZByx-yx得：机x一⑺n1 >=2—xnynN J—2yxX1voidpolyfit(intn,double*x,double*y,intpoly_n,doublep[]);voidgauss_solve(intdoubleA[],doublex口，voidpolyfit(intn,double*x,double*y,intpoly_n,doublep[]);voidgauss_solve(intdoubleA[],doublex口，doubleb[]);voidpolyfit(intn,doublex[],doubley[],intpolyn,doublep[])(inti,j;double*tempx,*tempy,*sumxx,*sumxy,*ata;tempx=(double*)calloc(n,sizeof(double));sumxx=(double*)calloc((poly_n*2+l),sizeof(double));tempy=(double*)calloc(n,sizeof(double));sumxy=(double*)calloc((poly_n+l),sizeof(double));ata=(double*)calloc((poly_n+l)*(poly_n+l),sizeof(double));for(i=0;i<n;i++){tempx[i]=l;tempy[i]=y[i];}for(i=0;i<2*poly_n+l;i++)(for(sumxx[i]=0,j=0;j<n;j++){sumxx[i]+=tempx[j];tempx[j]*=x[j])))for(i=0;i<poly_n+l;i++){for(sumxy[i]=0,j=0;j<n;j++){sumxy[i]+=tempy[j];tempy[j]*=x[j];}}polyfit(n,x,y,poly_n,a)拟合y=a0+a1*x+a2*x*2+ +apol)^n*xpolynn是数据个数xy是数据值poly_n是多项式的项数返回aO,al,a2, a[poly_n],系数比项数多一(常数项)为变量申请空间for(i=0;i<poly_n+l;i++)(for(j=0;j<poly_n+l;j++){ata[i*(poly_n+l)+j]=sumxx[i+j];})gauss_solve(poly_n+l,ata,p,sumxy);free(tempx);free(sumxx);free(tempy);free(sumxy);free(ata);)voidgauss_solve(intn,doubleA[],doublex[],doubleb[])(intdoublemax;for(k=0;k<n-l;k++)(max=fabs(A[k*n+k]); //findmaxmumr=k;for(i=k+l;i<n-l;i++)(if(max<fabs(A[i*n+i]))高斯消元法计算得到n次多项式的系数n:系数的个数ata:线性矩阵sumxy:线性方程组的Y值P：返回拟合的结果max=fabs(A[i*n+i]);r=i;if(r!=k){for(i=0;i<n;i++)//changearray:A[k]&A[r](max=A[k*n+i];A[k*n+i]=A[r*n+i];A[r*n+i]=max;)max=b[k]; //changearray:b[k]&b[r]b[k]=b[r];b[r]=max;)for(i=k+l;i<n;i++)(for(j=k+l;j<n;j++)A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k];b[i]-=A[i*n+k]*b[k]/A[k*n+k];))for(i=n-l;i>=0;x[i]/=A[i*n+i],i--)(for(j=i+l,x[i]=b[i];j<n;j++)x[i]-=A[i*n+j]*x[j]；}}实验结果:区MicrosoftVisualStudio调试控制台拟合系数，按升序排列如下：012345=-18.544118=6.725933=0.236626=-0.529331=-l.450258=0.999157Matlab实现：»loadload.txt»x=load(:,1)'；»y=load(：,2),;»par