zoukankan      html  css  js  c++  java
  • 多项式拟合的cpp实现

    当我们拥有一组散点图数据时,通常更愿意看到其走势。

    对现有数据进行拟合,并输出拟合优度是常用的方法之一。

    拟合结果正确性的验证,可以使用excel自带的功能。

    下面是c++代码的实现:

    #ifndef __Fit_h__
    #define __Fit_h__
    
    #include <vector>
    
    template<size_t Degree>
    class CFit
    {
    public:
    	CFit(std::vector<double>& xArr,std::vector<double>& yArr):m_xArr(xArr),m_yArr(yArr),m_ssr(0.0),m_sse(0.0),m_rmse(0.0){m_bFitYet = false;}
    
    	~CFit(){}
    protected:	
    	//- 高斯消除 
    	template<size_t realDegree>
    	static void GaussianElimination(double (&matrix)[realDegree+1][realDegree+2]
    	,double (&coeArr)[realDegree+1])
    	{
    		int i,j,k;
    		for (i = 0; i< realDegree; i++ )					//loop to perform the gauss elimination
    		{
    			for (k = i+1; k < (realDegree+1); k++)
    			{
    				double t = matrix[k][i]/matrix[i][i];
    				for (j=0;j<=(realDegree+1);j++)
    					matrix[k][j] -= t*matrix[i][j];			//make the elements below the pivot elements equal to zero or elimnate the variables
    			}
    		}
    		for (i = realDegree; i >= 0; i--)					//back-substitution
    		{													//x is an array whose values correspond to the values of x,y,z..
    			coeArr[i] = matrix[i][realDegree+1];    //make the variable to be calculated equal to the rhs of the last equation
    			for (j=0;j<(realDegree+1);j++)
    				if (j!=i)									//then subtract all the lhs values except the coefficient of the variable whose value                                   is being calculated
    					coeArr[i] -= matrix[i][j]*coeArr[j];
    			coeArr[i] = coeArr[i]/matrix[i][i];            //now finally divide the rhs by the coefficient of the variable to be calculated
    		}
    	}
    
    	/// 
    	/// rief 根据x获取拟合方程的y值
    	/// 
    eturn 返回x对应的y值
    	///
    	template<typename T>
    	double getY(const T x) const
    	{
    		double ans(0);
    		for (size_t i=0;i<(Degree+1);++i)
    		{
    			ans += m_coefficientArr[i]*pow((double)x,(int)i);
    		}
    		return ans;
    	}
    
    
    	/// 
    	/// rief 计算均值
    	/// 
    eturn 均值
    	///
    	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);
    	}
    
    	template<typename T>
    	void calcError(const T* x
    		,const T* y
    		,size_t length
    		,double& r_ssr
    		,double& r_sse
    		,double& r_rmse
    		)
    	{
    		T mean_y = Mean<T>(y,length);
    		T yi(0);
    		
    		for (size_t 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]));//残差平方和
    		}
    		r_rmse = sqrt(r_sse/(double(length)));
    	}
    
    
    	/**
    	* @brief    根据两组数据进行一元多项式拟合  
    	* @author  
    	* @param    [in]	int N,数据个数
    				[in]	const std::vector<double>& xArr,横坐标数据
    				[in]	const std::vector<double>& yArr,纵坐标数据
    	* @param    [out]	double (&coefficientArr)[Degree+1],拟合结果.一元多项式系数,从低到高
    	* @return   none
    	* @note     none
    	*/
    	static void PolynomialFit(int N,const std::vector<double>& xArr,const std::vector<double>& yArr,double (&coefficientArr)[Degree+1])
    	{
    		int i = 0,j = 0,k = 0;
    		//const int realDegree = Degree -1;
    		double X[2*Degree+1] = {0};              //Array that will store the values of sigma(xi),sigma(xi^2),sigma(xi^3)....sigma(xi^2n)
    		for (i=0;i<2*Degree+1;i++)
    		{
    			for (j=0;j<N;j++)
    				X[i] += pow(xArr[j],i);              //consecutive positions of the array will store N,sigma(xi),sigma(xi^2),sigma(xi^3)....sigma(xi^2n)
    		}
    		double Y[Degree+1] = {0};                //Array to store the values of sigma(yi),sigma(xi*yi),sigma(xi^2*yi)...sigma(xi^n*yi)
    		for (i=0;i<Degree+1;i++)
    		{    
    			for (j=0;j<N;j++)
    				Y[i] += pow(xArr[j],i)*yArr[j];		//consecutive positions will store sigma(yi),sigma(xi*yi),sigma(xi^2*yi)...sigma(xi^n*yi)
    		}
    
    		double B[Degree+1][Degree+2] = {0};			//B is the Normal matrix(augmented) that will store the equations
    		for (i=0;i<=Degree;i++)
    		{
    			for (j=0;j<=Degree;j++)
    			{
    				B[i][j] = X[i+j];					//Build the Normal matrix by storing the corresponding coefficients at the right positions except the last column of the matrix
    			}
    			B[i][Degree+1] = Y[i];				//load the values of Y as the last column of B(Normal Matrix but augmented)
    		}
    
    		GaussianElimination<Degree>(B,coefficientArr);
    	}
    public:
    	void PolyFit()
    	{
    		if (m_xArr.size() == m_yArr.size())
    		{
    			PolynomialFit(static_cast<int>(m_xArr.size()),m_xArr,m_yArr,m_coefficientArr);
    			m_bFitYet = true;
    			calcError(&m_xArr[0],&m_yArr[0],static_cast<int>(m_xArr.size()),m_ssr,m_sse,m_rmse);
    		}
    	}
    	//- 一元多项式计算 
    	double UnaryPolynomialCalc(double dX)
    	{
    		double dY = 0.0;
    		for (size_t ulDegree = 0; ulDegree <= Degree; ++ulDegree)
    		{
    			dY += pow(dX,(double)ulDegree) * m_coefficientArr[ulDegree];
    		}
    		return m_bFitYet ? dY : 0.0;
    	}
    
    	/// 
    	/// rief 剩余平方和
    	/// 
    eturn 剩余平方和
    	///
    	double getSSE(){return m_sse;}
    	/// 
    	/// rief 回归平方和
    	/// 
    eturn 回归平方和
    	///
    	double getSSR(){return m_ssr;}
    	/// 
    	/// rief 均方根误差
    	/// 
    eturn 均方根误差
    	///
    	double getRMSE(){return m_rmse;}
    	/// 
    	/// rief 确定系数,系数是0~1之间的数,是数理上判定拟合优度(goodness-of-fit)的一个量
    	/// 
    eturn 确定系数
    	///
    	double getR_square(){return 1-(m_sse/(m_ssr+m_sse));}
    
    	/// 
    	/// rief 根据阶次获取拟合方程的系数,
    	/// 如getFactor(2),就是获取y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n中a2的值
    	/// 
    eturn 拟合方程的系数
    	///
    	double getFactor(size_t i)
    	{
    		return (i <= Degree) ? m_coefficientArr[i] : 0.0;
    	}
    
    
    	
    private:
    	double m_coefficientArr[Degree+1];
    	const std::vector<double>& m_xArr;
    	const std::vector<double>& m_yArr;
    	bool m_bFitYet;//- 一元多项式计算时多项式拟合是否完成 [1/6/2017 wWX210786]
    
    	double m_ssr;                 ///<回归平方和
    	double m_sse;                 ///<(剩余平方和)
    	double m_rmse;                ///<RMSE均方根误差
    };
    
    
    #endif // __Fit_h__
    

      使用起来也很方便:

            double y[] = {7,16,6,18,6,6,10,8};
    		double x[] = {-109.71,-101.81,-103.83,-99.89,-90,-112.17,-93.5,-96.13};
    		std::vector<double> xArr(std::begin(x),std::end(x));
    		std::vector<double> yArr(std::begin(y),std::end(y));
    		typedef CFit<4> LineFit;
    		LineFit objPolyfit(xArr,yArr);
    		objPolyfit.PolyFit();                                
    		std::wstring coeArr[] = {L"",L"x",L"x²",L"xu00b3",L"x "};
    		CString info(_T("y = "));
    		for (int i=1;i>=0;i--)
    			info.AppendFormat(_T("+( %f%s )"),objPolyfit.m_coefficientArr[i],coeArr[i].c_str());
    
    		std::wcout << info.GetString() << "
    ";
    
    		//std::wcout << "斜率 = " << objPolyfit.getFactor(1) << "
    ";
    		//std::wcout << "截距 = " << objPolyfit.getFactor(0) << "
    ";
    		std::wcout << "goodness-of-fit = "<< objPolyfit.getR_square() << "
    ";
    

      

  • 相关阅读:
    第十三周编程总结
    第十二周总结
    第十一周课程总结
    第十周编程总结
    第七次实验报告及编程总结
    第六次实验报告及学习总结
    课程总结
    第十四周课程总结&实验报告
    第十三周课程总结
    第十二周总结
  • 原文地址:https://www.cnblogs.com/envoy/p/8492484.html
Copyright © 2011-2022 走看看