在进行曲线拟合时用的最多的是最小二乘法,其中以一元函数(线性)和多元函数(多项式)居多,下面这个类专门用于进行多项式拟合,可以根据用户输入的阶次进行多项式拟合,算法来自于网上,和GSL的拟合算法对比过,没有问题。此类在拟合完后还能计算拟合之后的误差:SSE(剩余平方和),SSR(回归平方和),RMSE(均方根误差),R-square(确定系数)。
1.fit类的实现
先看看fit类的代码:(只有一个头文件方便使用)
#ifndef CZY_MATH_FIT #define CZY_MATH_FIT #include/* 尘中远,于2014.03.20 主页:http://blog.csdn.net/czyt1988/article/details/21743595 参考:http://blog.csdn.net/maozefa/article/details/1725535 */ namespace czy{ /// /// \brief 曲线拟合类 /// class Fit{ std::vector factor; ///<拟合后的方程系数 double ssr; ///<回归平方和 double sse; ///<(剩余平方和) double rmse; /// fitedYs;///<存放拟合后的y值,在拟合时可设置为不保存节省内存 public: Fit():ssr(0),sse(0),rmse(0){factor.resize(2,0);} ~Fit(){} /// /// \brief 直线拟合-一元回归,拟合的结果可以使用getFactor获取,或者使用getSlope获取斜率,getIntercept获取截距 /// \param x 观察值的x /// \param y 观察值的y /// \param isSaveFitYs 拟合后的数据是否保存,默认否 /// template bool linearFit(const std::vector & x, const std::vector & y,bool isSaveFitYs=false) { return linearFit(&x[0],&y[0],getSeriesLength(x,y),isSaveFitYs); } template bool linearFit(const T* x, const T* y,size_t length,bool isSaveFitYs=false) { factor.resize(2,0); typename T t1=0, t2=0, t3=0, t4=0; for(int i=0; i ssr,this->sse,this->rmse,isSaveFitYs); return true; } /// /// \brief 多项式拟合,拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n /// \param x 观察值的x /// \param y 观察值的y /// \param poly_n 期望拟合的阶数,若poly_n=2,则y=a0+a1*x+a2*x^2 /// \param isSaveFitYs 拟合后的数据是否保存,默认是 /// template void polyfit(const std::vector & x ,const std::vector & y ,int poly_n ,bool isSaveFitYs=true) { polyfit(&x[0],&y[0],getSeriesLength(x,y),poly_n,isSaveFitYs); } template void polyfit(const T* x,const T* y,size_t length,int poly_n,bool isSaveFitYs=true) { factor.resize(poly_n+1,0); int i,j; //double *tempx,*tempy,*sumxx,*sumxy,*ata; std::vector tempx(length,1.0); std::vector tempy(y,y+length); std::vector sumxx(poly_n*2+1); std::vector ata((poly_n+1)*(poly_n+1)); std::vector sumxy(poly_n+1); for (i=0;i<2*poly_n+1;i++){ for (sumxx[i]=0,j=0;j ssr,this->sse,this->rmse,isSaveFitYs); } /// /// \brief 获取系数 /// \param 存放系数的数组 /// void getFactor(std::vector & factor){factor = this->factor;} /// /// \brief 获取拟合方程对应的y值,前提是拟合时设置isSaveFitYs为true /// void getFitedYs(std::vector & fitedYs){fitedYs = this->fitedYs;} /// /// \brief 根据x获取拟合方程的y值 /// \return 返回x对应的y值 /// template double getY(const T x) const { double ans(0); for (size_t i=0;i size_t getSeriesLength(const std::vector & x ,const std::vector & y) { return (x.size() > y.size() y.size() : x.size()); } /// /// \brief 计算均值 /// \return 均值 /// template static T Mean(const std::vector & v) { return Mean(&v[0],v.size()); } template static T Mean(const T* v,size_t length) { T total(0); for (size_t i=0;i void calcError(const T* x ,const T* y ,size_t length ,double& r_ssr ,double& r_sse ,double& r_rmse ,bool isSaveFitYs=true ) { T mean_y = Mean (y,length); T yi(0); fitedYs.reserve(length); for (int i=0; i void gauss_solve(int n ,std::vector & A ,std::vector & x ,std::vector & b) { gauss_solve(n,&A[0],&x[0],&b[0]); } template void gauss_solve(int n ,T* A ,T* x ,T* b) { int i,j,k,r; double max; for (k=0;k =0;x[i]/=A[i*n+i],i--) for (j=i+1,x[i]=b[i];j
为了防止重命名,把其放置于czy的命名空间中,此类主要两个函数:
1.求解线性拟合:
/// /// \brief 直线拟合-一元回归,拟合的结果可以使用getFactor获取,或者使用getSlope获取斜率,getIntercept获取截距 /// \param x 观察值的x /// \param y 观察值的y /// \param length x,y数组的长度 /// \param isSaveFitYs 拟合后的数据是否保存,默认否 /// templatebool linearFit(const std::vector & x, const std::vector & y,bool isSaveFitYs=false); template bool linearFit(const T* x, const T* y,size_t length,bool isSaveFitYs=false);
2.多项式拟合:/// /// \brief 多项式拟合,拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n /// \param x 观察