多项式拟合的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
}
}
///
/// \brief 根据x获取拟合方程的y值
/// \return 返回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;
}
///
/// \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);
}
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;
}
///
/// \brief 剩余平方和
/// \return 剩余平方和
///
double getSSE(){return m_sse;}
///
/// \brief 回归平方和
/// \return 回归平方和
///
double getSSR(){return m_ssr;}
///
/// \brief 均方根误差
/// \return 均方根误差
///
double getRMSE(){return m_rmse;}
///
/// \brief 确定系数，系数是0~1之间的数，是数理上判定拟合优度(goodness-of-fit)的一个量
/// \return 确定系数
///
double getR_square(){return 1-(m_sse/(m_ssr+m_sse));}
///
/// \brief 根据阶次获取拟合方程的系数，
/// 如getFactor(2),就是获取y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n中a2的值
/// \return 拟合方程的系数
///
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"x\u00b3",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() << "\n";
//std::wcout << "斜率 = " << objPolyfit.getFactor(1) << "\n";
//std::wcout << "截距 = " << objPolyfit.getFactor(0) << "\n";
std::wcout << "goodness-of-fit = "<< objPolyfit.getR_square() << "\n";
