最小二乘法多项式拟合的C++实现、MATLAB实现及验证
关于最小二乘法的原理之类的博客有很多，我在这里就不讲了，下面直接贴出是怎么实现的及验证：
【C++实现】：
1.新建一个基于对话框的工程文件，并建立如下界面，灰色的控件是Custom Control，其属性设置如下，左边的文本编辑框ID为：IDC_EDIT，右边一个文本编辑框ID为：IDC_EDIT1，上面的一个按钮是button1，下面的是button2
2.到CodeProject网站去下载High-speed Charting Control相关的文件，下载这个文件要注册等信息，如果嫌麻烦的话我下面会上传资源，可以免费下载的，下载这个文件后，解压，把文件夹名称改为ChartCtrl，方便后面使用，并将文件夹存放到自己新建工程的与源文件同文件夹下，如下
3.把文件夹下的头文件和源文件都添加到此工程中，项目，添加现有项，全部选中，确定。
4.把该文件夹包含到此工程的所在路径，因为此文件夹是我们自己添加的，如果我们在cpp文件中包含这些头文件的话，会提示错误：无法找到源文件。这就是没有把改文件夹添加到路径的原因。把刚才的文件夹添加到路径就可以了
5.在JiaoChengDlg.h中包含头文件，并在对话框类添加变量，叫m_ChartCtrl，在JiaoChengDlg.cpp中的CJiaoChengDlg::DoDataExchange函数里添加关联
CChartCtrl m_ChartCtrl; DDX_Control(pDX, IDC_ChartCtrl, m_ChartCtrl);6.编译运行可以得到如下界面
7.添加头文件Method.h，并在JiaoCheng.cpp中包含该头文件，代码如下：
#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,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
8.在JiaoChengDlg.h中添加一些全局变量，其中m_x，m_y存放的是你的离散点，m_size存放的是离散点的个数
std::vector<double> m_x,m_y,m_yploy;
size_t m_size;
CChartLineSerie *m_pLineSerie1;9.因为下面在JiaoChengDlg.cpp中用到了ChartLineSerie的指针，因此要在该源文件中包含头文件ChartLineSerie.h
#include "ChartLineSerie.h"10.在“离散点”按钮控件下添加如下代码，上面的arry1,arry2是我自己的离散点，你也可以放你自己的
void CJiaoChengDlg::OnBnClickedButton1()
{
// TODO: 在此添加控件通知处理程序代码
m_size=89;
// TODO: 在此添加控件通知处理程序代码
double arry1[89]={-10.4799576938179,-10.4427131117965,-10.3935955278671,-10.3426374470642,-10.2899673769581,-10.2358344513843,-10.1798183736441,-10.1224321068402,-10.0639831025028,-10.0038342606086,-9.94242846041357,-9.87925462149347,-9.81544120185450,-9.75025020464992,-9.68356278430574,-9.61685157675578,-9.54765066199345,-9.47795346582250,-9.40771353803564,-9.33531859017617,-9.26329619851466,-9.18988671334297,-9.11494202978480,-9.04046294563369,-8.96417718978457,-8.88696934573055,-8.80932993147811,-8.73133935992608,-8.65270157047910,-8.57190471187316,-8.49214268506696,-8.41160004102550,-8.33058324409408,-8.24944156448056,-8.16806366584155,-8.08715690190087,-8.00394960062074,-7.92220789617122,-7.84041190849538,-7.75754714348643,-7.67542521642965,-7.59267856421031,-7.51114791596066,-7.42876189264231,-7.34605306709827,-7.26368743725954,-7.18154018461744,-7.09998873001728,-7.01859514528303,-6.93781598917838,-6.85697571266515,-6.77520076634015,-6.69528804142964,-6.61602324004823,-6.53650092857705,-6.45906160005326,-6.37944789272735,-6.30212216338825,-6.22656177121927,-6.14921395947804,-6.07307680182649,-5.99931870060032,-5.92638025774982,-5.85311785355721,-5.78181144672067,-5.71015461474391,-5.64047367200112,-5.57310903359448,-5.50536121942027,-5.44003379673437,-5.37566110675296,-5.31203042530265,-5.25147688370231,-5.19275851059893,-5.13486204325950,-5.07871375008584,-5.02513462588582,-4.97420709775678,-4.92515882779729,-4.87937767370728,-4.83658274580924,-4.79691628400688,-4.76065314560399,-4.72800212324514,-4.70101203366723,-4.67786846033676,-4.66090507222426,-4.6518746405410,-4.64857668595693};
double arry2[89]={94.5486774453437,94.6091950327735,94.6863803006767,94.7635089756550,94.8404442356320,94.9167759923934,94.9930338215959,95.0687009656611,95.1431050512410,95.2171828575088,95.2903318550923,95.3632096312433,95.4345980252129,95.5050702910846,95.5750333666329,95.6427123978068,95.7106758431543,95.7771054939683,95.8418321849418,95.9065005376063,95.968815689355,96.0301934305458,96.0909289591861,96.1492869455995,96.2070222559737,96.2636020723178,96.3185372735893,96.3718634748319,96.4236645702007,96.4750637610061,96.5239981763049,96.5715019314410,96.6174967295572,96.6616735213237,96.7041828081614,96.7445702765672,96.7842651534391,96.8215358646009,96.8569247901061,96.8910333529152,96.9229949999752,96.9534418537805,96.9816427438242,97.0083214227293,97.0333430048479,97.0564980312684,97.0778307977434,97.0972405538661,97.1148432675382,97.1305528674368,97.1444814112080,97.1568278581921,97.1671840698493,97.1756464845404,97.1824017468487,97.1872169208711,97.1903323575983,97.1917119832175,97.1911851185799,97.1888997725816,97.1848552144178,97.1792118328244,97.1717406206538,97.1624583430699,97.1515834896052,97.1388060921078,97.1245723818418,97.1089216562317,97.0912669369991,97.0723896175601,97.0517896227558,97.0295347379153,97.0063753019352,96.9819354313352,96.9557746422944,96.9283997582167,96.9001946360233,96.8713352376738,96.8413962556383,96.8114152212890,96.7812273352373,96.7512414276797,96.7217213180292,96.6932375286215,96.6678835297399,96.6445093566755,96.6260818062986,96.6155504382452,96.6113860562677};
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);
for(size_t i =0;i<m_size;++i)
{
m_x[i] = arry1[i];
m_y[i] = arry2[i];
}
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("数据"));
}11.“多项式拟合”按钮控件下添加如下代码
void CJiaoChengDlg::OnBnClickedButton2()
{
// TODO: 在此添加控件通知处理程序代码
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_y,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);
}
12.编译运行，拟合，效果如下
【MATLAB实现】
clc; clear all; close all;
arry1=[-10.4799576938179,-10.4427131117965,-10.3935955278671,-10.3426374470642,-10.2899673769581,-10.2358344513843,-10.1798183736441,-10.1224321068402,-10.0639831025028,-10.0038342606086,-9.94242846041357,-9.87925462149347,-9.81544120185450,-9.75025020464992,-9.68356278430574,-9.61685157675578,-9.54765066199345,-9.47795346582250,-9.40771353803564,-9.33531859017617,-9.26329619851466,-9.18988671334297,-9.11494202978480,-9.04046294563369,-8.96417718978457,-8.88696934573055,-8.80932993147811,-8.73133935992608,-8.65270157047910,-8.57190471187316,-8.49214268506696,-8.41160004102550,-8.33058324409408,-8.24944156448056,-8.16806366584155,-8.08715690190087,-8.00394960062074,-7.92220789617122,-7.84041190849538,-7.75754714348643,-7.67542521642965,-7.59267856421031,-7.51114791596066,-7.42876189264231,-7.34605306709827,-7.26368743725954,-7.18154018461744,-7.09998873001728,-7.01859514528303,-6.93781598917838,-6.85697571266515,-6.77520076634015,-6.69528804142964,-6.61602324004823,-6.53650092857705,-6.45906160005326,-6.37944789272735,-6.30212216338825,-6.22656177121927,-6.14921395947804,-6.07307680182649,-5.99931870060032,-5.92638025774982,-5.85311785355721,-5.78181144672067,-5.71015461474391,-5.64047367200112,-5.57310903359448,-5.50536121942027,-5.44003379673437,-5.37566110675296,-5.31203042530265,-5.25147688370231,-5.19275851059893,-5.13486204325950,-5.07871375008584,-5.02513462588582,-4.97420709775678,-4.92515882779729,-4.87937767370728,-4.83658274580924,-4.79691628400688,-4.76065314560399,-4.72800212324514,-4.70101203366723,-4.67786846033676,-4.66090507222426,-4.65187464054103,-4.64857668595693];
arry2=[94.5486774453437,94.6091950327735,94.6863803006767,94.7635089756550,94.8404442356320,94.9167759923934,94.9930338215959,95.0687009656611,95.1431050512410,95.2171828575088,95.2903318550923,95.3632096312433,95.4345980252129,95.5050702910846,95.5750333666329,95.6427123978068,95.7106758431543,95.7771054939683,95.8418321849418,95.9065005376063,95.9688156893557,96.0301934305458,96.0909289591861,96.1492869455995,96.2070222559737,96.2636020723178,96.3185372735893,96.3718634748319,96.4236645702007,96.4750637610061,96.5239981763049,96.5715019314410,96.6174967295572,96.6616735213237,96.7041828081614,96.7445702765672,96.7842651534391,96.8215358646009,96.8569247901061,96.8910333529152,96.9229949999752,96.9534418537805,96.9816427438242,97.0083214227293,97.0333430048479,97.0564980312684,97.0778307977434,97.0972405538661,97.1148432675382,97.1305528674368,97.1444814112080,97.1568278581921,97.1671840698493,97.1756464845404,97.1824017468487,97.1872169208711,97.1903323575983,97.1917119832175,97.1911851185799,97.1888997725816,97.1848552144178,97.1792118328244,97.1717406206538,97.1624583430699,97.1515834896052,97.1388060921078,97.1245723818418,97.1089216562317,97.0912669369991,97.0723896175601,97.0517896227558,97.0295347379153,97.0063753019352,96.9819354313352,96.9557746422944,96.9283997582167,96.9001946360233,96.8713352376738,96.8413962556383,96.8114152212890,96.7812273352373,96.7512414276797,96.7217213180292,96.6932375286215,96.6678835297399,96.6445093566755,96.6260818062986,96.6155504382452,96.6113860562677];
%%%%%%%%%自写拟合函数
n=6;
TT=max(size(arry1));
X0=zeros(n+1,TT); %构造矩阵X0
for n0=1:n+1
X0(n0,:)=power(arry1(:),(n+1-n0));
end
X=X0';
ANSS=(X'*X)\X'*arry2';
x0=min(arry1):0.1:max(arry1);
y0=zeros(1,length(x0));%根据求得的系数初始化并构造多项式方程
for num=1:1:n+1
y0=y0+ANSS(num)*x0.^(n+1-num);
end
figure,plot(arry1,arry2,'*'),hold on
plot(x0,y0)【MATLAB验证】MATLAB里自带的多项式函数polyfit就是基于最小二乘法实现的
clc; clear all; close all;
arry1=[-10.4799576938179,-10.4427131117965,-10.3935955278671,-10.3426374470642,-10.2899673769581,-10.2358344513843,-10.1798183736441,-10.1224321068402,-10.0639831025028,-10.0038342606086,-9.94242846041357,-9.87925462149347,-9.81544120185450,-9.75025020464992,-9.68356278430574,-9.61685157675578,-9.54765066199345,-9.47795346582250,-9.40771353803564,-9.33531859017617,-9.26329619851466,-9.18988671334297,-9.11494202978480,-9.04046294563369,-8.96417718978457,-8.88696934573055,-8.80932993147811,-8.73133935992608,-8.65270157047910,-8.57190471187316,-8.49214268506696,-8.41160004102550,-8.33058324409408,-8.24944156448056,-8.16806366584155,-8.08715690190087,-8.00394960062074,-7.92220789617122,-7.84041190849538,-7.75754714348643,-7.67542521642965,-7.59267856421031,-7.51114791596066,-7.42876189264231,-7.34605306709827,-7.26368743725954,-7.18154018461744,-7.09998873001728,-7.01859514528303,-6.93781598917838,-6.85697571266515,-6.77520076634015,-6.69528804142964,-6.61602324004823,-6.53650092857705,-6.45906160005326,-6.37944789272735,-6.30212216338825,-6.22656177121927,-6.14921395947804,-6.07307680182649,-5.99931870060032,-5.92638025774982,-5.85311785355721,-5.78181144672067,-5.71015461474391,-5.64047367200112,-5.57310903359448,-5.50536121942027,-5.44003379673437,-5.37566110675296,-5.31203042530265,-5.25147688370231,-5.19275851059893,-5.13486204325950,-5.07871375008584,-5.02513462588582,-4.97420709775678,-4.92515882779729,-4.87937767370728,-4.83658274580924,-4.79691628400688,-4.76065314560399,-4.72800212324514,-4.70101203366723,-4.67786846033676,-4.66090507222426,-4.65187464054103,-4.64857668595693];
arry2=[94.5486774453437,94.6091950327735,94.6863803006767,94.7635089756550,94.8404442356320,94.9167759923934,94.9930338215959,95.0687009656611,95.1431050512410,95.2171828575088,95.2903318550923,95.3632096312433,95.4345980252129,95.5050702910846,95.5750333666329,95.6427123978068,95.7106758431543,95.7771054939683,95.8418321849418,95.9065005376063,95.9688156893557,96.0301934305458,96.0909289591861,96.1492869455995,96.2070222559737,96.2636020723178,96.3185372735893,96.3718634748319,96.4236645702007,96.4750637610061,96.5239981763049,96.5715019314410,96.6174967295572,96.6616735213237,96.7041828081614,96.7445702765672,96.7842651534391,96.8215358646009,96.8569247901061,96.8910333529152,96.9229949999752,96.9534418537805,96.9816427438242,97.0083214227293,97.0333430048479,97.0564980312684,97.0778307977434,97.0972405538661,97.1148432675382,97.1305528674368,97.1444814112080,97.1568278581921,97.1671840698493,97.1756464845404,97.1824017468487,97.1872169208711,97.1903323575983,97.1917119832175,97.1911851185799,97.1888997725816,97.1848552144178,97.1792118328244,97.1717406206538,97.1624583430699,97.1515834896052,97.1388060921078,97.1245723818418,97.1089216562317,97.0912669369991,97.0723896175601,97.0517896227558,97.0295347379153,97.0063753019352,96.9819354313352,96.9557746422944,96.9283997582167,96.9001946360233,96.8713352376738,96.8413962556383,96.8114152212890,96.7812273352373,96.7512414276797,96.7217213180292,96.6932375286215,96.6678835297399,96.6445093566755,96.6260818062986,96.6155504382452,96.6113860562677];
%%%%%%%%%自写拟合函数
n=6;
TT=max(size(arry1));
X0=zeros(n+1,TT); %构造矩阵X0
for n0=1:n+1
X0(n0,:)=power(arry1(:),(n+1-n0));
end
X=X0';
ANSS=(X'*X)\X'*arry2';
x0=min(arry1):0.1:max(arry1);
y0=zeros(1,length(x0));%根据求得的系数初始化并构造多项式方程
for num=1:1:n+1
y0=y0+ANSS(num)*x0.^(n+1-num);
end
figure,plot(arry1,arry2,'*'),hold on
plot(x0,y0)
%%%%%%%%%%MATLAB自带拟合函数
p = polyfit(arry1,arry2,6)
x1 = linspace(min(arry1),max(arry1),100);
y1 = polyval(p,x1);
Newy1=polyval(p,arry1);
SubOfNy1=Newy1-arry2;
SumOfDeta=sum(power(SubOfNy1,2))
figure,plot(arry1,arry2,'*'),hold on
plot(x1,y1,'r')【结论】
以上三种方法得到的结果都是一样的
下载相关文件链接：
****上传的ChartCtrl文件（不能免费上传，需要积分）
本文参考链接：
https://blog.****.net/czyt1988/article/details/21743595