MATLAB回归分析

来源：互联网发布：苹果7mac地址编辑：程序博客网时间：2024/06/05 21:54

数据预处理

1.计算每年每一个变量的基本统计量：均值，标准差，中位数；并用折线图给出原始数据和均值、中位数的变化趋势。

由所给的《中国企业商品价格指数数据99年至今》数据，利用MATLAB，很容易计算出结果，并画出折线图。

折线图如下：

2.给出各变量按季节变化趋势；每年各变量频数直方图的变化趋势，可做动态图。

利用MATLAB作图如下：

（此频数直方图是动态图）

3. 按原始数据，指出被解释变量和每解释变量之间的关系，画图；

利用MATLAB分别作图如下：

从图中可以看出，随着农产品价格指数的增长，总指数有线性增长的趋势，所以总指数与农产品价格指数有线性关系。

从图中可以看出，随着矿产品价格指数的增长，总指数有线性增长的趋势，所以总指数与矿产品价格指数有线性关系。

从图中可以看出，随着煤油电价格指数的增长，总指数有线性增长的趋势，所以总指数与煤油电价格指数有线性关系。

4.按年平均数、中位数预处理数据，指出被解释变量和每个解释变量之间的关系，画图；

利用MATLAB，分别作图如下：

从图中可以看出，对于每个解释变量（各个价格指数中均值），随着其增长，被解释变量（总指数均值）都有线性增长的趋势，所以被解释变量与每个解释变量分别有线性关系。

从图中可以看出，对于每个解释变量（各个价格指数均中位数），随着其增长，被解释变量（总指数中位数）都有线性增长的趋势，所以被解释变量与每个解释变量分别有线性关系。

（本节MATLAB代码见附录一）

一元回归分析

基本模型

下面先从简单入手，由上节数据预处理第三问求解得出，总指数y与农产品价格指数x之间存在线性关系，所以，可以建立以下一元线性回归模型：

y _i= β₀+β₁x_i+ε_i , i=1,2……，191

各ε_i独立同分布，其分布为N（0，σ²）

由数据（xi，yi）(i=1,2……，191)可获得β₀、β₁的估计称

为y关于x的回归方程，，为回归系数，

εi是随机误差。

模型求解

利用MATLAB统计工具箱，得到模型的回归系数估计值及其置信区间（置信水平α=0.05）、检验统计量R²，F，p，s²。

回归系数

回归系数估计值

43.5577

[39.1754 ,47.9401]

0.5606

[ 0.5184 ,0.6028]

R²= 0.7839 F=685.6542 p<0.0001 s² =4.4997

残差图如下：

剔除异常值后，再对模型进行求解，结果见下表：

回归系数

回归系数估计值

46.6171

[43.0129,50.2213]

0.5323

[0.4976,0.5670]

R²= 0.8391 F=917.7627 p<0.0001 s² =2.9046

回归方程为：= 46.6171+0.5323x

结果分析

R²=0.8391是指因变量y的83.91%可由模型确定。F=917.7627远大于临界值，p远小于α，说明回归方程显著。

同时，回归系数置信区间不包含0，说明回归变量对因变量y的影响是显著的。

综上，该模型是可用的。

(本节代码见附录二)

多元回归分析

基本模型

由上节数据预处理第三问求解得出，总指数y与农产品价格指数x1、矿产品价格指数x2,煤油电价格指数x3之间均存在线性关系，故而，建立如下多元线性回归模型：

yi= β₀+β₁X1i+β₂X2i+β₃X3i+ε_i , i=1,2……，191

各ε_i独立同分布，其分布为N（0，σ²）

由数据（x1i，x2i，x3i，yi）(i=1,2……，191)可获得β₀、β₁β₂、β₃的估计称

为y关于x1,x2,x3的回归方程,，为回归系数,εi是随机误差。

模型求解

利用MATLAB统计工具箱，得到模型的回归系数估计值及其置信区间（置信水平α=0.05）、检验统计量R²，F，p，s²。结果见下表：

回归系数

回归系数估计值

回归系数置信区间

32.5824

[30.6793,34.4855]

0.3742

[0.3485,0.4001]

0.1372

[0.1087,0.1656]

0.1503

[0.1312,0.1693]

R² = 1.0000 F=2008.4 p<0.0001 s² = 0.6000

残差图如下：

剔除异常值后，再对模型进行求解，结果见下表：

回归系数

回归系数估计值

回归系数置信区间

32.6886

[30.9173,34.4598]

0.3771

[0.3532,0.4010]

0.1313

[0.1051,0.1575]

0.1524

[0.1347,0.1701]

R²= 1.0000 F=2291.3 p<0.0001 s² = 0.5000

回归方程：

= 32.6886+0.3771x1+0.1313x2+0.1524x3

结果分析

R²近似等于1，说明回归方程非常显著。F=2291.3远远大于临界值，p远远小于α，故而同样说明回归方程显著。

同时，回归系数置信区间均不包含0，说明各个回归变量对因变量y的影响是显著的。

综上，该模型是非常好的。

（本节代码见附录三）

时间序列模型

问题分析

题目所给的原始数据是以时间为序的，称时间序列。许多经济数据在时间上有一定的滞后性，实际上，在对时间序列做回归分析时，随机误差实际上，在对时间序列做回归分析时，随机误差ε_i有可能存在相关性，违背模型关于ε_i（对时间）相互独立的基本假设，同一变量的顺序观测值之间出现相关现象（称为自相关）是很自然的。然而一旦数据中存在这种自相关序列，如果仍采用普通的回归模型直接处理，将会出现不良结果，其预测也就失去了意义。为此，我们必须先来诊断数据是否存在自相关，如果存在，就要考虑自相关关系，建立新的回归模型。

残差e_i = yi -可以作为εi估计值，画出e_i– e_i-1的散点图，能够从直观上判断ε_i的自相关性。由前面多元线性回归模型得到残差e_i，作出e_i– e_i-1的散点图

如下：

可以看到：大部分点都落到了第一、第三象限中，表明ε_i存在正的自相关。

为了对ε_i的自相关性作定量诊断，并在诊断后得到新的结果，我们考虑如下的模型

yi= β₀+β₁X1i+β₂X2i+β₃X3i+ε_i,ε_i = ρε_i-1+u_i

其中，ρ是相关系数，u_i相互独立且服从均值为0的正态分布。

若ρ= 0，则退化为普通的回归模型；若ρ>0，则ε_i存在正的自相关；若ρ<0，则ε_i存在负的自相关。

下面做定量检验，即Durbin-Watson检验（简称D-W检验），计算DW统计量：

DW =

经过简单的运算可知，当n较大时，

正是自相关系数ρ的估计值，于是 = 1 – DW/2。若在1附近，即DW在0或4附近，ε_i的自相关性很强。

基本模型

要根据DW的具体数值确定ε_i是否存在自相关，应该在给定的检验水平下，依照样本容量和回归变量数目，查D-W分布表，得到检验的临界值d_L和d_U，继而作出判断^[1]。

计算可得DW= 0.3014<d_L， =0.8493，α= 0.05，n=183，k=4，查表知，ε_i存在自相关。

作变换

然后，利用变换后的数据重新估计。

模型求解

利用MATLAB工具箱计算，得到结果如下表：

回归系数

回归系数估计值

回归系数置信区间

6.0977

[5.6641,6.5313]

0.2714

[0.2438,0.2989]

0.1875

[0.1610,0.2139]

0.1271

[0.1071,0.1471]

R² = 0.9124 F=618.2362 p<0.0001 s² = 0.1071

回归方程：

把、、、还原为原始变量、、、得到结果为

结果分析

对上述模型在进行一次自相关检验，即诊断随机误差u_i是否还存在自相关。计算出DW = 1.4430，α= 0.05，n=182，k=4 ，查表知，d_u<DW<4-d_u。即随机误差u_i不存在自相关。

基于时间序列模型的回归模型与前面的多元线性回归预测值，如下图：

基于时间序列模型的回归模型与前面的多元线性回归残差图，如下图：

从两图中易知，基于时间序列，即加入自相关后的模型更贴近实际，预测结果更好，也就是说其更适合。

参考文献

[1]何晓群，刘文卿.应用回归分析.清华大学出版社，2001.

[2]姜启源，谢金星.数学模型(第四版).高等教育出版社.2012.

附录一

load data

year = 1999 : 2014;

for i = 1 : 4

for j = 1 : 15

jun_zhi(j,i) = mean(data((j-1)*12+1 : j*12 ,i+1));

biao_zhun_cha(j,i) = std(data((j-1)*12+1 : j*12 ,i+1));

zhong_wei_shu(j,i) = median(data((j-1)*12+1 : j*12 ,i+1));

end

for i = 1 : 4

jun_zhi(16,i) = mean(data(181:191,i+1));

biao_zhun_cha(16,i) = std(data(181:191,i+1));

zhong_wei_shu(16,i) = median(data(181:191,i+1));

end

%完成数据的初始化

%开始绘制原始数据折线图

plot(data(:,1),data(:,2:5),'-*')

legend('总产品','农产品','矿产品','煤油电')

title('原始数据折线图')

xlabel('时间')

ylabel('商品价格指数')

%完成原始数据折线图的绘制

%开始绘制均值和中位数的折线图

figure;

plot(year,jun_zhi,'-*')

legend('总产品','农产品','矿产品','煤油电')

title('商品价格指数均值折线图')

xlabel('时间')

ylabel('均值')

figure;

plot(year,zhong_wei_shu,'-*')

legend('总产品','农产品','矿产品','煤油电')

title('商品价格指数中位数折线图')

xlabel('时间')

ylabel('中位数')

%完成均值和中位数的折线图的绘制

%开始均值按季度变化绘图

ones(63,4);

for i=1:63

jun_zhi_season(i,:) = (data(i*3-2,2:5)...

+data(i*3-1,2:5)+data(i*3,2:5))./3;

end

x_season = 1:63;

figure;

plot(x_season,jun_zhi_season,'-*');

xlabel('季度(1-63代表顺序的63个季度)');

ylabel('季度均值');

legend('总产品','农产品','矿产品','煤油电');

title('各变量按季节变化趋势');

%完成均值按季度变化图形的绘制

% 动态频数直方图

data_dymic = sort(data(:,2:5));

temp = data_dymic(191,:);

data_dymic(3:191,:) = data_dymic(2:190,:);

data_dymic(2,:)= temp;

i=1;

figure;

while(i<=190)

subplot(2,2,1);

hist(data_dymic(1:i+1,1),5);

axis([90 112 0 70]);

title('总指数频数直方图');

subplot(2,2,2);

hist(data_dymic(1:i+1,2),5);

axis([90 125 0 70]);

title('农产品价格指数频数直方图');

subplot(2,2,3);

hist(data_dymic(1:i+1,3),5);

axis([80 125 0 70]);

title('矿产品价格指数频数直方图');

subplot(2,2,4);

hist(data_dymic(1:i+1,4),5);

axis([80 130 0 65]);

title('煤油电价格指数频数直方图');

pause(0.01);

i = i+1;

end

%结束动态图的绘制

%开始总指数其他各个指数之间关系的散点图并拟合的绘制

all = data(:,2);

agriculture = data(:,3);

mining = data(:,4);

coe = data(:,5);

figure;

A = polyfit(agriculture,all,1); %拟合

all2 = polyval(A,agriculture);

plot(agriculture,all,'b*',agriculture,all2,'r-')

title('总指数与农产品价格指数的关系图');

xlabel('农产品价格指数');

ylabel('总指数');

figure;

A = polyfit(mining,all,1);

all2 = polyval(A,mining);

plot(mining,all,'b*',mining,all2,'r-')

title('总指数与矿产品价格指数的关系图');

xlabel('矿产品价格指数');

ylabel('总指数');

figure;

A = polyfit(coe,all,1);

all2 = polyval(A,coe);

plot(coe,all,'b*',coe,all2,'r-')

title('总指数与煤油电价格指数的关系图');

xlabel('煤油电价格指数');

ylabel('总指数');

%结束总指数其他各个指数之间关系的散点图并拟合的绘制

%开始总指数与各指数均值、中位数关系的散点图并拟合

figure;

subplot(2,2,1);

B = polyfit(jun_zhi(:,2),jun_zhi(:,1),1);

all3 = polyval(B,jun_zhi(:,2));

plot(jun_zhi(:,2),jun_zhi(:,1),'b*',jun_zhi(:,2),all3,'r-')

title('总指数均值与农产品指数均值');

xlabel('农产品指数均值');

ylabel('总指数均值');

subplot(2,2,2);

B = polyfit(jun_zhi(:,3),jun_zhi(:,1),1);

all3 = polyval(B,jun_zhi(:,3));

plot(jun_zhi(:,3),jun_zhi(:,1),'b*',jun_zhi(:,3),all3,'r-')

title('总指数均值与矿产品指数均值');

xlabel('矿产品指数均值');

ylabel('总指数均值');

subplot(2,2,3);

B = polyfit(jun_zhi(:,4),jun_zhi(:,1),1);

all3 = polyval(B,jun_zhi(:,4));

plot(jun_zhi(:,4),jun_zhi(:,1),'b*',jun_zhi(:,4),all3,'r-')

title('总指数均值与煤油电指数均值');

xlabel('煤油电指数均值');

ylabel('总指数均值');

figure;

subplot(2,2,1);

B =polyfit(zhong_wei_shu(:,2),zhong_wei_shu(:,1),1);

all3 = polyval(B,zhong_wei_shu(:,2));

plot(zhong_wei_shu(:,2),zhong_wei_shu(:,1),'b*',...

zhong_wei_shu(:,2),all3,'r-')

title('总指数中位数与农产品指数中位数');

xlabel('农产品指数中位数');

ylabel('总指数中位数');

subplot(2,2,2);

B =polyfit(zhong_wei_shu(:,3),zhong_wei_shu(:,1),1);

all3 = polyval(B,zhong_wei_shu(:,3));

plot(zhong_wei_shu(:,3),zhong_wei_shu(:,1),'b*',...

zhong_wei_shu(:,3),all3,'r-')

title('总指数中位数与矿产品指数中位数');

xlabel('矿产品指数中位数');

ylabel('总指数中位数');

subplot(2,2,3);

B =polyfit(zhong_wei_shu(:,4),zhong_wei_shu(:,1),1);

all3 = polyval(B,zhong_wei_shu(:,4));

plot(zhong_wei_shu(:,4),zhong_wei_shu(:,1),'b*',...

zhong_wei_shu(:,4),all3,'r-')

title('总指数中位数与煤油电指数中位数');

xlabel('煤油电指数中位数');

ylabel('总指数中位数');

%第一题结束

附录二

load data

all = data(:,2);

agriculture = data(:,3);

X = [ ones(191,1),agriculture ];

[b,bint,r,rint,stats] = regress(all,X);

bint

stats

rcoplot(r,rint) %残差图

for i =191:-1:1 %剔除异常

if(abs(rint(i,1))+abs(rint(i,2))==abs(rint(i,1)+rint(i,2)))

all(i) = [];

agriculture(i) = [];

end

X = [ ones(178,1),agriculture ];

[b,bint,r,rint,stats] = regress(all,X);

bint

stats

附录三

load data

all = data(:,2);

agriculture = data(:,3);

mining = data(:,4);

coe = data(:,5);

X = [ ones(191,1),agriculture,mining,coe ];

[b,bint,r,rint,stats] = regress(all,X);

bint

stats

rcoplot(r,rint) %残差图

for i =191:-1:1 %剔除异常**必须从后往前剔除

if(abs(rint(i,1))+abs(rint(i,2))==abs(rint(i,1)+rint(i,2)))

all(i) = [];

agriculture(i) = [];

mining(i) = [];

coe(i) = [] ;

end

X = [ ones(183,1),agriculture,mining,coe ];

[b,bint,r,rint,stats] = regress(all,X);

bint

stats

附录四

et =all-(b(1)*ones(183,1)+b(2)*agriculture+b(3)*mining+b(4)*coe);

plot(et(1:182),et(2:183),'+');

hold on

x_time =linspace(0,0,10);

y_time = linspace(-2,2,10);

plot(x_time,y_time)

plot(y_time,x_time)

hold off

s1 = 0;

s2 = 0;

for i = 2:183

s1 = s1+et(i)*et(i-1);

s2 = s2+et(i)*et(i);

end

DW = 2*(1-s1/s2)

p_time = s1/s2

%数据变换

all_time= all(2:183)-p_time*all(1:182);

agriculture_time =agriculture(2:183)-p_time*agriculture(1:182);

mining_time=mining(2:183)-p_time*mining(1:182);

coe_time= coe(2:183)-p_time*coe(1:182);

%回归

X_time = [ones(182,1),agriculture_time,mining_time,coe_time ];

[b_time,bint_time,r_time,rint_time,stats_time]=...

regress(all_time,X_time);

b_time

bint_time

stats_time

%%%%%%%%%%%%%%%%%%%%

et2 = all_time-(b_time(1)*ones(182,1)+b_time(2)*...

agriculture_time+b_time(3)*mining_time+b_time(4)*coe_time);

s1 = 0;

s2 = 0;

for i = 2:182

s1 = s1+et2(i)*et2(i-1);

s2 = s2+et2(i)*et2(i);

end

DW = 2*(1-s1/s2)

figure;

yuce_multi = b(1)*ones(183,1)+b(2)*agriculture+b(3)*mining+b(4)*coe;

yuce_time(1) = 95.88;

yuce_time(2:183) =6.0977*ones(182,1)+0.8493*all(1:182)+0.2714...

*agriculture(2:183)-0.2305*agriculture(1:182)+0.1875*mining...

(2:183)-0.1592*mining(1:182)+0.1271*coe(2:183)-0.1079*coe(1:182);

count = 1:183;

plot(count,all,'o',count,yuce_time,'*',count,yuce_multi,'+')

legend('总指数','时间序列','多元线性回归');

figure;

et_time = all - yuce_time';

plot(count,et_time,'*',count,et,'+')

legend('时间序列','多元线性回归');

hold on

x_time =linspace(0,183,10);

y_time = linspace(0,0,10);

plot(x_time,y_time)

hold off

0 0