急求最小二乘法算法程序

Question

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

Answer 1

告诉你一个算法：刚学过数值计算
   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具体代码就不用我写了吧

Answer 2

呵呵，刚好有一个，就贴给你了：
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];
}
}
}

Answer 3

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

Answer 4

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

Answer 5

我给你吧
#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;
  }