数模利用Matlab和SPSS软件进行线性回归

§3.利用Matlab和SPSS进行线性回归分析回归分析是处理两个及两个以上变量间线性依存关系的统计方法。可以通过软件Matlab和SPSS实现。1.利用Matlab软件实现在Matlab中，可以直接调用命令实现回归分析，（1）[b,bint,r,rint,stats]=regress(y,x)，其中b是回归方程中的参数估计值，bint是b的置信区间，r和rint分别表示残差及残差对应的置信区间。stats包含三个数字，分别是相关系数，F统计量及对应的概率p值。（2）recplot(r,rint)作残差分析图。（3）rstool(x,y)一种交互式方式的句柄命令。以下通过具体的例子来说明。例，现有多个样本的因变量和自变量的数据，下面我们利用Matlab，通过回归分析建立两者之间的回归方程。%一元回归x=[1097128415021394130315551917205121112286231120032435262529483,553372];%因变量时间序列数据y=[698872988807738102513161539156117651762196019022013244627362825];%自变量时间序列数据X=[ones(size(x')),x'],pause[b,bint,r,rint,stats]=regress(y',X,0.05),pause%调用一元回归分析函数rcoplot(r,rint)%画出在置信度区间下误差分布。%多元回归分析%输入各种因变量数据x1=[5.52.58332.98946.55.55653.58647.57]';x2=[3155675038713056427360445039557040506259]';x3=[108127812125851112610106111199]';x4=[86916151781041671264414681311]';%输入自变量数据y=[79.3200.1163.1200.1146.0177.730.9291.9160339.4159.686.3237.5107.2155201.4100.2135.8223.3195]';X=[ones(size(x1)),x1,x2,x3,x4];[b,bint,r,rint,stats]=regress(y,X)%回归分析Q=r'*rsigma=Q/18rcoplot(r,rint);pauseX1=[x1,x2,x3,x4];stepwise(X1,y,[1,2,3])%逐步回归%X2=[ones(size(x1)),x2,x3];%X3=[ones(size(x1)),x1,x2,x3];%X4=[ones(size(x1)),x2,x3,x4];%[b1,b1int,r1,r1int,stats1]=regress(y,X2)%[b2,b2int,r2,r2int,stats2]=regress(y,X3);%[b3,b3int,r3,r3int,stats3]=regress(y,X4);x=[143145146147149150153154155156157158159160162164]';X=[ones(16,1)x];Y=[8885889192939395969897969899100102]';[b,bint,r,rint,stats]=regress(Y,X)rcoplot(r,rint)基于Matlab软件回归及其残差、交互验证分析Matlab

2009-10-0712:05

阅读8

评论0字号：大中小%读取数据

cs=load('d:/lee0921/cl.txt')

%回归分析

y=cs(:,2);

%x=[ones(size(cs(:,3))),cs(:,3)];%包含常数项

x=cs(:,3);

alpha=0.05;

%[b,bint,r,rint]=regress(y,x,alpha)%这里X需要人工加入常数项

%b为回归系数、bint为回归系数置信区间、r为残差、rint为残差置信区间、stats为R^2、F、F显著性概率、误差项方差估计

stats=regstats(y,x);%给出回归常用统计量

%regstats(y,X2)%给出回归常用统计量调出交互式对话框

%绘制散点图及其回归曲线

plot(x,y,'or');%绘制散点图红色

lsline()%添加回归线

%refine(0,3)%添加参考线

title('MainTitle');%添加标题、副标题、坐标轴

xlabel('a','fontsize',12,'fontweight','b')

ylabel('b','fontsize',12,'fontweight','b')

axis([min(x)+5max(x)+5min(y)+5max(y)+5]);%给轴1添加标注和标志控制坐标最大最小值

text(65,70,'\leftarrowheight','FontSize',18);%右箭头为rightarrowyhat=stats.yhat;

sdres=stats.standres;

%绘制多幅图--残差图、条形图、qq图、ACF图

figure;

subplot(2,4,1)

plot(sdres,'ob');

refline(0,0);

axis([0length(y)min(sdres)+0.25max(sdres)+0.25]);

xlabel('index');ylabel('Residuals');

%残差条形图、密度曲线估计

subplot(2,4,2)

%hist(sdres);

%holdon;

%[f,xi]=ksdensity(sdres);%采用核密度估计

%plot(xi,f);

%holdoff;xmax=1.1*max(sdres);

xmin=1.1*min(sdres)

bins=floor(xmin):0.4:ceil(xmax);

h=bar(bins,histc(sdres,bins)/length(sdres),'histc');

%h=bar(bins,histc(sdres,bins),'histc');

set(h,'FaceColor',[.9.9.9]);

xlabel('BlockMaximum');ylabel('ProbabilityDensity');

xlim([xmin-1xmax+1]);

holdon;

[f,xi]=ksdensity(sdres);%采用核密度估计

line(xi,f);

holdoff;

%待落实关于直方图和密度曲线绘制

%画QQ图

subplot(2,4,3)

qqplot(sdres);

title('');

boxon;

%自相关图

subplot(2,4,4)

[ACF,Lags,Bounds]=autocorr(sdres,[],2);%计算自相关系数值

autocorr(y,[],2)%绘制自相关图

title('');

subplot(2,4,5)

plot(sdres,'ob');

refline(0,0);

axis([0length(y)min(sdres)+0.25max(sdres)+0.25]);

xlabel('index');ylabel('Residuals');

%残差条形图、密度曲线估计

subplot(2,4,6)

%hist(sdres);

%holdon;

%[f,xi]=ksdensity(sdres);%采用核密度估计

%plot(xi,f);

%holdoff;xmax=1.1*max(sdres);

xmin=1.1*min(sdres)

bins=floor(xmin):0.4:ceil(xmax);

h=bar(bins,histc(sdres,bins)/length(sdres),'histc');

%h=bar(bins,histc(sdres,bins),'histc');

set(h,'FaceColor',[.9.9.9]);

xlabel('BlockMaximum');ylabel('ProbabilityDensity');

xlim([xmin-1xmax+1]);

holdon;

[f,xi]=ksdensity(sdres);%采用核密度估计

line(xi,f);

holdoff;%待落实关于直方图和密度曲线绘制

%画QQ图

subplot(2,4,7)

qqplot(sdres);

title('');

boxon;

%自相关图

subplot(2,4,8)

[ACF,Lags,Bounds]=autocorr(sdres,[],2);%计算自相关系数值

autocorr(y,[],2)%绘制自相关图

title('');

%交互验证--crossvalidation并绘制散点图1与折线图2

cv=diag([ones(length(x),1),x]*stats.beta_i);

%图1--散点图plot(cv,y,'o');

lsline;

xlabel('Crossvalidation');ylabel('Observeddata');%图2--折线图：观测值、交互验证折线图，分别用实线、虚线给出图形

figure

plot(y,'-b')%红色实线,坐标设置

xlabel('Year');ylabel('Observeddata');

holdon;

plot(cv,'--r');%蓝色虚线

legend('Observed','Crossvalidated',2);

legend('boxoff')%不绘制盒子

holdoff;%采用核密度估计绘制密度曲线

%figure

%xmax=1.1*max(sdres);

%xmin=1.1*min(sdres)

%bins=floor(xmin):0.4:ceil(xmax);

%h=bar(bins,histc(sdres,bins)/length(sdres),'histc');

%%h=bar(bins,histc(sdres,bins),'histc');

%set(h,'FaceColor',[.9.9.9]);

%xlabel('BlockMaximum');ylabel('ProbabilityDensity');

%xlim([xmin-1xmax+1]);

%holdon;

%[f,xi]=ksdensity(sdres);%采用核密度估计

%line(xi,f);

%holdoff;

数学实验十一回归分析实验目的：1．

了解回归分析的基本原理，掌握MATLAB的实现方法；2．

联系用回归分析方法解决实际问题。实验内容：1）用切削机床加工时，为实时地调整机床需测定刀具的磨损程度，先每隔一小时测量刀具的厚度得到以下的数据，试建立刀具厚度关于切削时间的回归模型，对模型和回归系数进行检验，预测15小时后刀具的厚度。计算分析：首先对原始数据进行观察，确定回归模型，然后通过计算最终确定模型和模型参数，并对模型和回归系数进行检验。在进行回归分析之前，现将刀具厚度和时间的关系绘图，程序如下：x1=0:1:10;y=[30.6 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5]';plot(x1,y,'+')结果输出（右图）：从原始数据上看，可以建立线性回归模型。下面计算回归系数：alpha=0.05;x=[ones(11,1),x1'];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)结果输出：b=29.5455-0.3291bint=28.976930.1140-0.4252-0.2330stats=0.869660.00180.0000因为p<alpha，拒绝，回归模型成立。观察残差分布，发现第一个数据（0，30.6）残差的置信区间不包括零点，应视为异常点。将其剔除后重新计算，可得：alpha=0.05;x1=1:1:10;y=[29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5]';x=[ones(10,1),x1'];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)

b=29.0533-0.2588bint=28.833429.2732-0.2942-0.2233stats=0.9726283.55990.0000再次计算，发现原始数据中的第二个数据（1，29.1）残差的置信区间也不包括零点，仍该视为异常点。剔除再次计算，可得：alpha=0.05;x1=2:1:10;y=[28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5]';x=[ones(9,1),x1'];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)

b=28.8667-0.2333bint=28.779628.9537-0.2467-0.2200stats=1.0e+003*0.00101.71500.0000这次的数据残差的置信区间全部包括零点，无异常点。对比分析：对比剔除前后的变化，发现置信区间明显缩小，相关系数和F都明显变大。表明异常点的剔除有利于更好的建立模型。输出最终计算结果和图形（其中b为最后输出是的b）：x1=0:1:10;y=[30.6 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5]';x0=0:0.1:12;y0=b(1)+b(2)*x0;plot(x1,y,'+',x0,y0)y15=b(1)+b(2)*15得：y15=25.3667

回归模型:预测结果：15小时后，刀具的厚度将为25.3667cm。结果分析：1．模型检验与异常点剔除不能混为一谈。从两次剔除异常点前后p值没有明显变化（准确的说是一开始的p值就表明了这是一个线性回归模型）可以看出，异常点的存在并不能从根本上动摇和改变合理的模型，而不合理的模型也不会通过几个异常点的剔除就变得合理的。从显著性水平也可以证明这一点，即残差置信区间的显著性水平和模型检验的显著性水平并没有关系，是相互独立的。2．异常点剔除次序。这也是本题目带来的一个问题，两个异常点不是一次剔除的。在剔除第一个点之前，第二个点还是合理的，但剔除后就不合理了。这是可以给出合理解释的，从原始数据会出的图上可以看出，第一二个点都是不合理的，但由于它们都在开始的位置，影响回归模型发生的变化使得第二个点恰好在回归模型上；而第一个点被剔除后，第二个点再没有点来屏蔽，其异常性就显露出来。这也就是为什么我们看到残差区间图中第二点并不像是异常点的原因。3．从上一个讨论可以得到一个思路，就是在异常点剔除时，首先要对原始数据加以分析，从中选出可能是异常的点，伴随残差置信区间图一起分析，才能避免分析过程中可能出现的片面性。4．对于异常点的解释：这可能是刀具还没有进入稳定工作状态，致使开始时磨损过快。5．模型选择的考虑：本题目选择了线性模型，这是通过观察数据得到的，而通过对数据关系的分析也可以得到这一点。在刀具没有大幅度磨损之前，磨损速度应该是一定的，可以建立起线性模型，这是客观决定的。⑶得到某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型并进行检验，预测平均收入为1000、价格为6时的商品需求量。需求量10075807050659010011060收入10006001200500300400130011001300300价格5766875439计算分析：观察原始数据，试分析回归模型x1=[1000 600 1200 500 300 400 1300 1100 1300 300];x2=[5 7 6 6 8 7 5 4 3 9];y=[100 75 80 70 50 65 90 100 110 60];plot(x1,y,'+')pauseplot(x2,y,'+')

从数据上看，y与x1、x2均有线性关系，与x1的线性关系更明显一些，可以计算一下：alpha=0.05;x1=[1000 600 1200 500 300 400 1300 1100 1300 300];x2=[5 7 6 6 8 7 5 4 3 9];y=[100 75 80 70 50 65 90 100 110 60]';x=[ones(10,1),x1',x2'][b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats结果输出：b=111.69180.0143-7.1882bint=56.0503167.3334-0.01200.0406-13.2306-1.1458stats=0.894429.65330.0004从结果看，这个线性回归模型可以接受：数据有较强的线性相关性（r=0.8944），同时也能通过一般程度的假设检验（p=0.0004），没有异常点。当然也不可否认，这个模型并不理想，第二个参数包括了零点，可以考虑用二次函数改进。二次函数改进：x1=[1000 600 1200 500 300 400 1300 1100 1300 300];x2=[5 7 6 6 8 7 5 4 3 9];y=[100 75 80 70 50 65 90 100 110 60]';x=[x1',x2']rstool(x,y)beta,rmse,residuals①Linear（线性）：beta=111.69180.0143-7.1882rmse=7.2133residuals=9.95235.0477-5.7188-5.7109-8.4750-2.0929-4.33681.33441.28678.7133②PureQuadratic（纯二次）：beta=110.53130.1464-26.5709-0.00011.8475rmse=4.5362residuals=5.2724-0.7162-4.5158-1.9390-3.33153.45663.4843-3.4452-0.09761.8320③Interaction（交叉）：beta=96.99750.0304-5.0692-0.0026rmse=7.5862residuals=10.89816.1178-4.4214-4.0155-9.3533-1.4260-4.33701.2410-1.19496.4911④FullQuadratic（完全二次）：beta=-106.60950.326121.2990-0.0200-0.0001-0.7609rmse=4.4179residuals=3.5434-2.0192-1.97960.8256-4.92303.08911.9182-3.45020.27502.7207《计量地理学》（徐建华，高等教育出版社，2005）配套实习指导PAGEPAGE19剩余标准差对比：线性纯二次交叉完全二次7.21334.53627.58624.4179从比较结果可得，完全二次模型的剩余标准差最小，应该建立完全二次模型。其参数为：beta=-106.60950.326121.2990-0.0200-0.0001-0.7609回归模型：预测结果：平均收入为1000，价格为6时的商品需求量为89.3601。结果分析：1．模型选择。本题目的重点在于选择最合适的模型，之所以不说合理而说合适，是因为在本题中，几种模型都很合理，但与原始数据的符合性上还要找出最合适的。首先是线性模型，正如前面分析，线性模型虽然不是最合适的，但在很大程度上符合了原始数据，这表明了商品需求量与收入和价格主要体现了线性关系，在要求不高时完全可以用线性模型近似。但因为线性模型忽略了二次关系和相互关系，在较高要求的回归分析中是不够的，这就需要进一步用多项式回归来确定。计算表明，完全二次模型是最合适的。2．对模型的分析：模型的系数表明，需求量与收入和价格很大程度上是线性关系，同时二次项将消除一部分正面的增长，交叉项表明两者之间有相互影响关系。3．疑惑：按照上述分析，线性关系应该是主要关系，那么几种模型中得线性系数应该不会有太大变化，但事实上存在很大变化，甚至会出现符号的变化。这是很不能让人理解的，因为这将表明线性关系不是合理的关系，只是一种巧合使得数据看起来有一些线性关系。如果这样的话，讨论1和2便没有说服力了。⑺一矿脉有13个相邻样本点，人为地设定一原点，现测得各样本点对原点的距离x，与该样本点处某种金属含量y的一组数据如下，画出散点图观察二者的关系，试建立合适的回归模型，如二次曲线、双曲线、对数曲线等。x23457810y106.42109.20109.58109.50110.00109.93110.49x111415151819

y110.59110.60110.90110.76111.00111.20

计算分析：观察原始数据绘图。x1=[2 3 4 5 7 8 10 11 14 15 15 18 19]';y=[106.42 109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';plot(x1,y,'+')从图像上看，如果第一个点数据剔除，线性关系比较明显，但并不能排除其他模型。下面就对几种模型都加以计算比较。①线性模型计算程序：alpha=0.05;x1=[2 3 4 5 7 8 10 11 14 15 15 18 19]';y=[106.42 109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';x=[ones(13,1),x1];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)结果输出：b=108.25810.1742bint=107.2794109.23670.08910.2593stats=0.648420.28660.0009线性相关系数较小，线性回归模型在alpha>0.0009成立

第一个点为异常点（仅指线性模型下），予以剔除，再计算得：计算程序：alpha=0.05;x1=[3 4 5 7 8 10 11 14 15 15 18 19]';y=[109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';x=[ones(12,1),x1];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)结果输出：b=109.06680.1159bint=108.8264109.30720.09580.1360stats=0.9428164.80600.0000剔除第一个点后线性系数和p值都变得好了很多。没有异常点。线性模型：对该模型求剩余标准差：rmse=sqrt(sum((y-b(1)-b(2)*x1).^2)/10)得：rmse=0.1635②二次曲线：考虑第一个点偏离太多，予以剔除后重新计算可得（此乃参考几个模型所得之结论）：计算程序：x=[3 4 5 7 8 10 11 14 15 15 18 19];y=[109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20];[p,S]=polyfit(x,y,2);p结果输出：p=-0.00430.2102108.6718二次模型：对该模型求剩余标准差：[Y,delta]=polyconf(p,x,S);rmse=sqrt(sum((y-Y).^2)./10),得：rmse=0.1231③双曲线模型：双曲线模型类似于，可以通过将x的倒数代换转化为线性模型来求。计算程序：alpha=0.05;x1=[2 3 4 5 7 8 10 11 14 15 15 18 19]';y=[106.42 109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';x=[ones(13,1),1./x1];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)结果输出：b=111.4405-9.0300bint=111.1068111.7743-10.6711-7.3889stats=0.9302146.67330.0000有两个异常点，剔除后可得（第一次计算后又剔除了一个点，这里不予列出，只写最后的程序）：alpha=0.04;x1=[5 7 8 10 11 14 15 15 18 19]';y=[109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';x=[ones(10,1),1./x1];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)

b=111.5653-10.9938bint=111.2882111.8424-13.5873-8.4002stats=0.9309107.76230.0000双曲线模型：对该模型求剩余标准差：rmse=sqrt(sum((y-b(1)-b(2)./x1).^2)/8)得：rmse=0.1487④对数曲线：类似于双曲线模型：计算程序：alpha=0.05;x1=[2 3 4 5 7 8 10 11 14 15 15 18 19]';y=[106.42 109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';x=[ones(13,1),log(x1)];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)结果输出：b=106.71131.5663bint=105.6382107.78441.08282.0499stats=0.822150.82850.0000剔除异常点，重新计算可得：alpha=0.05;x1=[3 4 5 7 8 10 11 14 15 15 18 19]';y=[109.20 109.58 109.50 110.00 109.93 110.49 110.59 110.60 110.90 110.76 111.00 111.20]';x=[ones(12,1),log(x1)];[b,bint,r,rint,stats]=regress(y,x,alpha);b,bint,stats,rcoplot(r,rint)结果输出：b=107.97621.0496bint=107.6403108.31210.90371.1956stats=0.9625256.70140.0000对数模型：对该模型求剩余标准差：rmse=sqrt(sum((y-b(1)-b(2)*log(x1)).^2)/10)得：rmse=0.1324结果比较：通过对几个模型的比较可得，二次模型的剩余标准差最小，可以采用。不过几个模型的差别很小。线性模型二次模型双曲线模型对数模型0.16350.12310.14870.1324结果分析：1．对第一个点的讨论。纵观四个模型，第一个点都属于异常点，需要剔除。像这样的剔除是否合理呢？是不是为了模型的美观就可以剔除用不上的点呢？对于这个问题，我认为不合理的点固然要剔除，但同时点数的减少又将使得样本的容量变小，信度降低，这就需要使用者的判断。向本题中的第一个数据，很明显不符合任何模型，严重干扰回归分析，可以判断为是异常点，予以剔除。2．模型的选择。这已经是一个老话题了，前两道题都提到了这一点。本题目的特点在于，因为对矿物分布和地质知识的缺乏，不能从理论上加以分析，只能从数据本身出发，加以分析。这就隐藏了很多问题，因为没有理论支持，这几种模型都是可行的（无非都是对纯粹的数据进行回归），而是否合理并没有判断标准，事实上，符合得最好的模型将会是插值模型，但这就失2.利用SPSS软件实现在SPSS软件中，同样可以简单的实现回归分析，因为回归分析包含了线性回归与曲线拟合两部分内容，首先来看线性回归分析过程（LINEAR）。一元回归分析例如，教材中上表3.1.2中数据，把降水量（P）看作因变量，把纬度（Y）看作自变量，在平面直角坐标系中作出散点图（课本上图3.1.1），发现它们之间呈线性相关关系，因此，可以用一元线性回归方程近似地描述它们之间的数量关系。步骤如下：在菜单中选择Regression==>liner，系统弹出线性回归对话框如下：各项选项按钮解释如下：【Dependent框】用于选入回归分析的应变量。即在这个例子中选择年降水量。【Block按钮组】由Previous和Next两个按钮组成，用于将下面Independent框中选入的自变量分组。由于多元回归分析中自变量的选入方式有前进、后退、逐步等方法，如果对不同的自变量选入的方法不同，则用该按钮组将自变量分组选入即可。下面的例子会讲解其用法。【Independent框】用于选入回归分析的自变量。即在这个例子中选择纬度。【Method下拉列表】用于选择对自变量的选入方法，有Enter（强行进入法）、Stepwise（逐步法）、Remove（强制剔除法）、Backward（向后法）、Forward（向前法）五种。该选项对当前Independent框中的所有变量均有效。【SelectionVariable框】选入一个筛选变量，并利用右侧的Rules按钮建立一个选择条件，这样，只有满足该条件的记录才会进入回归分析。【CaseLabels框】选择一个变量，其取值将作为每条记录的标签。最典型的情况是使用记录ID号的变量。【WLS>>钮】可利用该按钮进行权重最小二乘法的回归分析。单击该按钮会扩展当前对话框，出现WLSWeight框，在该框内选入权重变量即可。【Statistics钮】弹出Statistics对话框，用于选择所需要的描述统计量。有如下选项：RegressionCoefficients复选框组：定义回归系数的输出情况，选中Estimates可输出回归系数B及其标准误，t值和p值，还有标准化的回归系数Beta；选中ConfidenceIntervals则输出每个回归系数的95%可信区间；选中CovarianceMatrix则会输出各个自变量的相关矩阵和方差、协方差矩阵。以上选项默认只选中Estimates。Residuals复选框组：用于选择输出残差诊断的信息，可选的有Durbin-Watson残差序列相关性检验、超出规定的n倍标准误的残差列表。ModelFit复选框：模型拟合过程中进入、退出的变量