C++最小二乘法拟合-(线性拟合和多项式拟合)(一)

2014-11-24 11:12:24 · 作者: · 浏览: 3


在进行曲线拟合时用的最多的是最小二乘法,其中以一元函数(线性)和多元函数(多项式)居多,下面这个类专门用于进行多项式拟合,可以根据用户输入的阶次进行多项式拟合,算法来自于网上,和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 拟合后的数据是否保存,默认否
///
template
                                       
                                        
bool 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 观察