C++最小二乘法拟合-（线性拟合和多项式拟合）
来源：互联网 发布：淘宝卖家纸条 编辑：程序博客网 时间：2024/04/29 00:53
在进行曲线拟合时用的最多的是最小二乘法，其中以一元函数（线性）和多元函数（多项式）居多，下面这个类专门用于进行多项式拟合，可以根据用户输入的阶次进行多项式拟合，算法来自于网上，和GSL的拟合算法对比过，没有问题。此类在拟合完后还能计算拟合之后的误差：SSE(剩余平方和)，SSR(回归平方和)，RMSE(均方根误差)，R-square(确定系数)。
1.fit类的实现
先看看fit类的代码：（只有一个头文件方便使用）
[cpp] view plaincopy
- #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.求解线性拟合：
[cpp] view plaincopy
- ///
- /// \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.多项式拟合：
[cpp] view plaincopy
- ///
- /// \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新建对话框文件，
对话框资源文件如图所示：
加入下面的这些变量：
[cpp] view plaincopy
- std::vector<double> m_x,m_y,m_yploy;
- const size_t m_size;
- CChartLineSerie *m_pLineSerie1;
- CChartLineSerie *m_pLineSerie2;
由于m_size是常量，因此需要在构造函数进行初始化，如：
[cpp] view plaincopy
- ClineFitDlg::ClineFitDlg(CWnd* pParent /*=NULL*/)
- : CDialogEx(ClineFitDlg::IDD, pParent)
- ,m_size(512)
- ,m_pLineSerie1(NULL)
[cpp] view plaincopy
- 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是随机数生成函数，实现如下：
[cpp] view plaincopy
- 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;
- }
运行程序，如图所示
线性拟合的使用如下：
[cpp] view plaincopy
- 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里面，可以通过相关方法调用，方法在头文件中都有详细说明
运行结果如图所示：
多项式拟合的使用如下：
[cpp] view plaincopy
- 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项式拟合，最后的布尔变量是标定是否需要把拟合的结果点保存起来，保存点会根据观察的x值计算拟合的y值，保存结果点会花费更多的内存，如果拟合后需要绘制，设为true会更方便，如果只需要拟合的方程，可以设置为false
- 拟合完后，拟合的相关结果保存在czy::Fit里面，可以通过相关方法调用，方法在头文件中都有详细说明
代码：
[cpp] view plaincopy
- 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就添加“+”号
结果如下：
源代码下载：
C++最小二乘法拟合-（线性拟合和多项式拟合）
0 0
- C++最小二乘法拟合-（线性拟合和多项式拟合）
- C++最小二乘法拟合-（线性拟合和多项式拟合）
- C++最小二乘法拟合-（线性拟合和多项式拟合）
- C++最小二乘法拟合-（线性拟合和多项式拟合）
- 最小二乘法多项式拟合
- 线性拟合(最小二乘法)
- 最小二乘法解线性拟合
- 最小二乘法的线性拟合
- 最小二乘法线性拟合
- 最小二乘法线性拟合计算器
- 最小二乘法的基本原理和多项式拟合
- 最小二乘法详解（线性拟合与非线性拟合）
- 最小二乘法的多项式拟合代码
- 最小二乘法求多项式拟合曲线
- 最小二乘法拟合多项式曲线原理
- 最小二乘法直线拟合(C++)
- pytho作线性拟合、多项式拟合、对数拟合
- 最小二乘法拟合
- 服务端操作MYSQL数据库
- Host '*' is not allowed to connect to this MariaDB server
- JS 跨域访问的几种方式
- 微风 微蓝 之 微风
- 2016年书单和电影
- C++最小二乘法拟合-（线性拟合和多项式拟合）
- 网络设备常用协议总结二(锐捷命令)
- HDU 4179Difficult Routes spfa存了反向边
- 生成二维码
- 编译boost程序出现如下错误fatal error LNK1104: cannot open file 'libboost_system-vc100-mt-gd-1_54.lib'的解决方法
- httpclient介绍
- 数据库工具类DBUtil
- Java resource loading explained.
- android开发中，有关TABhost中各个子Activity跳转的问题（帮大忙）
