[数值算法]曲线线拟合的最小二乘法
[数值算法]曲线线拟合的最小二乘法.
给出一些离散的点通过一条近似的曲线来表示这些点所代表的函数,这就是曲线拟合的目的.最小二乘法是和种简便而又实用的曲线拟合方法.
它除了能直接拟合形如y=a*x+b形的函数外.
还能对
1/y=a+b/x,(a>0)
y=a*exp(b*x),(a>0)
a=a*exp(b/x)(x>0,a>0),
y=a+b*lg(x),
y=a*x^b(a>0),
y=1/(a+b*exp(-x))(a>0);
等式子变换成线性表达式后进行拟合.
相关的理论请参考相关的数值算法的书籍，我这里只给出关键的函数及主程序段，其余相关的细节就不再一一罗列了.
void min2ByFixMethod1(Type* xArr,Type* yArr,int dataNum,FILE* outputFile)
{
Type* x_xArr,* x_yArr;/*make up data*/
Type c[2][2]={{0,0},{0,0}};/*medium data variable*/
Type d[2]={0,0};
Type tween;/*fatcor's mother*/
Type a,b;/*argu ...*/
assertF(xArr!=NULL,"in min2ByFixMethod1,xArr is NULL");
assertF(yArr!=NULL,"in min2ByFixMethod1,yArr is NULL");
assertF(outputFile!=NULL,"in min2ByFixMethod1,outputFile is NULL");
/*Memory apply*/
x_xArr=(Type*)malloc(sizeof(Type)*dataNum);
x_yArr=(Type*)malloc(sizeof(Type)*dataNum);
assertF(x_xArr!=NULL,"in min2ByFixMethod1,x_xArr is NULL");
assertF(x_yArr!=NULL,"in min2ByFixMethod1,x_yArr is NULL");
/*Other aid data make*/
arrEachBy(xArr,xArr,&x_xArr,dataNum);
arrEachBy(xArr,yArr,&x_yArr,dataNum);
c[0][0]=sumArr(x_xArr,dataNum);
c[0][1]=sumArr(xArr,dataNum);
c[1][0]=sumArr(xArr,dataNum);
c[1][1]=(Type)dataNum;
d[0]=sumArr(x_yArr,dataNum);
d[1]=sumArr(yArr,dataNum);
tween=c[0][0]*c[1][1]- c[0][1]*c[1][0];
assertF(fabs(tween)>=0.00001,"in min2ByFixMethod1,fabs(tween)<0.00001");
a=(d[0]*c[1][1]-c[0][1]*d[1])/tween;
b=(c[0][0]*d[1]-c[1][0]*d[0])/tween;
printf("the fix fun is %f+%fx.\n",b,a);
fprintf(outputFile,"the fix fun is %f+%fx.\r\n",b,a);
}
Type sumArr(Type* inArr,int len)
{
int i=0;/*iterator num*/
Type ans=0;
assertF(inArr!=NULL,"in sumArr,inArr is NULL\n");
for(i=0;i<len;i++)
ans+=inArr[i];
return ans;
}
void arrEachBy(Type* firstArr,Type* secondArr,Type** ansArr,int len)
{
int i=0;/*iterator num*/
/*assertion*/
assertF(firstArr!=NULL,"in arrEachBy,firstArr is NULL\n");
assertF(secondArr!=NULL,"in arrEachBy,secondArr is NULL\n");
assertF(ansArr!=NULL,"in arrEachBy,ansArr is NULL\n");
for(i=0;i<len;i++)
(*ansArr)[i]=firstArr[i]*secondArr[i];
}
/*Min2ByMethod Algorithm test program*/
#include "Global.h"
#include "Ulti.h"
#include "MyAssert.h"
#include "Fix.h"
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
char *inFileName="inputData.txt";
/*
input data specification
len,
x0,x1,...,xn;
y0,y1,...,yn;
*/
char *outFileName="outputData.txt";
#define DEBUG 1
void main(int argc,char* argv[])
{
FILE *inputFile;/*input file*/
FILE *outputFile;/*output file*/
double startTime,endTime,tweenTime;/*time callopsed info*/
/*The read in data*/
int len;
Type* xList,* yList;
int i;/*iterator index*/
/*input file open*/
if(argc>1)strcpy(inFileName,argv[1]);
assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");
printf("input file open success\n");
/*outpout file open*/
if(argc>2)strcpy(outFileName,argv[2]);
assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");
printf("output file open success\n");
fscanf(inputFile,"%d,",&len);
/*Memory apply*/
xList=(Type*)malloc(sizeof(Type)*len);
yList=(Type*)malloc(sizeof(Type)*len);
/*Read info data*/
for(i=0;i<len-1;i++)
fscanf(inputFile,"%f,",&xList[i]);
fscanf(inputFile,"%f;",&xList[i]);
for(i=0;i<len-1;i++)
fscanf(inputFile,"%f,",&yList[i]);
fscanf(inputFile,"%f;",&yList[i]);
printf("read in data info,len:%d,\n",len);
showArrListFloat(xList,0,len);
showArrListFloat(yList,0,len);
#if DEBUG
printf("\n*******start of test program******\n");
printf("now is runnig,please wait...\n");
startTime=(double)clock()/(double)CLOCKS_PER_SEC;
/******************Core program code*************/
min2ByFixMethod1(xList,yList,len,outputFile);
/******************End of Core program**********/
endTime=(double)clock()/(double)CLOCKS_PER_SEC;
tweenTime=endTime-startTime;/*Get the time collapsed*/
/*Time collapsed output*/
printf("the collapsed time in this algorithm implement is:%f\n",tweenTime);
fprintf(outputFile,"the collapsed time in this algorithm implement is:%f\r\n",tweenTime);
printf("\n*******end of test program******\n");
#endif
free(xList);
free(yList);
printf("program end successfully,\n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer.\n");
getchar();/*Screen Delay Control*/
return;
}
//测试数据1
7,
0.6,1.3,1.64,1.8,2.1,2.3,2.44;
7.05,12.2,14.4,15.2,17.4,19.6,20.2;
//输出:
the fix fun is 2.708288+7.150409x.
//测试数据2
2,
0,1;
3,5;
the fix fun is 3.000000+2.000000x.
