c++ 多项式拟合算法

算法 专栏收录该内容
2 篇文章 0 订阅
#ifndef CZY_MATH_FIT
#define CZY_MATH_FIT
#include <vector>
/*
多项式拟合
*/
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, int 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, typename TD>
		void polyfit(const std::vector<typename T>& x
			, const std::vector<typename TD>& y
			, int poly_n
			, bool isSaveFitYs = true)
		{
			polyfit(&x[0], &y[0], getSeriesLength(x, y), poly_n, isSaveFitYs);
		}
		template<typename T, typename TD>
		void polyfit(const T* x, const TD* y, int 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 (int i = 0; i < (int)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, typename TD>
		int getSeriesLength(const std::vector<typename T>& x
			, const std::vector<typename TD>& 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, int length)
		{
			T total(0);
			for (int i = 0; i < length; ++i)
			{
				total += v[i];
			}
			return (total / length);
		}
		/// 
		/// \brief 获取拟合方程系数的个数
		/// \return 拟合方程系数的个数
		///
		int getFactorSize(){ return factor.size(); }
		/// 
		/// \brief 根据阶次获取拟合方程的系数,
		/// 如getFactor(2),就是获取y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n中a2的值
		/// \return 拟合方程的系数
		///
		double getFactor(int i){ return factor.at(i); }
	private:
		template<typename T, typename TD>
		void calcError(const T* x
			, const TD* y
			, int length
			, double& r_ssr
			, double& r_sse
			, double& r_rmse
			, bool isSaveFitYs = true
			)
		{
			TD mean_y = Mean<TD>(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
用法:
void CDTData::Polyfit_Fun(WSDATA &wData, int nOrder, int nDataCnt)
{
	vector<double> vecX;
	vector<double> vecY;
	int nBeginPoints = 0; // 开始拟合点
	//保存所有点用于拟合
	for (int i = 0; i < nDataCnt - nBeginPoints; ++i)
	{
		vecX.push_back(i);
		vecY.push_back(wData[i + nBeginPoints]);
	}

	int poly_n = nOrder;	//拟合阶数

	//拟合
	czy::Fit fit;
	fit.polyfit(vecX, vecY, poly_n);

	//获取所有拟合后的值
	vector<double> vecYPloy;
	fit.getFitedYs(vecYPloy);

	//减去拟合值
	for (int i = nBeginPoints; i < nDataCnt; ++i)
	{
		wData[i] = wData[i] - vecYPloy[i - nBeginPoints];
	}
}

忘记从哪里弄的了

  • 0
    点赞
  • 0
    评论
  • 11
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

参与评论 您还未登录,请先 登录 后发表或查看评论
#ifndef FUNCTION_H_ #define FUNCTION_H_ #include #include #include "polyfit.h" #include using namespace std; dxs::dxs() { ifstream fin("多项式拟合.txt"); fin>>n; x=new float[n]; y=new float[n]; for(int i=0;i>x[i]; } for(i=0;i>y[i]; } cout<>nn; m=nn+1; u=new float*[m]; for(i=0;i<m;i++) { u[i]=new float[m+1]; }//创建m行,m+1列数组 } void dxs::dfine() { for(int i=0;i<m;i++) { for(int j=0;j<m+1;j++) { u[i][j]=0; } } for(i=0;i<m;i++) { for(int j=0;j<m;j++) { for(int k=0;k<n;k++) { u[i][j]=u[i][j]+pow(x[k],j+i); } } } for(i=0;i<m;i++) { for(int k=0;k<n;k++) { u[i][m]=u[i][m]+pow(x[k],i)*y[k]; } } } void dxs::show() { for(int i=0;i<m;i++) { for(int j=0;j<m+1;j++) { cout<<u[i][j]<<" ";//<<endl; } cout<<endl; } ////显示具有m行m+1列u数组的各元素值 } void dxs::select_main(int k,float **p,int m) { double d; d=*(*(p+k)+k); //cout<<d; int l=k; int i=k+1; for(;i fabs(d)) { d=*(*(p+i)+k); l=i; } else continue; } if(d==0) cout<<"错误"; else { if(k!=l) { for(int j=k;j<m+1;j++) { double t; t=*(*(p+l)+j); *(*(p+l)+j)=*(*(p+k)+j); *(*(p+k)+j)=t; } } } } void dxs::gaosi() { for(int k=0;k<m;k++) { select_main(k,u,m);//调用列主元函数 for(int i=1+k;i<m;i++) { // *(*(p+i)+k)=(float) *(*(p+i)+k) / *(*(p+k)+k); u[i][k]=(float) u[i][k] / u[k][k]; } for(i=k+1;i<m;i++) { for(int j=k+1;j=0;i--) { float a=0; for(int j=i+1;j<m;j++) { //a=a + (*(*(p+i)+j) * *(*(p+j)+m)); a=a+u[i][j] * u[j][m]; } //*(*(p+i)+n-1)= (*(*(p+i)+n-1) - a) / *(*(p+i)+i); u[i][m]= (u[i][m] -a) / u[i][i]; } cout<<"方程组的解为:"<<endl; for(i=0;i<m;i++) { cout<<"a"<<i+1<<"="; cout<<u[i][m]<<endl; // l[i]=*(*(p+i)+n-1); } cout<<"y="<<u[0][m]; for(i=1;i<m;i++) { cout<<showpos<<u[i][m]<<"x"; if(i!=1)cout<<"^"<<noshowpos<<i; } cout<<endl; } dxs::~dxs() { delete[]x,y; delete []*u; } #endif
©️2022 CSDN 皮肤主题:鲸 设计师:meimeiellie 返回首页

打赏作者

IT_Kyle

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值