第一次上csdn提问，各位关照了，关于最小二乘回归法拟和线性方程

[BOBSHEN]

对具有线性关系的一组物理量 x,y,z 测出一系列测试值 xi,yi,zi(i=0..n),用最
小二乘法求"最佳"线性方程 z = A + B x + C y 的参数A, B, C。
大学好像学过，记不清了。希望大家能把解法讲讲仔细，就当我是小学生：），给出
算法更好，谢谢。
另：不知有没有通解，上面给出的是三元方程，二元可以作为它的一个特例（z=0）,那么如果有n元的解，不是大家都可以套上了,不知我的理解对不对。
刚才有个朋友帮我问了，但他草草结贴了，我还有些不明白，请各位继续关注，分我可以开贴另加！怎么最高只能放100分？

[寻开心]

不是这样的
平面方程是 ax + by + cz - d = 0;
其中 (a，b，c) 表示的是平面的法向量 d表示的是原点到平面的距离
任何一点(x0 y0 z0) 到平面的距离就是 (ax0+by0+cz0-d)/sqrt(a*a+b*b+c*c)
如果法向量是单位的话，分母就不要了。
对于一组Pi（xi,yi,zi），他们到平面的距离就是：
axi+byi+czi-d
平方距离就是：
(axi+byi+czi-d)*(axi+byi+czi-d)
所有点到这个平面的平方距离的和就是：
E（ (axi+byi+czi-d)*(axi+byi+czi-d) ） 。。。。。 （f1）
这里用E代表对i从1到n的求和符号

我们的任务是使得公式f1的数值最小，根据微积分当中的知识我们知道，在最小点的时候
E对a，b，c和d的偏导数都0，注意d也是一个变量
这样我们就可以得到方程组：
对a的偏导 E((axi+byi+czi+d)*xi) ＝ 0
对b的偏导 E((axi+byi+czi+d)*yi) ＝ 0
对c的偏导 E((axi+byi+czi+d)*zi) ＝ 0
对d的偏导 E((axi+byi+czi+d)*1) ＝ 0
变化一下，就是：
a*E(xi*xi) + b*E(xi*yi) + c*E(xi*zi) + d*E(xi) = 0
a*E(yi*xi) + b*E(yi*yi) + c*E(yi*zi) + d*E(yi) = 0
a*E(zi*xi) + b*E(zi*yi) + c*E(zi*zi) + d*E(zi) = 0
a*E(xi) + b*E(yi) + c*E(zi) + d*E(1） = 0
显然其中E(*)都是常量
这个四元一次方程组求解，就可以得到a，b，c和d的数值了。

[寻开心]

其实其它的逆合方案也都是按照这种思路作出来的
无非是采用了不同的距离公式和不同的逆合函数而已。
在这个例子当中，逆合公式就是平面方程，距离采用的是 平方距
这种逆合在求偏导的时候，由于使用的是一个二次的距离函数，所以导数好算。
如果是曲线逆合的时候这个导数就麻烦一些，主要是要解的方程组就不一定是线性的了。

[寻开心]

dcyu(Dd）说的对。
对于多维的最终化解的都是方程组，最小二乘模型下方程组都是线性的，可以使用矩阵来求解，曲线的时候就不同了。

[BOBSHEN]

非常感谢 happy__888 我先慢慢消化，有问题一回再问。
另：我开给分贴可不可以？CSDN我规矩我不大董，我觉得分数限制的有点太低了。

[BOBSHEN]

问题又来了：
抛开前面关于原理的讨论，是不是对于一组数据，我只要套用上面的公式构造方程组求解就可以了？对于多元的情况，只要他是线性的就都可以套用上面的公式，只是多了几个未知数和方程而已？

[dcyu]

看来没必要去扛书了，有人都答了，偶记得是列出一个数据矩阵3*n的，然后和转置相乘得到一个矩阵3*3的，再...，后面的忘了。

你再用同样的问题发一次贴，然后叫别人再回答一遍就行了。

[dyf2001]

上面说得很好，不过我感觉有个问题就是这是在理想状况下，数据都在假定平面附近，但如果有特殊点距离理想平面很远应该剔除该点，不然结果不可能是最优解

[BOBSHEN]

TO dcyu(Dd):
谢谢，你再帮忙翻翻书也不错啊，就当温温课了。我也想知道用矩阵怎么解，是不是最终也回到
楼上总结的方程组上去了，换一种思路解决问题也不错。分不成问题

另：happy__888请取接分！

http://expert.csdn.net/Expert/topic/1965/1965060.xml?temp=.6579248
http://expert.csdn.net/Expert/topic/1965/1965047.xml?temp=.7074243

[BOBSHEN]

TO dyf2001:
老兄说的没错，是个应该注意的问题，不过我的数据再采集之前已经用上限、下限以及跳变率校验过了，所以异常数据大都被剔除掉了。

[寻开心]

异常大的仅仅是一个方面
你还可以通过多次逆合方法来提高精度
办法是先用大量数据做第一次逆合
然后提出距离上次逆合结果距离比较大的点（自己定义一个最大距离的数值），然后用剩余的数据重新逆合。
这样可以提高一些精度。
对于最小二乘逆合而言，做两三次这样的操作就可以了。

[寻开心]

哈。
要是再多一些这样的肯给分的好楼主，俺可就可以当猩猩喽！

[BOBSHEN]

我kao!!!放分按钮在什么地方！！！烧香连庙门都找不到，太菜！

to happy__888:
我上面的回复请你再看看
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
问题又来了：
抛开前面关于原理的讨论，是不是对于一组数据，我只要套用上面的公式构造方程组求解就可以了？对于多元的情况，只要他是线性的就都可以套用上面的公式，只是多了几个未知数和方程而已？

[寻开心]

晕！
为了要当猩猩，只好再详细一点回答啦。
a*E(xi*xi) + b*E(xi*yi) + c*E(xi*zi) + d*E(xi) = 0
a*E(yi*xi) + b*E(yi*yi) + c*E(yi*zi) + d*E(yi) = 0
a*E(zi*xi) + b*E(zi*yi) + c*E(zi*zi) + d*E(zi) = 0
a*E(xi) + b*E(yi) + c*E(zi) + d*E(1） = 0
没有理解“E(*)都是常量”这句话吧。
对于一组数据（x1, y1, z1）, (x2, y2, z2), ... , （xn, yn,zn）
E(xi*xi) ＝ x1*x1 + x2*x2 + .... + xn*xn
E(xi*yi) = x1*y1 + x2*y2 + .... + xn*yn
其它的类似
E(1) = 1 + 1 + ... + 1 = n
这样对于给定的一组数据，上面的四元一次方程组就出来了，直接求解方程组，就可以得到了。

多元的就是解更多参数的方程组。

不好意思，要分心切：
在数据结构和算法栏目当中，点击你发的帖子后面的管理按钮，就可以给分喽。

[BOBSHEN]

>>这样对于给定的一组数据，上面的四元一次方程组就出来了，直接求解方程组，就可以得到了。

其实我要的就是这句话！
对下面的过程等有空的时候在补补微积分课吧
>> 对a的偏导 E((axi+byi+czi+d)*xi) ＝ 0
>>转化成
>> a*E(xi*xi) + b*E(xi*yi) + c*E(xi*zi) + d*E(xi) = 0

搞了半天，一个放出去了，下面这个死活放不出去，第一次提示只有"只有非0分贴。。。"第二次连给分按钮都没了。（悄悄说：csdn的页面作的台繁杂了）。
君子爱分，取之有道，怎么接这个贴？
http://expert.csdn.net/expert/Topicview2.asp?id=1965047

[寻开心]

>> 对a的偏导 E((axi+byi+czi+d)*xi) ＝ 0
>>转化成
>> a*E(xi*xi) + b*E(xi*yi) + c*E(xi*zi) + d*E(xi) = 0

就是提取公因子啊。

第一个公式展开就是：
(a*x1 + b*y1 + c*z1 +d)*x1 + (a*x2 + b*y2 + c*z2 +d)*x2 + ...
+ (a*xn + b*yn + c*zn +d)*xn
=
a * x1*x1 + b*y1*x1 + c*z1*x1 + d*x1 +
a * x2*x2 + b*y2*x2 + c*z2*x2 + d*x2 +
...
a * xn*xn + b*yn*xn + c*zn*xn + d*xn
=
a * (x1*x1 + x2*x2 + ... + xn*xn) +
b * (y1*x1 + y2*x2 + ... + yn*xn) +
c * (z1*x1 + z2*x2 + ... + zn*xn) +
d * (x1 + x2 + ... + xn )
=
a * E(xi*xi) + b*E(yi*xi) + c*E(zi*xi) + d*E(xi)

其实你记住 求和符号和乘法符号是可以交换的就可以了

[寻开心]

刚才回答了一半，没写完呢。

最小二乘的原理就是先虚拟一个最优解，然后根据点到直线或者平面的表达式，计算得到测试点到虚拟的平面的平方距离。把这些平方距离取和，得到一个表达式，也就是这些离散点到这个虚拟的直线或者是平面的平方距离的和D(XI,YI,ZI)。
接下来的任务是如何具体落实这个直线或者平面，使得这个平方和最小。
根据微积分当中的知识，我们对这个距离函数取导数，导数数值为0的就代表了极值点，在这个条件下这个极值点肯定是最小值而不是最大值，这个导数应该是偏导数，对直线或者是平面的生成参数的各个偏导数值都是0。

根据这个原理我们就可以自己推导出最小二乘公式的。

[BOBSHEN]

呵呵，csdn好像有点欺生啊，怎么没人理？
分太少吗？我说了可以开贴再加的！

[HUNTON]

对两个物理量的还可以有公式，三个的应该就没有直接的公式了，就按照happy__888说的通法做吧

[寻开心]

都个量的要解多元方程组，比较麻烦，没有通解是自然的。

[boylez]

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

double getd(double [3][4],int);
main(){
int i,n;
int xarray[10],yarray[10];
int sumx,sumy,sumxy,sumx2,sumx2y,sumx3,sumx4;
double matrix[3][4];
double a,b,c,x,delta;
double d[4];
printf("\npress any key to start.\n");
for(i=0;i<10;i++){
printf("Input x&y(0&0 end):");
scanf("%d%d",&xarray[i],&yarray[i]);
if(xarray[i]==0&&yarray[i]==0) break;
}
sumx=sumy=sumxy=sumx2=sumx2y=sumx3=sumx4=0;
for(n=0;n<i;n++){
sumx+=xarray[n];
sumy+=yarray[n];
sumxy+=xarray[n]*yarray[n];
sumx2+=xarray[n]*xarray[n];
sumx2y+=xarray[n]*xarray[n]*yarray[n];
sumx3+=xarray[n]*xarray[n]*xarray[n];
sumx4+=xarray[n]*xarray[n]*xarray[n]*xarray[n];
}
matrix[0][0]=i;
matrix[0][1]=matrix[1][0]=sumx;
matrix[0][2]=matrix[2][0]=matrix[1][1]=sumx2;
matrix[1][2]=matrix[2][1]=sumx3;
matrix[2][2]=sumx4;
matrix[0][3]=sumy;
matrix[1][3]=sumxy;
matrix[2][3]=sumx2y;
d[0]=getd(matrix,3);
for(n=1;n<4;n++){
d[n]=getd(matrix,n-1);
}
c=d[1]/d[0];
b=-d[2]/d[0];
a=d[3]/d[0];
printf("\n\n%fx^2",a);
if(b>0) printf("+");
if(b) printf("%fx",b);
if(b&&c>0) printf("+");
if(c) printf("+%f",c);
printf("=0\n\n");
if(b*b-4*a*c>=0){
delta=pow(b*b-4*a*c,0.5);
x=(b*(-1)+delta)/(2*a);
printf("\nx=%f",x);
x=(b*(-1)-delta)/(2*a);
printf("\nx=%f",x);
}
else{
printf("error!!!\nb^2-4ac<0\n");
exit(0);
}
printf("\nend!press any key to end.");
getch();
return 0;
}
double getd(double origrinalmatrix[3][4],int n)
{
int i,j,k;
double matrix[3][3];
double result;
for(i=0;i<3;i++){
for(j=0,k=0;j<3&&k<4;j++,k++){
if(k==n){
j--;
continue;
}
else{
matrix[i][j]=origrinalmatrix[i][k];
}
}
}
result=(matrix[0][0]*matrix[1][1]*matrix[2][2]-matrix[0][0]*matrix[2][1]*matrix[1][2]-matrix[0][1]*matrix[1][0]*matrix[2][2]+matrix[0][1]*matrix[1][2]*matrix[2][0]+matrix[0][2]*matrix[1][0]*matrix[2][1]-matrix[0][2]*matrix[1][1]*matrix[2][0]);
return result;
}