用最小二乘法拟合直线的问题

[dkbrain]

我用最小二乘法拟合直线y=a0+a1x,但是在计算a1、a0的时候，我用偏差为最小，然后求偏导，得出来的公式进行计算。可算法没法得到垂直于x轴的直线。
我的算法根据的是法方程：
（φ0,φ0)a0+(φ0,φ1)a1=(φ0,f)
(φ1,φ0)a0+(φ1,φ1)a1=(φ1,f)

大家给个最小二乘法直线拟合的伪代码；不要源代码的，算法也可以。

BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A)
{
 //X,Y --  X,Y两轴的坐标
 //M   --  结果变量组数
 //N   --  采样数目
 //A   --  结果参数

 register long i,j,k;
 double Z,D1,D2,C,P,G,Q;
 CDoubleArray B,T,S;
 B.SetSize(N);
 T.SetSize(N);
 S.SetSize(N);
 if(M>N)M=N;
 for(i=0;i<M;i++)
  (*A)[i]=0;
 Z=0;
 B[0]=1;
 D1=N;
 P=0;
 C=0;
 for(i=0;i<N;i++)
 {
  P=P+(*X)[i]-Z;
  C=C+(*Y)[i];
 }
 C=C/D1;
 P=P/D1;
 (*A)[0]=C*B[0];
 if(M>1)
 {
  T[1]=1;
  T[0]=-P;
  D2=0;
  C=0;
  G=0;
  for(i=0;i<N;i++)
  {
   Q=(*X)[i]-Z-P;
   D2=D2+Q*Q;
   C=(*Y)[i]*Q+C;
   G=((*X)[i]-Z)*Q*Q+G;
  }
  C=C/D2;
  P=G/D2;
  Q=D2/D1;
  D1=D2;
  (*A)[1]=C*T[1];
  (*A)[0]=C*T[0]+(*A)[0];
 }
 for(j=2;j<M;j++)
 {
  S[j]=T[j-1];
  S[j-1]=-P*T[j-1]+T[j-2];
  if(j>=3)
  {
   for(k=j-2;k>=1;k--)
    S[k]=-P*T[k]+T[k-1]-Q*B[k];
  }
  S[0]=-P*T[0]-Q*B[0];
  D2=0;
  C=0;
  G=0;
  for(i=0;i<N;i++)
  {
   Q=S[j];
   for(k=j-1;k>=0;k--)
    Q=Q*((*X)[i]-Z)+S[k];
   D2=D2+Q*Q;
   C=(*Y)[i]*Q+C;
   G=((*X)[i]-Z)*Q*Q+G;
  }
  C=C/D2;
  P=G/D2;
  Q=D2/D1;
  D1=D2;
  (*A)[j]=C*S[j];
  T[j]=S[j];
  for(k=j-1;k>=0;k--)
  {
   (*A)[k]=C*S[k]+(*A)[k];
   B[k]=T[k];
   T[k]=S[k];
  }
 }
 return TRUE;
}

谁能说下上面代码的算法啊？
还有就是这个函数如何调用的？
为什么我用下面的代码去调用的时候会出错？

CDoubleArray XXX,YYY,canshu;
	XXX.SetSize(m);
	YYY.SetSize(n);
	//canshu.SetSize(2);
	for(i=0,m=0,n=0;i<P.nHeight;i++)
	{
		for(j=0;j<P.nWidth;j++)
		{
			if(list3[i][j]==1)
			{
				XXX[m]=j;
				YYY[n]=i;
				m++;n++;
			}
		}
	}
	long M,N;
	M=2;
	N=(long)m;
	mUtil.CalculateCurveParameter(&XXX,&YYY,M,N,&canshu);

编译的时候会出现下面的提示：
ImageProcessView.obj : error LNK2019: 无法解析的外部符号 "public: int __thiscall CUtilities::CalculateCurveParameter(class CArray<double,double> *,class CArray<double,double> *,long,long,class CArray<double,double> *)" (?CalculateCurveParameter@CUtilities@@QAEHPAV?$CArray@NN@@0JJ0@Z)，该符号在函数 "public: void __thiscall CImageProcessView::OnTest(void)" (?OnTest@CImageProcessView@@QAEXXZ) 中被引用


我有定义全局typedef CArray<double,double>CDoubleArray;
也有加头文件，由于我对CArray不熟，所以不知道该怎么去定义以及设置CalculateCurveParameter函数的变量；希望各位帮帮忙看下。不甚感激

[whoo]

在 ImageProcessView.cpp 里面加上一句 ：

#include "Utilities.h"

[dkbrain]

头文件是有加的。

[dkbrain]

2楼的算法和我的算法是一样的啊，当需要拟合的线垂直于X轴的时候，是拟合不出来的。
各位有没有什么好的想法啊；给数学推导过程就可以了。

[dkbrain]

为什么没什么人回应啊？

[hhwei1985]

int main(int argc, char *argv[])

{

//输入数据

cout<<"输入数据项数：";

int num =0;

cin>> num;

valarray<double> data_x(num);

valarray<double> data_y(num);



while (num)

{

cout<< "输入第"<<num<<"项的X:";

cin>>data_x[num-1];

cout<<"输入第"<<num<<"项的Y:";

cin>> data_y[num-1];

num --;

}

cout<<"输入完毕"<<endl;



//计算系数

double A = 0.0;

double B = 0.0;

double C = 0.0;

double D = 0.0;



A = data_y.sum ();

B = data_x.sum ();

C = (data_x * data_y).sum();

D = (data_x * data_x).sum();



double k = (A * B - data_x.size () * C) / (B * B - data_x.size () * D);

double b = (B * C - A * D) / (B * B - data_x.size () * D);



//输出

cout<<"斜率k:"<<k<<endl;

cout<<"截距b:"<<b<<endl;



system("PAUSE");

return 0;

}

[Snow_Ice11111]

没仔细看你的代码，不过垂直于X轴的直线的斜率为无穷大，在用程序算时，如果没有专门针对这个处理的话，会出现除0的错误(数值溢出)，所以通常此时是得到不到正确结果的，当然也就得不到方程了。以前在学校时，我也做过类似的题目，都是要专门写代码处理这种特殊情况的。

[jian_tian_yang]

帮楼主顶一下

[jian_tian_yang]

帮楼主顶一下

[dkbrain]

各位大哥大姐帮帮忙啊，很急的我这个。

[zhangyan_wt]

关键不在拟合算法，而在于从根上拟合用的方程就是有问题的，y=ax+b，当然描述不了所有的直线(x=b)，所以一旦遇到必然抓瞎。
除非使用10楼的办法预先把这种情况排除掉，否则算法无论怎样改进也没用的。

或者还可以换一个能够描述所有直线的方程，比如：r=xcos(theta)+ysin(theta)，然后去拟合(r,theta)

[dkbrain]

我用的方法就是10楼所说的算法，只是处理点（1，1）；（1，2）；（1，3）；（1，4）；（2，2）的时候的结果就不是我所期望的一条近似垂直于x=1的直线了。
上面的点是我随便写的，直接用那个方程算过，结果居然是：y=3-0.5x;这个结果显然是不然接受的。

[neohope]

别着急，你将求解y=a0+a1x

改为求解 a2y-a1x-a0=0或者a2y-a1x=a0

然后,你函数的主要部分是不必变化的，唯一的改变是首先判断a2究竟是不是0
如果不是0，那么a2=1，调用原来的方法
如果是0，那么就a2=0,a1=1,只要求解a0就好了，一个除法而已
不同的是，要输出3个变量

第二个问题出在CDoubleArray上面，你typedef CArray <double,double>CDoubleArray; 添加了几遍？而且好象是typedef CArray <double,&double>CDoubleArray;
你要是着急，直接用数组先解决了再说

[arong1234]

没听说过用法方程求最小二乘解的不知道你怎么算的
首先：对于x=b这种类型的直线要在最小二乘之前计算的
你可以计算所有点的x座标的方差，如果小于一个给定的常数，就认为它是垂直x轴的，那么b就可以是所有x的均值

在上述条件不满足， 你才用最小二乘求解，我不知道你的算法是怎么来的，我的方法一般是：
e=sum((axi+b-yi)^2) 误差等于所有点对代进去后的误差的平方
把函数e对a,b分别求偏导数得到

de/da=2sum(xi(axi+b-yi)) = 2sum(xi*xi)a + 2sum(xi)b - 2sum(xi*yi) = 0
de/db=2sum(yi(axi+b-yi)) = 2sum(xi)a + 2Nb - 2sum(yi) = 0

因此这个就就是方程

sum(xi*xi)a + sum(xi) b = 2sum(xi*yi)
sum(xi)a + N* b = 2sum(yi)

的解，二元一次方程组，非常容易解

[dkbrain]

我的现在的处理方法是同时用：y=a0+a1x和a1y=a0+x两条直线来拟合，让后根据斜率来判断选用哪一条，虽然可以实现目标，但是这种总是不好，要写两个拟合函数，并且还是做条件判断，运算速度也不好，因为我要处理的数据点非常多，达到上几百万个点，而且要拟合很多次，所以对运算速度要求比较高，各位同仁有什么好的办法啊？

[dkbrain]

谢谢14楼的，给个建议很好，但有个小问题，因为三个参数就极值不太好处理啊。
二个参数就极值很方便，各自就骗到就OK，但是三个参数不行。

[yjingzeming]

我也找答案啊