你的近似解析表达式为y=at+bt^2+ct^2
是不是想写成为y=at+bt^2+ct^3
但是实际拟合出来的表达式为y=a[3]+a[2]t+a[1]t^2+a[0]t^3会有个常数项的。
简单的讲，所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn}，通过调整该函数中若干待定系数f(λ1, λ2,…,λ3), 使得该函数与已知点集的差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非线性拟合或者非线性回归。表达式也可以是分段函数，这种情况下叫作样条拟合。
曲线拟合：
#include <stdioh>
#include <stdlibh>
#include <malloch>
#include <mathh>
Smooth(double x,double y,double a,int n,int m,double dt1,double dt2,double dt3);
void main()
{
int i ,n ,m ;
double x,y,a,dt1,dt2,dt3,b;
n = 12;// 12个样点
m = 4; //3次多项式拟合
b = 0; //x的初值为0
/分别为x,y,a分配存贮空间/
x = (double )calloc(n,sizeof(double));
if(x == NULL)
{
printf("内存分配失败\n");
exit (0);
}
y = (double )calloc(n,sizeof(double));
if(y == NULL)
{
printf("内存分配失败\n");
exit (0);
}
a = (double )calloc(n,sizeof(double));
if(a == NULL)
{
printf("内存分配失败\n");
exit (0);
}
for(i=1;i<=n;i++)
{
x[i-1]=b+(i-1)5;
/每隔5取一个点，这样连续取12个点/
}
y[0]=0;
y[1]=127;
y[2]=216;
y[3]=286;
y[4]=344;
y[5]=387;
y[6]=415;
y[7]=437;
y[8]=451;
y[9]=458;
y[10]=402;
y[11]=464;
/x[i-1]点对应的y值是拟合已知值/
Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /调用拟合函数/
for(i=1;i<=m;i++)
printf("a[%d] = %10f\n",(i-1),a[i-1]);
printf("拟合多项式与数据点偏差的平方和为：\n");
printf("%10e\n",dt1);
printf("拟合多项式与数据点偏差的绝对值之和为：\n");
printf("%10e\n",dt2);
printf("拟合多项式与数据点偏差的绝对值最大值为：\n");
printf("%10e\n",dt3);
free(x); /释放存储空间/
free(y); /释放存储空间/
free(a); /释放存储空间/
}
Smooth(double x,double y,double a,int n,int m,double dt1,double dt2,double dt3)//(x,y,a,n,m,dt1,dt2,dt3 )
//double x; /实型一维数组，输入参数，存放节点的xi值/
//double y; /实型一维数组，输入参数，存放节点的yi值/
//double a; /双精度实型一维数组，长度为m。返回m一1次拟合多项式的m个系数/
//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 ,t ,b分配存贮空间/
s = (double )calloc(n,sizeof(double));
if(s == NULL)
{
printf("内存分配失败\n");
exit (0);
}
t = (double )calloc(n,sizeof(double));
if(t == NULL)
{
printf("内存分配失败\n");
exit (0);
}
b = (double )calloc(n,sizeof(double));
if(b == NULL)
{
printf("内存分配失败\n");
exit (0);
}
z = 0;
for(i=1;i<=n;i++)
z=z+x[i-1]/n; /z为各个x的平均值/
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]=cb[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+qq;
c=y[i-1]q+c;
g=(x[i-1]-z)qq+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[1]=ct[1];
a[0]=ct[0]+a[0];
}
for(j=3;j<=m;j++)
{
s[j-1]=t[j-2];
s[j-2]=-pt[j-2]+t[j-3];
if(j>=4)
for(k=j-2;k>=2;k--)
s[k-1]=-pt[k-1]+t[k-2]-qb[k-1];
s[0]=-pt[0]-qb[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+qq;
c=y[i-1]q+c;
g=(x[i-1]-z)qq+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[j-1]=cs[j-1];
t[j-1]=s[j-1];
for(k=j-1;k>=1;k--)
{
a[k-1]=cs[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+dtdt;
dt2=dt2+fabs(dt);
}
/释放存储空间/
free(s);
free(t);
free(b);
return(1);
}
曲线拟合的最小二乘法
如有实验数据如下
求x=55处的近似值 并画出图形
x 0 10 20 30 40 50 60 70 80 90
y 68 671 664 656 646 618 610 608 604 60
#include "stdioh"
#include "graphicsh"
double s(double x[],double y[],double x1)
{
double a,b,tmp1=0,tmp2=0,tmp3=0,tmp4=0,result;
int i=0;
for(;i<10;i++)
{
tmp1+=x[i]x[i];
tmp2+=y[i];
tmp3+=x[i];
tmp4+=x[i]y[i];
}
b=(tmp1tmp2-tmp3tmp4)/(10tmp1-tmp3tmp3);
a=(10tmp4-tmp2tmp3)/(10tmp1-tmp3tmp3);
result=ax1+b;
return result;
}
void hzbz(int x,int y)
{
int i,j;
line(0,y,1024,y);
moveto(x-200,y);
j=x-200;
for(i=0;i<100;i++)
{
if(i%5!=0)
lineto(j,y-5);
else
lineto(j,y-8);
j=j+10;
moveto(j,y);
}
line(x,768,x,0);
moveto(x,y+200);
j=y+640;
for(i=0;i<300;i++)
{
if(i%5!=0)
lineto(x+5,j);
else
lineto(x+8,j);
j=j-10;
moveto(x,j);
}
line(x,0,x-8,6);
line(x,0,x+8,6);
line(640,y,634,y-8);
line(640,y,634,y+8);
}
void main()
{
double x[]={0,10,20,30,40,50,60,70,80,90};
double y[]={68,671,664,656,646,618,610,608,604,60};
int graphdriver=DETECT,graphmode;
int i;
printf("x=55\ny=%lf\n",s(x,y,55));
getch();
initgraph(&graphdriver,&graphmode, "");
setbkcolor(BLUE);
hzbz(200,200);
moveto(200,200);
for(i=0;i<200;i++)
lineto(i+200,-s(x,y,i)+200);
getch();
closegraph();
}
//线性拟合的话好办
#include<stdioh>
#define N 10
int main( void )
{
double X[N] = {1,2,3,4,5,6,7,8,9,10};
double Y[N] = {10,9,8,7,6,5,4,3,2,1};
double xx, y, x, xy;
double k, b;
int i;
xx = y = x = xy = 0;
for( i = 0; i < N; ++i )
{
x += X[i];
xx += X[i] X[i];
xy += X[i] Y[i];
y += Y[i];
}
k = ( N xy - x y ) / ( N xx - x x );
b = ( xx y - x xy ) / ( N xx - x x );
printf( "y = %lf x + %lf\n", k, b );
return 0;
}
解决：
#include <iostreamh>
void main()
{
double X[21] = {000,0056,0112,0168,0224,0280,0336,0392,0448,0504,0560, 0616,0672,0728,0784,084,0896,0952,01008,01064,112} ;
double Y[21] = {000,166,331,496,66,822,982,114,1294,1443,1586,1722,185, 1969,2079,2179,2271,2353,2425,2487,254};
double midx;
double midy;
midx=0;
midy=0;
///求平均值
for(int i=0;i<21;++i)
{
midx+=(X[i]/210);
midy+=(Y[i]/210);
}
///求斜率
double lineK;
lineK=0;
double tempMu;
double tempZi;
tempMu=0;
tempZi=0;
for(int v=0;v<21;++v)
{
tempMu+=((X[v]-midx)(X[v]-midx));
tempZi+=((X[v]-midx)(Y[v]-midy));
}
if(tempMu==0)
{
cout <<"fail";
}
lineK=tempZi/tempMu;
////求截距
double lineB;
lineB=midy-lineKmidx;
////////////
cout <<"直线方程："<<endl;
cout <<"Y="<<lineK<<"X+"<<lineB;
cout <<endl;
///////
/////////求误差
double delta;
for(int q=0;q<21;++q)
{
delta=Y[q]-lineKX[q]-lineB;
cout <<delta <<" ";
}
}
另外，我发现倒数第二和第三组XY的误差特别大。所以我怀疑
X[21]的数据应该是：
double X[21] = {000,0056,0112,0168,0224,0280,0336,0392,0448,0504,0560, 0616,0672,0728,0784,084,0896,0952,1008,1064,112} ;
最小二乘法是统计学里的东西额。。怎么又幂函数？我猜想你是要logn复杂度的乘方快速幂吧：
//a^n mod m的logn算法
//之所以取I64d，是因为两个大数相乘，在模除之前可能会爆int，当mod超过5W的时候。
__int64 multimod(__int64 x,__int64 n,__int64 mod)
{
__int64 tmp=x,res=1LL;
while(n)
{
if(n&1LL)
{
res=tmp;
res%=mod;
}
tmp=tmp;
tmp%=mod;
n>>=1LL;
}
return res;
}
以上是我参加ACM比赛时用的模板，__int64是一种整数数据类型，范围是int的平方大小，因为乘方一般结果很大，所以比赛时常常要求运算完成后取模。如果不需要取模的话，只要将res%=mod 和 tmp%=mod两句注释掉即可
以上就是关于求C或C++语言编写的用最小二乘法进行曲线拟合全部的内容，包括:求C或C++语言编写的用最小二乘法进行曲线拟合、几个最小二乘的曲线拟合的C/C++语言代码、哪位大神可以用C语言编个程序，实现用最小二乘法，求回归线方程（暂定所用点数为10组），拜谢了！等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）