参考链接:
【MATLAB】MATLAB 线性拟合小结 —— polyfit_vera吧噗_新浪博客
数据拟合与插值方法_tracezl的博客-CSDN博客_数据插值方法
【MATLAB】MATLAB 线性拟合小结 —— REGRESS多元线性回归(用最小二乘估计法)_vera吧噗_新浪博客
1多元线性回归:regress
函数表达形式为 y=β0+β1x+ε (ε是随机误差,满足E(ε)=0,var(ε)=σ2)
- [b,bint,r,rint,stats]= regress(y,X,alpha) 调用regress函数作一元线性回归
- 输入X为自变量,n行p列。X应该包含一个全“1”的列,这样则该模型包含常数项。
否则会警告没有P值:
R-square and the F statistic are not well-defined unless X has a column of ones.
利用X=[ones(size(x,1)) ,x]实现。y为因变量,n行1列。alpha为估算置信区间的精度。
-
- b为回归系数向量,从低到高的列向量,而polyfit中为从高到低。yhat = xdata*b; % 计算y的估计值
- bint:b的95%置信区间矩阵,置信区间不大,说明有效性较好;若含零点,说明结果无效。
- r:残差向量(因变量的真实值减去估计值)
- rint:返回残差的95%置信区间,它是一个2×n的矩阵,第1列为置信下限,第2列为置信上限。该矩阵可以用来诊断异常(即发现奇异观测值)。如果第i组观测的残差的置信区间RINT(i,:)所定区间没有包含0,则第i个残差在默认的5%的显著性水平比我们所预期的要大,这可说明第i个观测值是个奇异点(即说明该点可能是错误而无意义的,如记录错误等)
可用rcoplot(r,rint)画残差,红色的表示超出期望值的数据。圆圈代表残差的值,竖线代表置信区间的范围。
- stats:向量,STATS中的4个值分别为:R2(判定系数),F(总模型的F测验值),P(总模型F的概率值P(F>Fz)),MSq(离回归方差或误差方差的估计值)。
- 判定系数(the Coefficient of the Determination)R2:是判断回归模型拟合程度的一个指标,其取值范围为[0, 1];判定系数越大说明回归模型的拟合程度越高,回归方程越显著。
- F>F(1-α)(k, n-k-1)时拒绝H0,F越大,说明回归方程越显著。
与F对应的概率P<α时拒绝H0,回归模型成立。
- MSq:由于最小二乘法中不求误差方差σ2,其误差平方和Msq定义为SSR/自由度,其中SSR为Regression Sum of Squares。
- stats = reglm(y,X,model,varnames) 多重线性回归分析
- 产生线性回归分析的方差分析表和参数估计结果,并以表格形式显示在屏幕上。参数X是自变量观测值矩阵,它是n行p列的矩阵。y是因变量观测值向量,它是n行1列的列向量
- stats:返回一个包括了回归分析的所有诊断统计量的结构体变量stats。
- 可选的model参数来控制回归模型的类型
- 'linear' 带有常数项的线性模型(默认情况)
- 'interaction' 带有常数项、线性项和交叉项的模型
- 'quadratic' 带有常数项、线性项、交叉项和平方项的模型
- 'purequadratic' 带有常数项、线性项和平方项的模型
eg:
scatter(x,y,'k.');
hold on
[xx,yy]=prepareCurveData(x,y); %去除空值
xx=[ones(length(xx),1),xx]; %添加常数项
[b,bint,r,rint,stats]= regress(yy,xx,0.05); %x自变量,y因变量,置信水平0.05
% rcoplot(r,rint); %残差图,红色为超出期望值的数据
y2=b(1)+b(2)*x; %拟合线
plot(x,y2,'r-');
text(2,15,['y=',num2str(roundn(b(2),-2)),'*x+',num2str(roundn(b(1),-2))])
if stats(3)<0.01
str='^{**}';
elseif stats(3)<0.05
str='^*';
else
str=[];
end
text(2,14,['R^2=',num2str(roundn(stats(1),-2)),str])
text(2,13,['P=',num2str(stats(3))])
axis([0 16 0 16])
legend('wind speed')
xlabel('lac')
ylabel('PL')
box on
hold off
2多项式拟合:polyfit线性最小二乘拟合
函数表达形式为:f(x)=a1*xm+...+am*x+am+1
- 函数命令为:[a,s]=polyfit(x,y,m)
- x,y为对应的自变量,m为需要拟合的最高次幂,产生一个m次的多项式
- 矩阵S用于生成预测值的误差估计。 S是一个结构体数组(struct),用来估计预测误差,包含了R,df和normr。
- R:polyfit函数中,先根据输入的x构建范德蒙矩阵V,然后进行QR分解,得到的上三角矩阵。
- df:自由度, df=length(y)-(n+1)。df>0时,为超定方程组的求解,即拟合点数比未知数(p(1)~p(n+1))多。
- normr:标准偏差、残差范数,normr=norm(y-V*p),此处的p为求解之后的数值。
- [a,s,mu]=polyfit(x,y,m)
mu=[mean(X); std(X)],mean(X)求X每一列的均值,std(X)求X的标准差。
此处加了mu会改变拟合结果,未知原因。建议不加mu。
- 相关系数R2:rsquare=1-s.normr^2/norm(y-mean(y))^2
- y=polyval(a,x); 根据拟合的函数得出x对应的因变量的值
- [y,delta]=polyconf(a,x,s,alpha) 对ployfit拟合曲线进行评价和置信区间估计
计算回归多项式在x处的预测值y及预测值的显著性为1-alpha的置信区间Y±DELTA,alpha缺省时为0.05
- eg:
- 多项式拟合polyfit
scatter(ts,ta,'.')
hold on
[xx,yy]=prepareCurveData(ts,ta); %去除空值
[a1,b1]=polyfit(xx,yy,1); %拟合
x1=0:1:30;
y1=polyval(a1,x1);
plot(x1,y1,'r-');
rsquare=1-b1.normr^2/norm(yy-mean(yy))^2; %求R2
text(2,15,['y=',num2str(roundn(a1(1),-2)),'*x+',num2str(roundn(a1(2),-2))])
text(2,12,['R2=',num2str(roundn(rsquare,-2))])
axis([0 20 0 20])
xlabel('ts')
ylabel('ta')
box on
hold off
- 绘制置信区间
x=0:30;
y=xlsread('C:\Users\Desktop\data.xlsx');
[p,s]=polyfit(x,y,4);
y1= polyval(p,x);
[yfit,dy] = polyconf(p,x,s,'predopt','curve');
%绘制置信范围
xconf=[x,fliplr(x)]; %一个来回
yconf=[yfit-dy,fliplr(yfit+dy)] %条带宽度delta
h1=fill(xconf,yconf,'r','FaceColor',[1 0.8 0.8],'EdgeColor','none'); %设定颜色
hold on
h2=plot(x,y1,'b','linewidth',2);
h3=plot(x,y,'r.','markersize',15);
点状置信区间plot(x,y,'k+',x,Y,'r * ',x,yy+DELTA,'r o ',x,yy-DELTA,'r o ')
3稳健回归:robust
- [b,stats]=robustfit(x,y,’wfun’,tune,’const’) 受异常值的影响较小
- b返回系数估计向量;stats返回各种参数估计;’wfun’指定一个加权函数;tune为调协常数;’const’的值为’on’(默认值)时添加一个常数项;为’off ’时忽略常数项。
4自定义函数拟合 nlinfit()或cfun
- [beta,r,J]=nlinfit(X,y,model,beta0) 多元非线性回归
- beta返回函数model中的待定常数;r表示残差;J表示雅可比矩阵。X,y为数据;model自定义函数;beta0待定常数初值。
eg:
mymodel = inline('beta(1)*x(:,1).^3 + beta(2)*x(:,2).^2 + beta(3)*log2(x(:,3)','beta','x');
%x(:1)是取x的第一列数据,相应的x(:,2)取第二列数据
%beta(1),beta(2),beta(3)为所求系数
x = [x1,x2,x3]; %自变量x1,x2,x3
beta0= [a1,b1,c1,d1,e1,f1,m1,n1,z1]; %系数初始值
beta = nlinfit(x,y,myfun,temp); %回归得到系数beta
- fit 自定义拟合函数
eg:
f(t)=a*t+b
f=fittype('a*t+b','independent','t','coefficients',{'a','b'}); %fittype是自定义拟合函数
cfun=fit(b(:,1),b(:,2),f) %根据自定义拟合函数f来拟合数据x,y
x1=0:1:30;
y1=cfun(x1);
plot(x1,y1,'b-');
5曲线拟合工具箱cftool
直接在matlab命令行输入cftool命令即可进入cftool窗口,输入命令回车后就得到界面
进行曲线拟合的时候,最基本的操作包括点击“data...”按钮和“fitting...”按钮。
2.1 data...按钮
首先,到命令窗口中为点击"data..."做好准备:
- >> x=[8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20];
- >> y=[0.6,0.62,0.64, 0.65, 0.66, 0.67, 0.68, 0.68, 0.69, 0.66, 0.65, 0.65,0.64];
- >>
然后,点击“data...”按钮,选x横坐标和选y纵坐标变量:
1)在弹出的Data窗口中,在X Data选项中选择x,在Y Data中选择y。其中X, Y都是刚刚在命令行所准备的数据变量名。
2)经过此步骤后点击“Create data set”按钮。这个时候Data对话框背后的对话框中已经有散点图了。
不要害怕没有出结果,请果断的点击Data中的close按钮吧。
2.2 fitting...按钮
在正确设置Data对话框之后,在关闭Data对话框之后,就可以点击fitting...按钮了。会弹出这么一个对话框:
1)在Fitting对话框中点击New fit后在Fit name为此次曲线取名字:“cftool 曲线拟合”,当有多条曲线需要同时绘制时,每次都需要点击New fit,然后选择对应的Data set值即对应相应的数据变量就可以绘制多条曲线了。即同时拟合多条曲线。选择data set。
2)选择相应类型的曲线拟合(到其它地方粘贴过来的),点击Type of fit 下拉菜单:
Custom Equations:用户自定义的函数类型
· Exponential:指数逼近,有2种类型, a*exp(b*x) 、 a*exp(b*x) + c*exp(d*x)
· Fourier:傅立叶逼近,有7种类型,基础型是 a0 + a1*cos(x*w) + b1*sin(x*w)
· Gaussian:高斯逼近,有8种类型,基础型是 a1*exp(c1-((x-b1)/c1)^2)
· Interpolant:插值逼近,有4种类型,linear、nearest neighbor、cubic spline、shape-preserving
· Polynomial:多形式逼近,有9种类型,linear ~、quadratic ~、cubic ~、4-9th degree ~
· Power:幂逼近,有2种类型,a*x^b 、a*x^b + c
· Rational:有理数逼近,分子、分母共有的类型是linear ~、quadratic ~、cubic ~、4-5th degree ~;此外,分子还包括constant型
· Smoothing Spline:平滑逼近(翻译的不大恰当,不好意思)
· Sum of Sin Functions:正弦曲线逼近,有8种类型,基础型是 a1*sin(b1*x + c1)
· Weibull:只有一种,a*b*x^(b-1)*exp(-a*x^b)
3)最后再选择一种曲线拟合下的具体类型,点击Apply按钮就可以看见所拟合的曲线了。
如果发现曲线和数据点不是很“磨合”就需要重新选择曲线类型了。
上图是点击了应用按钮之后的状态,在Results中得来的就是拟合数据的曲线函数关系式。
点击Appy后的运行结果为:
2.3 CFTOOL菜单使用
在上面显示的曲线拟合结果中,可以明显的看见横坐标和纵坐标都不是从原点(0,0)开始设置的坐标。那么我们可以使用
Curve Fitting Tool 对话框中Tools菜单中的
然后Tools中的各个选项都有什么特殊的功能呢,呵呵,亲,你随便点点用用就知道了:
1)New Cunstom Equation:表示使用自定义的函数来模拟图中显示的散点数据的函数式子。
2)Legend:就是图中显示着“y vs.x cftool 曲线拟合”的那个框。
3)Grid:跟指令grid一样的功能,是曲线背景充满虚线框。
4)Axis Limit Control:设置横轴、纵轴坐标的范围,这个按钮就可以使横纵坐标从你想要的初值到终值了。