切切以待（数据处理）

[Ngadgets]

x={2001,2002,2003,2004,2005};
y={500.6,442.4,428.6,370.1,343.1};
预测当 x=2006,2007,2008时，y的值；

用数据拟合的方法，C++能实现么？
哪位学长学姐能给出方法？

[pnz_fox]

路过～

[wonxlei]

不会，帮顶 ^_^

[jixingzhong]

最小二乘法 二次拟合

[jixingzhong]

用c++写成的最小二乘法的源代码
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <process.h>
#define N 5//N个点
#define T 3 //T次拟合
#define W 1//权函数
#define PRECISION 0.00001
float pow_n(float a,int n)
{
int i;
if(n==0)
return(1);
float res=a;
for(i=1;i<n;i++)
{
res*=a;
}
return(res);
}
void mutiple(float a[][N],float b[][T+1],float c[][T+1])
{
float res=0;
int i,j,k;
for(i=0;i<T+1;i++)
for(j=0;j<T+1;j++)
{
res=0;
for(k=0;k<N;k++)
{
res+=a[i][k]*b[k][j];
c[i][j]=res;
}
}
}
void matrix_trans(float a[][T+1],float b[][N])
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<T+1;j++)
{
b[j][i]=a[i][j];
}
}
}
void init(float x_y[][2],int n)
{
int i;
printf("请输入%d个已知点：\n",N);
for(i=0;i<n;i++)
{
printf("(x%d y%d):",i,i);
scanf("%f %f",&x_y[i][0],&x_y[i][1]);
}
}
void get_A(float matrix_A[][T+1],float x_y[][2],int n)
{
int i,j;
for(i=0;i<N;i++)
{
for(j=0;j<T+1;j++)
{
matrix_A[i][j]=W*pow_n(x_y[i][0],j);
}
}
}
void print_array(float array[][T+1],int n)
{
int i,j;
for(i=0;i<n;i++)
{
for(j=0;j<T+1;j++)
{
printf("%-g",array[i][j]);
}
printf("\n");
}
}
void convert(float argu[][T+2],int n)
{
int i,j,k,p,t;
float rate,temp;
for(i=1;i<n;i++)
{
for(j=i;j<n;j++)
{
if(argu[i-1][i-1]==0)
{
for(p=i;p<n;p++)
{
if(argu[p][i-1]!=0)
break;
}
if(p==n)
{
printf("方程组无解!\n");
exit(0);
}
for(t=0;t<n+1;t++)
{
temp=argu[i-1][t];
argu[i-1][t]=argu[p][t];
argu[p][t]=temp;
}
}
rate=argu[j][i-1]/argu[i-1][i-1];
for(k=i-1;k<n+1;k++)
{
argu[j][k]-=argu[i-1][k]*rate;
if(fabs(argu[j][k])<=PRECISION)
argu[j][k]=0;
}
}
}
}
void compute(float argu[][T+2],int n,float root[])
{
int i,j;
float temp;
for(i=n-1;i>=0;i--)
{
temp=argu[i][n];
for(j=n-1;j>i;j--)
{
temp-=argu[i][j]*root[j];
}
root[i]=temp/argu[i][i];
}
}
void get_y(float trans_A[][N],float x_y[][2],float y[],int n)
{
int i,j;
float temp;
for(i=0;i<n;i++)
{
temp=0;
for(j=0;j<N;j++)
{
temp+=trans_A[i][j]*x_y[j][1];
}
y[i]=temp;
}
}
void cons_formula(float coef_A[][T+1],float y[],float coef_form[][T+2])
{
int i,j;
for(i=0;i<T+1;i++)
{
for(j=0;j<T+2;j++)
{
if(j==T+1)
coef_form[i][j]=y[i];
else
coef_form[i][j]=coef_A[i][j];
}
}
}
void print_root(float a[],int n)
{
int i,j;
printf("%d个点的%d次拟合的多项式系数为:\n",N,T);
for(i=0;i<n;i++)
{
printf("a[%d]=%g,",i+1,a[i]);
}
printf("\n");
printf("拟合曲线方程为:\ny(x)=%g",a[0]);
for(i=1;i<n;i++)
{
printf(" + %g",a[i]);
for(j=0;j<i;j++)
{
printf("*X");
}
}
printf("\n");
}
void process()
{
float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_formu[T+1][T+2],y[T+1],a[T+1];
init(x_y,N);
get_A(matrix_A,x_y,N);
printf("矩阵A为：\n");
print_array(matrix_A,N);
matrix_trans(matrix_A,trans_A);
mutiple(trans_A,matrix_A,coef_A);
printf("法矩阵为：\n");
print_array(coef_A,T+1);
get_y(trans_A,x_y,y,T+1);
cons_formula(coef_A,y,coef_formu);
convert(coef_formu,T+1);
compute(coef_formu,T+1,a);
print_root(a,T+1);
}
void main()
{
process();
}
]]>
</Content>
<PostDateTime>2007-4-19 19:23:57</PostDateTime>
</Reply>
<Reply>
<PostUserNickName></PostUserNickName>
<rank>一级(初级)</rank>
<ranknum>user1</ranknum>
<credit>100</credit>
<ReplyID>40389872</ReplyID>
<TopicID>5478010</TopicID>
<PostUserId>1526752</PostUserId>
<PostUserName>jiangxc2004</PostUserName>
<Point>0</Point>
<Content>
<![CDATA[
你可以改一下
不从终端输入,直接在程序中给出参数
请输入5个已知点：
(x0 y0):-2 -0.1
(x1 y1):-1 0.1
(x2 y2):0 0.4
(x3 y3):1 0.9
(x4 y4):2 1.6
矩阵A为：
1 -2 4 -8
1 -1 1 -1
1 0 0 0
1 1 1 1
1 2 4 8
法矩阵为：
5 0 10 0
0 10 0 34
10 0 34 0
0 34 0 130
5个点的3次拟合的多项式系数为:
a[1]=0.408571, a[2]=0.391667, a[3]=0.0857143, a[4]=0.00833333,
拟合曲线方程为:
y(x)=0.408571 + 0.391667*X + 0.0857143*X*X + 0.00833333*X*X*X

这样就可以直接调用process()函数了!
二次拟合的话就把宏 T 成2;
拟合点的数目 N 也可以修改!
也可以去到注释的部分进行返回值的调用!

[jixingzhong]

看看上面的例子，
然后自己调试一个程序出来吧

基本工作都已经好了

[G_ZhaoKun]

将x和y定义成一条曲线(y = f(x)),曲线中每两点定义成直线

[Ngadgets]

jixingzhong(瞌睡虫·星辰) 同志：
很高兴收到你的原代码，我发现对于五个数的拟合很适应。
但是，
对于9位数的拟合就不行了。
如，当我3次拟合时（T=9），
输入与输出为：

请输入9个已知点：
(x0 y0):1997
197.4
(x1 y1):1998
297.7
(x2 y2):1999
412.8
(x3 y3):2000
547.0
(x4 y4):2001
579.8
(x5 y5):2002
547.5
(x6 y6):2003
527.0
(x7 y7):2004
492.3
(x8 y8):2005
437.0
矩阵A为：
119973.98801e+0067.96405e+009
119983.992e+0067.97602e+009
119993.996e+0067.98801e+009
120004e+0068e+009
120014.004e+0068.01201e+009
120024.008e+0068.02402e+009
120034.01201e+0068.03605e+009
120044.01602e+0068.0481e+009
120054.02003e+0068.06015e+009
法矩阵为：
9180093.60361e+0077.21084e+010
180093.60361e+0077.21084e+0101.4429e+014
3.60361e+0077.21084e+0101.4429e+0142.88725e+017
7.21084e+0101.4429e+0142.88725e+0175.77745e+020
9个点的3次拟合的多项式系数为:
a[1]=-33002.3,a[2]=6.83082,a[3]=0.00356185,a[4]=6.89068e-007,
拟合曲线方程为:
y(x)=-33002.3 + 6.83082*X + 0.00356185*X*X + 6.89068e-007*X*X*X
AND
1997: 331.347
1998: 360.655
1999: 389.979
2000: 419.319
2001: 448.673
2002: 478.043
2003: 507.429
2004: 536.83
2005: 566.246
2006: 595.678
2007: 625.125
2008: 654.588
2009: 684.066
2010: 713.559
Press any key to continue

（其中输出窗口中AND后面的内容为本人在原代码中加入了for循环实现的）

[Ngadgets]

很明显，原数据形成的趋势为先升高后降低，类似为2次曲线，但当我采用2次拟合（T=9）,以达到理想拟合时。
输入与输出为：


请输入9个已知点：
(x0 y0):1997
197.4
(x1 y1):1998
297.7
(x2 y2):1999
412.8
(x3 y3):2000
547.0
(x4 y4):2001
579.8
(x5 y5):2002
547.5
(x6 y6):2003
527.0
(x7 y7):2004
492.3
(x8 y8):2005
437.0
矩阵A为：
119973.98801e+006
119983.992e+006
119993.996e+006
120004e+006
120014.004e+006
120024.008e+006
120034.01201e+006
120044.01602e+006
120054.02003e+006
法矩阵为：
9180093.60361e+007
180093.60361e+0077.21084e+010
3.60361e+0077.21084e+0101.4429e+014
9个点的2次拟合的多项式系数为:
a[1]=-39378.2,a[2]=10.5023,a[3]=0.00469825,
拟合曲线方程为:
y(x)=-39378.2 + 10.5023*X + 0.00469825*X*X
AND
1997: 331.55
1998: 360.822
1999: 390.103
2000: 419.394
2001: 448.694
2002: 478.003
2003: 507.322
2004: 536.65
2005: 565.988
2006: 595.334
2007: 624.691
2008: 654.057
2009: 683.432
2010: 712.816
Press any key to continue


很明显，不符合数据的实际变化规律。
学生我还未来得及看透你的代码奥秘，但我很希望你能作出合理的修正和解释。
不胜感激。