在进行曲线拟合时用的最多的是最小二乘法,其中以一元函数(线性)和多元函数(多项式)居多,下面这个类专门用于进行多项式拟合,可以根据用户输入的阶次进行多项式拟合,算法来自于网上,和GSL的拟合算法对比过,没有问题。此类在拟合完后还能计算拟合之后的误差:SSE(剩余平方和),SSR(回归平方和),RMSE(均方根误差),R-square(确定系数)。
1.fit类的实现
先看看fit类的代码:(只有一个头文件方便使用)
#ifndef CZY_MATH_FIT #define CZY_MATH_FIT #include <vector> /* 尘中远,于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<double> factor; ///<拟合后的方程系数 double ssr; ///<回归平方和 double sse; ///<(剩余平方和) double rmse; ///<RMSE均方根误差 std::vector<double> 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<typename T> bool linearFit(const std::vector<typename T>& x, const std::vector<typename T>& y,bool isSaveFitYs=false) { return linearFit(&x[0],&y[0],getSeriesLength(x,y),isSaveFitYs); } template<typename T> 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<length; ++i) { t1 += x[i]*x[i]; t2 += x[i]; t3 += x[i]*y[i]; t4 += y[i]; } factor[1] = (t3*length - t2*t4) / (t1*length - t2*t2); factor[0] = (t1*t4 - t2*t3) / (t1*length - t2*t2); ////////////////////////////////////////////////////////////////////////// //计算误差 calcError(x,y,length,this->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<typename T> void polyfit(const std::vector<typename T>& x ,const std::vector<typename T>& y ,int poly_n ,bool isSaveFitYs=true) { polyfit(&x[0],&y[0],getSeriesLength(x,y),poly_n,isSaveFitYs); } template<typename T> 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<double> tempx(length,1.0); std::vector<double> tempy(y,y+length); std::vector<double> sumxx(poly_n*2+1); std::vector<double> ata((poly_n+1)*(poly_n+1)); std::vector<double> sumxy(poly_n+1); for (i=0;i<2*poly_n+1;i++){ for (sumxx[i]=0,j=0;j<length;j++) { sumxx[i]+=tempx[j]; tempx[j]*=x[j]; } } for (i=0;i<poly_n+1;i++){ for (sumxy[i]=0,j=0;j<length;j++) { sumxy[i]+=tempy[j]; tempy[j]*=x[j]; } } for (i=0;i<poly_n+1;i++) for (j=0;j<poly_n+1;j++) ata[i*(poly_n+1)+j]=sumxx[i+j]; gauss_solve(poly_n+1,ata,factor,sumxy); //计算拟合后的数据并计算误差 fitedYs.reserve(length); calcError(&x[0],&y[0],length,this->ssr,this->sse,this->rmse,isSaveFitYs); } /// /// \brief 获取系数 /// \param 存放系数的数组 /// void getFactor(std::vector<double>& factor){factor = this->factor;} /// /// \brief 获取拟合方程对应的y值,前提是拟合时设置isSaveFitYs为true /// void getFitedYs(std::vector<double>& fitedYs){fitedYs = this->fitedYs;} /// /// \brief 根据x获取拟合方程的y值 /// \return 返回x对应的y值 /// template<typename T> double getY(const T x) const { double ans(0); for (size_t i=0;i<factor.size();++i) { ans += factor[i]*pow((double)x,(int)i); } return ans; } /// /// \brief 获取斜率 /// \return 斜率值 /// double getSlope(){return factor[1];} /// /// \brief 获取截距 /// \return 截距值 /// double getIntercept(){return factor[0];} /// /// \brief 剩余平方和 /// \return 剩余平方和 /// double getSSE(){return sse;} /// /// \brief 回归平方和 /// \return 回归平方和 /// double getSSR(){return ssr;} /// /// \brief 均方根误差 /// \return 均方根误差 /// double getRMSE(){return rmse;} /// /// \brief 确定系数,系数是0~1之间的数,是数理上判定拟合优度的一个量 /// \return 确定系数 /// double getR_square(){return 1-(sse/(ssr+sse));} /// /// \brief 获取两个vector的安全size /// \return 最小的一个长度 /// template<typename T> size_t getSeriesLength(const std::vector<typename T>& x ,const std::vector<typename T>& y) { return (x.size() > y.size() ? y.size() : x.size()); } /// /// \brief 计算均值 /// \return 均值 /// template <typename T> static T Mean(const std::vector<T>& v) { return Mean(&v[0],v.size()); } template <typename T> static T Mean(const T* v,size_t length) { T total(0); for (size_t i=0;i<length;++i) { total += v[i]; } return (total / length); } /// /// \brief 获取拟合方程系数的个数 /// \return 拟合方程系数的个数 /// size_t getFactorSize(){return factor.size();} /// /// \brief 根据阶次获取拟合方程的系数, /// 如getFactor(2),就是获取y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n中a2的值 /// \return 拟合方程的系数 /// double getFactor(size_t i){return factor.at(i);} private: template<typename T> 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<T>(y,length); T yi(0); fitedYs.reserve(length); for (int i=0; i<length; ++i) { yi = getY(x[i]); r_ssr += ((yi-mean_y)*(yi-mean_y));//计算回归平方和 r_sse += ((yi-y[i])*(yi-y[i]));//残差平方和 if (isSaveFitYs) { fitedYs.push_back(double(yi)); } } r_rmse = sqrt(r_sse/(double(length))); } template<typename T> void gauss_solve(int n ,std::vector<typename T>& A ,std::vector<typename T>& x ,std::vector<typename T>& b) { gauss_solve(n,&A[0],&x[0],&b[0]); } template<typename T> void gauss_solve(int n ,T* A ,T* x ,T* b) { int i,j,k,r; double max; for (k=0;k<n-1;k++) { max=fabs(A[k*n+k]); /*find maxmum*/ r=k; for (i=k+1;i<n-1;i++){ if (max<fabs(A[i*n+i])) { max=fabs(A[i*n+i]); r=i; } } if (r!=k){ for (i=0;i<n;i++) /*change array:A[k]&A[r] */ { max=A[k*n+i]; A[k*n+i]=A[r*n+i]; A[r*n+i]=max; } } max=b[k]; /*change array:b[k]&b[r] */ b[k]=b[r]; b[r]=max; for (i=k+1;i<n;i++) { for (j=k+1;j<n;j++) A[i*n+j]-=A[i*n+k]*A[k*n+j]/A[k*n+k]; b[i]-=A[i*n+k]*b[k]/A[k*n+k]; } } for (i=n-1;i>=0;x[i]/=A[i*n+i],i--) for (j=i+1,x[i]=b[i];j<n;j++) x[i]-=A[i*n+j]*x[j]; } }; } #endif
为了防止重命名,把其放置于czy的命名空间中,此类主要两个函数:
1.求解线性拟合:
/// /// \brief 直线拟合-一元回归,拟合的结果可以使用getFactor获取,或者使用getSlope获取斜率,getIntercept获取截距 /// \param x 观察值的x /// \param y 观察值的y /// \param length x,y数组的长度 /// \param isSaveFitYs 拟合后的数据是否保存,默认否 /// template<typename T> bool linearFit(const std::vector<typename T>& x, const std::vector<typename T>& y,bool isSaveFitYs=false); template<typename T> 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 观察值的x /// \param y 观察值的y /// \param length x,y数组的长度 /// \param poly_n 期望拟合的阶数,若poly_n=2,则y=a0+a1*x+a2*x^2 /// \param isSaveFitYs 拟合后的数据是否保存,默认是 /// template<typename T> void polyfit(const std::vector<typename T>& x,const std::vector<typename T>& y,int poly_n,bool isSaveFitYs=true); template<typename T> void polyfit(const T* x,const T* y,size_t length,int poly_n,bool isSaveFitYs=true);
这两个函数都用模板函数形式写,主要是为了能使用于float和double两种数据类型
2.fit类的MFC示范程序
下面看看如何使用这个类,以MFC示范,使用了开源的绘图控件Hight-Speed Charting,使用方法见http://blog.csdn.net/czyt1988/article/details/8740500新建对话框文件,
对话框资源文件如图所示:
加入下面的这些变量:
std::vector<double> m_x,m_y,m_yploy; const size_t m_size; CChartLineSerie *m_pLineSerie1; CChartLineSerie *m_pLineSerie2;
由于m_size是常量,因此需要在构造函数进行初始化,如:
ClineFitDlg::ClineFitDlg(CWnd* pParent /*=NULL*/) : CDialogEx(ClineFitDlg::IDD, pParent) ,m_size(512) ,m_pLineSerie1(NULL)
CChartAxis *pAxis = NULL; pAxis = m_chartCtrl.CreateStandardAxis(CChartCtrl::BottomAxis); pAxis->SetAutomatic(true); pAxis = m_chartCtrl.CreateStandardAxis(CChartCtrl::LeftAxis); pAxis->SetAutomatic(true); m_x.resize(m_size); m_y.resize(m_size); m_yploy.resize(m_size); for(size_t i =0;i<m_size;++i) { m_x[i] = i; m_y[i] = i+randf(-25,28); m_yploy[i] = 0.005*pow(double(i),2)+0.0012*i+4+randf(-25,25); } m_chartCtrl.RemoveAllSeries();//先清空 m_pLineSerie1 = m_chartCtrl.CreateLineSerie(); m_pLineSerie1->SetSeriesOrdering(poNoOrdering);//设置为无序 m_pLineSerie1->AddPoints(&m_x[0], &m_y[0], m_size); m_pLineSerie1->SetName(_T("线性数据")); m_pLineSerie2 = m_chartCtrl.CreateLineSerie(); m_pLineSerie2->SetSeriesOrdering(poNoOrdering);//设置为无序 m_pLineSerie2->AddPoints(&m_x[0], &m_yploy[0], m_size); m_pLineSerie2->SetName(_T("多项式数据"));
rangf是随机数生成函数,实现如下:
double ClineFitDlg::randf(double min,double max) { int minInteger = (int)(min*10000); int maxInteger = (int)(max*10000); int randInteger = rand()*rand(); int diffInteger = maxInteger - minInteger; int resultInteger = randInteger % diffInteger + minInteger; return resultInteger/10000.0; }
运行程序,如图所示
线性拟合的使用如下:
void ClineFitDlg::OnBnClickedButton1() { CString str,strTemp; czy::Fit fit; fit.linearFit(m_x,m_y); str.Format(_T("方程:y=%gx+%g\r\n误差:ssr:%g,sse=%g,rmse:%g,确定系数:%g"),fit.getSlope(),fit.getIntercept() ,fit.getSSR(),fit.getSSE(),fit.getRMSE(),fit.getR_square()); GetDlgItemText(IDC_EDIT,strTemp); SetDlgItemText(IDC_EDIT,strTemp+_T("\r\n------------------------\r\n")+str); //在图上绘制拟合的曲线 CChartLineSerie* pfitLineSerie1 = m_chartCtrl.CreateLineSerie(); std::vector<double> x(2,0),y(2,0); x[0] = 0;x[1] = m_size-1; y[0] = fit.getY(x[0]);y[1] = fit.getY(x[1]); pfitLineSerie1->SetSeriesOrdering(poNoOrdering);//设置为无序 pfitLineSerie1->AddPoints(&x[0], &y[0], 2); pfitLineSerie1->SetName(_T("拟合方程"));//SetName的作用将在后面讲到 pfitLineSerie1->SetWidth(2); }
需要如下步骤:
- 声明Fit类,用于头文件在czy命名空间中,因此需要显示声明命名空间名称czy::Fit fit;
- 把观察数据输入进行拟合,由于是线性拟合,可以使用LinearFit函数,此函数把观察量的x值和y值传入即可进行拟合
- 拟合完后,拟合的相关结果保存在czy::Fit里面,可以通过相关方法调用,方法在头文件中都有详细说明
运行结果如图所示:
多项式拟合的使用如下:
void ClineFitDlg::OnBnClickedButton2() { CString str; GetDlgItemText(IDC_EDIT1,str); if (str.IsEmpty()) { MessageBox(_T("请输入阶次"),_T("警告")); return; } int n = _ttoi(str); if (n<0) { MessageBox(_T("请输入大于1的阶数"),_T("警告")); return; } czy::Fit fit; fit.polyfit(m_x,m_yploy,n,true); CString strFun(_T("y=")),strTemp(_T("")); for (int i=0;i<fit.getFactorSize();++i) { if (0 == i) { strTemp.Format(_T("%g"),fit.getFactor(i)); } else { double fac = fit.getFactor(i); if (fac<0) { strTemp.Format(_T("%gx^%d"),fac,i); } else { strTemp.Format(_T("+%gx^%d"),fac,i); } } strFun += strTemp; } str.Format(_T("方程:%s\r\n误差:ssr:%g,sse=%g,rmse:%g,确定系数:%g"),strFun ,fit.getSSR(),fit.getSSE(),fit.getRMSE(),fit.getR_square()); GetDlgItemText(IDC_EDIT,strTemp); SetDlgItemText(IDC_EDIT,strTemp+_T("\r\n------------------------\r\n")+str); //绘制拟合后的多项式 std::vector<double> yploy; fit.getFitedYs(yploy); CChartLineSerie* pfitLineSerie1 = m_chartCtrl.CreateLineSerie(); pfitLineSerie1->SetSeriesOrdering(poNoOrdering);//设置为无序 pfitLineSerie1->AddPoints(&m_x[0], &yploy[0], yploy.size()); pfitLineSerie1->SetName(_T("多项式拟合方程"));//SetName的作用将在后面讲到 pfitLineSerie1->SetWidth(2); }
步骤如下:
- 和线性拟合一样,什么Fit变量
- 输入观察值,同时输入需要拟合的阶次,这里输入2阶,就是2项式拟合
- 拟合完后,拟合的相关结果保存在czy::Fit里面,可以通过相关方法调用,方法在头文件中都有详细说明
代码:
for (int i=0;i<fit.getFactorSize();++i) { if (0 == i) { strTemp.Format(_T("%g"),fit.getFactor(i)); } else { double fac = fit.getFactor(i); if (fac<0) { strTemp.Format(_T("%gx^%d"),fac,i); } else { strTemp.Format(_T("+%gx^%d"),fac,i); } } strFun += strTemp; }
是用于生成方程的,由于系数小于时,打印时会把负号“-”显示,而正数时却不会显示正号,因此需要进行判断,如果小于0就不用添加“+”号,如果大于0就添加“+”号
结果如下:
源代码下载: