拟合算法
Posted skdev
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了拟合算法相关的知识，希望对你有一定的参考价值。
3.1 直线拟合
原理：已知一组数据Xi,Yi。Xi与Yi具有线性关系；利用最小二乘法，有如下关系
Φ(a,b)= (aXi+b-Yi)2 （1）
要使得Φ(a,b)取得最小值，那么
Φ(a,b)/ a=0, Φ(a,b)/ b=0
是极小值点的必要条件，于是(1)可化为
2Xi(aXi+b-Yi)=0
2i(aXi+b-Yi)=0
对上式求解可得：
a=[(n+1) XiYi- Xi Yi]/[(n+1) Xi2-( Xi)2]
b=[ Xi2 Yi- Xi XiYi]/[(n+1) Xi2-( Xi)2]
typedef struct tagLpoint
{
double x;
double y;
}Lpoint;
int LineMin2(HDC hDC,Lpoint *lp,int n)
{
int i;
double a,b,sumX,sumX2,sumY,sumXY;
sumX=0.0;
sumX2=0.0;
sumY=0.0;
sumXY=0.0;
for(i=0;i<n;i++)
{
sumX+=lp[i].x;
sumX2+=lp[i].x*lp[i].x;
sumY+=lp[i].y;
sumXY+=lp[i].x*lp[i].y;
}
a=(n*sumXY-sumX*sumY)/(n*sumX2-sumX*sumX);
b=(sumX2*sumY-sumX*sumXY)/(n*sumX2-sumX*sumX);
MoveToEx(hDC,0,-b/a,NULL);
LineTo(hDC,600,(600-b)/a); //600表示从高度600开始画
SelectObject(hDC,CreatePen(0,1,RGB(255,0,0)));
MoveToEx(hDC,0,b,NULL);
LineTo(hDC,700,700*a+b);//0-700的x值
return 1;
}
3.2 二次以上拟合
typedef struct tagLpoint
{
double x;
double y;
}Lpoint;
//检查方程组是否有解
int Check(double s[10][10],int R,int k)
{
int i;
i=k;
do i=i+1;
while((i<=R)&&(abs(s[i][k])<=ER));
if(i>=R)
{
MessageBox(NULL,"方程组无解或无惟一解，程序终止",NULL,NULL);
return 0;
}
return 1;
}
///把最大元素作为对角线元数
int Pivot(double s[10][10],double c[10],int R,int k)
{
int l,j;
double swap;
l=k;
for(j=k;j<R;j++)
if(abs(s[j][k])>abs(s[l][k]))
l=j;//即找出最大绝对值的下标
if(l==k)
return 0;
swap=c[k];
c[k]=c[l];
c[l]=swap;
for(j=k;j<R;j++)
{
swap=s[k][j];
s[k][j]=s[l][j];
s[l][j]=swap;
}
return 1;
}
//高斯消元法解线性方程组
int GussianSalve(double s[10][10],double c[10],int R,double a[10])
{
int i,j,k,ret;
double f;
for(k=0;k<R-1;k++)
{
ret=Check(s,R,k);
if(ret==0)
return 0;
Pivot(s,c,R,k);
for(i=k+1;i<R;i++)
{
f=s[i][k]/s[k][k];
c[i]=c[i]-c[k]*f;
for(j=k+1;j<R;j++)
s[i][j]=s[i][j]-s[k][j]*f;
}
}
a[R-1]=c[R-1]/s[R-1][R-1];
for(i=R-2;i>=0;i--)
{
f=0.0;
for(j=i+1;j<R;j++)
f=f+s[i][j]*a[j];
a[i]=(c[i]-f)/s[i][i];
}
return 1;
}
//多项式最小二乘法算法
/*名称：PolyMin
功能：对一组数据进行多项式拟合
参数：hDC绘画句柄，tempasc时个包含x,y值的数组，R表示几次如2表示一次，意思是要求ab两个参数，n表示数组个数
*/
void PolyMin(HDC hDC,Lpoint *tempasc,int R,int n)
{
int i,j,k;
double s[10][10],c[10],a[10];
double **px,temp,te;
px=new double*[n];
for(i=0;i<n;i++)
px[i]=new double[10];
for(i=0;i<R;i++)
{
c[i]=0.0;
for(j=0;j<R;j++)
s[i][j]=0.0;
}
for(i=0;i<n;i++)
{
px[i][0]=1.0;
for(j=1;j<R;j++)
px[i][j]=px[i][j-1]*tempasc[i].x;
}
for(i=0;i<R;i++)
{
for(j=0;j<n;j++)
c[i]=c[i]+tempasc[j].y*px[j][i];
for(j=0;j<R;j++)
for(k=0;k<n;k++)
s[i][j]=s[i][j]+px[k][i]*px[k][j];
}
GussianSalve(s,c,R,a);
SelectObject(hDC,CreatePen(0,1,RGB(255,0,0)));
MoveToEx(hDC,40,400-50*a[0],NULL);
for(te=0;te<31;te+=0.25)
{
temp=a[R-1];
for(j=R-1;j>0;j--)
temp=temp*te+a[j-1];
LineTo(hDC,40+te*50,400-temp*50);//从坐标点40,500处开始画
}
for(i=0;i<n;i++)
delete []px[i];
delete []px;
}
3.3 拉格朗日插值曲线
double LagBasChg(Lpoint *p,int n,double ly,int k)
{
int i;
double lf;
lf=1.0;
for(i=0;i<n;i++)
if(i!=k)
lf=lf*(ly-p[i].x)/(p[k].x-p[i].x);
return lf;
}
void LagrangeInterpChg(HDC hDC,int color,Lpoint *p,int n)
{
int i;
double ly,te;
HPEN hPen;
ly=0.0;
for(i=0;i<n;i++)
ly=ly+p[i].pasc*LagBasChg(p,n,0.0,i);
hPen=CreatePen(0,1,color);
SelectObject(hDC,hPen);
MoveToEx(hDC,40,500-ly,NULL);
for(te=0;te<=30;te+=0.25)
{
ly=0.0;
for(i=0;i<N;i++)
ly=ly+p[i].y*LagBasChg(p,n,te,i);
LineTo(hDC,40+te*20,500-ly);
}
DeleteObject(hPen);
}
以上是关于拟合算法的主要内容，如果未能解决你的问题，请参考以下文章