急求最小二乘法算法程序

[adashm]

我要对二元一次方程y=ax1+bx2进行a,b系数估计,谁能给出源程序,急等

[Featured]

找本《数值计算》的书，光盘里面多半有这玩意

[zengfancy]

告诉你一个算法：刚学过数值计算
y/x1 = a + b*(x2/x1) => m = a + b*n //m,n是小数数组
下面求解矩阵方程：
|sigma(1) sigma(n) | | a | | sigma(m) |
| | | | = | |
|sigma(n) sigma(n^2)| | b | | sigma(m*n) |

sigma是求和的意思，for example:
int m[5] = {1,2,3,4,5}, int n[5] = {6,7,8,9,10},
sigma(m)=15, sigma(n)=6+7+8+9+10

具体代码就不用我写了吧

[dfyang]

呵呵，刚好有一个，就贴给你了：
void CDataBaseDoc::CalculateCurveParameter(CDoubleArray* X, CDoubleArray* Y, CDoubleArray* a, int M, int N)
{
//X,Y -- X,Y两轴的坐标
//M -- 结果变量组数
//N -- 采样数目
//a -- 结果参数

int 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];
}
}
}

[fibbery]

bs楼上，把弟弟们都给惯坏了。给出思路就行了嘛！

[dfyang]

下不为例
typedef CArray<double,double> CDoubleArray;

[Mycro]

对，给出思路就行了。。

这样，不好。。。

动手能力很重要。。

[adashm]

谢谢啊,我急着用,下次一定自己编,思路我有啊,可是这个部分今天就的完成,还有别的要做呢

[adashm]

弱弱的问,CDoubleArray 可以是多维数组吗

[adashm]

我还是不明白,这个能算多元的吗

[dfyang]

这个可做曲线拟合，你将M设成2就是直线了
多维数组不行

[ldiqing]

我给你吧
#include "math.h"
//x[]是x坐标数组，y[]是y坐标数组，n是数据点个数，a[]你和结果，即多项式的系数数组；
//m为拟合阶数。
void iapcir(double x[],double y[],int n,double a[],int m,double dt[])
{ int i,j,k;
double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
for (i=0; i<=m-1; i++) a[i]=0.0;
if (m>n) m=n;
if (m>20) m=20;
z=0.0;
for (i=0; i<=n-1; i++) z=z+x[i]/(1.0*n);
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
for (i=0; i<=n-1; 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.0; t[0]=-p;
d2=0.0; c=0.0; g=0.0;
for (i=0; i<=n-1; i++)
{ q=x[i]-z-p; d2=d2+q*q;
c=c+y[i]*q;
g=g+(x[i]-z)*q*q;
}
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-1; 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.0; c=0.0; g=0.0;
for (i=0; i<=n-1; i++)
{ q=s[j];
for (k=j-1; k>=0; k--)
q=q*(x[i]-z)+s[k];
d2=d2+q*q; c=c+y[i]*q;
g=g+(x[i]-z)*q*q;
}
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];
}
}
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
for (i=0; i<=n-1; i++)
{ q=a[m-1];
for (k=m-2; k>=0; k--)
q=a[k]+q*(x[i]-z);
p=q-y[i];
if (fabs(p)>dt[2]) dt[2]=fabs(p);
dt[0]=dt[0]+p*p;
dt[1]=dt[1]+fabs(p);
}
return;
}