曲线拟合问题：最小二乘法算法

[duangexin521]

现在想对一个曲线进行拟合，在网上看到一些算法。试验结果参数系数是算出来了。但是返回把结果放进去验证 得出来的结果相差甚远。请问哪个高手有相关资料。请指点一下。
void run()
{
int i,n,m;

double Yvalue[21] = {0.00,0.056,0.112,0.168,0.224,0.280,0.336,0.392,0.448,0.504,0.560,
0.616,0.672,0.728,0.784,0.84,0.896,0.952,1.008,1.064,1.12};

double CRAvalue[21] = {0.00,1.66,3.31,4.96,6.6,8.22,9.82,11.4,12.94,14.43,15.86,17.22,18.5,
19.69,20.79,21.79,22.71,23.53,24.25,24.87,25.4};

double *x,*y,*a,dt1,dt2,dt3;

double aaa[6];

n = 21;
m = 6;

/*分别为x,y,a分配存贮空间*/

x = new double[n];
y = new double[n];
a = new double[n];

/*取21个点来进行拟合*/
for(i=0;i<n;i++)
{
x[i] = Yvalue[i];
y[i] = CRAvalue[i];

}

/*调用拟合函数,求出各多项式的参数值,并求出拟合多项式与数据点偏差的平方和dt1,绝对值之和dt2,绝对值最大值dt3*/
Smooth(x,y,a,n,m,&dt1,&dt2,&dt3);

for(i=0;i<6;i++)
{
aaa[i] = *a;
++a;
}

double aa = 0.056;

double bb;

bb = a[0]+a[1]*aa+a[2]*aa*aa+a[3]*aa*aa*aa+a[4]*aa*aa*aa*aa+a[5]*aa*aa*aa*aa*aa;


if (!x)
{
delete []x;
x = NULL;
}

if (!y)
{
delete []y;
y = NULL;
}

if(!a)
{
delete []a;
a = NULL;
}

}


/*-------------------------------------------------------------------------

double *x; / *实型一维数组，输入参数，存放节点的xi值* /
double *y; / *实型一维数组，输入参数，存放节点的yi值* /
double *a; / *双精度实型一维数组，长度为m。返回m一1次拟合多项式的m个系数* /
int n; / *整型变量，输入参数，给定数据点的个数* /
int m; / *整型变量，输入参数，拟合多项式的项数* /
double *dt1; / *实型变量，输出参数，拟合多项式与数据点偏差的平方和* /
double *dt2; / *实型变量，输出参数，拟合多项式与数据点偏差的绝对值之和* /
double *dt3; / *实型变量，输出参数，拟合多项式与数据点偏差的绝对值最大值* /

---------------------------------------------------------------------------*/
Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3 )
{
int i,j,k;
double *s,*t,*b,z,d1,p,c,d2,g,q,dt;


s = new double[n];
t = new double[n];
b = new double[n];


z = 0;
for(i=1;i<=n;i++)
z=z+x[i-1]/n;

b[0]=1;
d1=n;
p=0;
c=0;
for(i=1;i<=n;i++)
{
p=p+x[i-1]-z;
c=c+y[i-1];
}

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=1;i<=n;i++)
{
q=x[i-1]-z-p;
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-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=3;j<=m;j++)
{
s[j-1]=t[j-2];
s[j-2]=-p*t[j-2]+t[j-3];
if(j>=4)
for(k=j-2;k>=2;k--)
s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];
s[0]=-p*t[0]-q*b[0];
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=s[j-1];
for(k=j-1;k>=1;k--)
q=q*(x[i-1]-z)+s[k-1];
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;

a[j-1]=c*s[j-1];
t[j-1]=s[j-1];
for(k=j-1;k>=1;k--)
{
a[k-1]=c*s[k-1]+a[k-1];
b[k-1]=t[k-1];
t[k-1]=s[k-1];
}

}
*dt1=0;
*dt2=0;
*dt3=0;
for(i=1;i<=n;i++)
{
q=a[m-1];
for(k=m-1;k>=1;k--)
q=q*(x[i-1]-z)+a[k-1];
dt=q-y[i-1];
if(fabs(dt)>*dt3)
*dt3=fabs(dt);
*dt1=*dt1+dt*dt;
*dt2=*dt2+fabs(dt);
}


if(!s)
{
delete []s;
s = NULL;
}

if(!t)
{
delete []t;
t = NULL;
}

if(!b)
{
delete []b;
b = NULL;
}

//free(s);
//free(t);
//free(b);
return(1);
}

[duangexin521]

自己已经解决！

[实达诚实]

mark

[会思考的草]

我比较喜欢用Chebyshev多项式去拟合，能保证精度，多项式的次数也不会太高。
另外用Matlab做这个快多了。

[hityrj]

我想用最小二乘法拟合出我想要的公式类型，不是普通的多项式，比如说Y=（B-A）/(pow(C/X,D)+1),怎么才能算出这些系数呢……请大家多多指教……我刚学C++

[hityrj]

我想用最小二乘法拟合出我想要的公式类型，不是普通的多项式，比如说Y=（B-A）/(pow(C/X,D)+1),怎么才能算出这些系数呢……请大家多多指教……我刚学C++

[luxiankao]

楼主您是如何解决的方便把代码贴出来吗，谢谢！