多元线性回归matlab实现

% 其解决问题的大致方法、步骤如下：
% 1.根据研究目的收集数据和预分析，收集一组包含因变量和自变量的数据；
% 2.根据散点图是否具有线性关系建立基本回归模型，选定因变量和自变量之间的模型，即一个数学式子，利用数据按照最小二乘准则计算模型中的系数；
% 3.利用统计分析方法对不同的模型进行比较，找出与数据拟合得最好的模型；
% 4.检验得到的模型是否适合于这组数据；
% 5.利用模型对因变量作出预测或解释。
% 每次修改输入输出数据，m、n和回归的循环
%…………………………………………
clc
%第一步：原始数据处理与线性分析
% Y:
%公路客运量(万人)
glkyl=[12815 15543 19326 22864 26150 28468 30882 39375 45759 49589 52560 48726 51083 56495 62767 83606 92090 101370 107317 108654 111847 112872 116997 126007 128980];
%公路货运量(万吨)
glhyl=[2690 2998 3012 3042 3616 3728 3988 9397 17680 19426 24128 24354 22879 24162 28957 36439 40593 45052 47400 45224 45338 45815 47151 55705 63532];
% X:
%市区人数(万人)
sqrs=[47.8 52.2 59 63.1 68.5 70 72 79.2 84.7 88.6 91 93 97.5 103.7 110 123 129.6 132 137.6 141 145 155.5 157 163.1 165.9];
%市区机动车数(万辆)
sqjdcs=[1.2 1.5 1.7 1.8 2.1 2.7 2.9	3.2	3.4	3.7	4.3	4.4	4.5	4.7	5 5.2 5.4 5.7 5.9 6.2 6.3 6.7 7.2 7.5 7.9];
%市区公路面积(万平方公里)
sqglmj=[0.2 0.25 0.25 0.3 0.45 0.5 0.5 0.7 0.7 0.75 0.8 0.8 0.85 1.1 1.25 1.3 1.3 1.5 1.55 1.75 1.8 1.8 2.05 2.1 2.3];

% 作出因变量Y与各自变量的样本散点图
% 作散点图的目的主要是观察因变量Y与各自变量间是否有比较好的线性关系，以便选择恰当的数学模型形式
subplot(2,3,1),plot(sqrs,glkyl,'g*'),
subplot(2,3,2),plot(sqjdcs,glkyl,'k+'),
subplot(2,3,3),plot(sqglmj,glkyl,'ro'),
subplot(2,3,4),plot(sqrs,glhyl,'g*'),
subplot(2,3,5),plot(sqjdcs,glhyl,'k+'),
subplot(2,3,6),plot(sqglmj,glhyl,'ro')
%若这些点大致分布在一条直线旁边，有比较好的线性关系，可以采用线性回归
%%

bb=[sqrs;sqjdcs;sqglmj]';  %输入数据矩阵
cc=[glkyl;glhyl]';         %输出数据矩阵Y
m=3                        %自变量个数
n=25                       %样本数量n
V=[];                      %将来的系数矩阵β
sh=[];
X=[ones(n,1) , bb];       %自变量矩阵（注意：第一列全为1）
% ones(n,1)表示有常数项，如果无常数项可以去掉
alpha=0.0005
%%
%第二步：回归法求回归系数[b,bint,r,rint,stats] = regress(y,X,alpha)
% b：回归系数，bint：回归系数的区间估计，r：残差，rint：置信区间
% stats：用于检验回归模型的统计量。有4个数值：判定（决定）系数r^2，F统计量观测值，检验的p的值，误差方差的估计
% alpha：是显著性水平p（缺省时默认0.05）
for n=1:2    % 因为有两个y，所以回归两次
    [b,bint,r,rint,ss]=regress(cc(:,n),X,alpha); %regress函数
    V(:,n)=b;              %系数矩阵
    sh(:,n)=r;             %残差矩阵
    ss(n,:)=ss             %检验的相关值
end

%% 
%模型检验:利用检验统计量R，F，p　的值判断该模型是否可用。
% 模型合格的要求：
% 相关系数R绝对值在0.8~1范围内。
if ss(1,1)<0.81 | ss(2,1)<0.81
    disp("NO:相关系数<0.81")
end
% F>F(m,n-m-1),后者用finv(p,m,n)计算，其中m为自变量的数量，n为样本数量。
F=finv(alpha,m,n)
if ss(1,2)<F | ss(2,2)<F
    disp("NO:F<F(m,n-m-1)")
end
% p<alpha
if ss(1,3)>F | ss(2,3)>F
    disp("NO:p>alpha")
end
% 绘制残差图：以观测值序号为横坐标，残差为纵坐标所得到的散点图称为时序残差图
rcoplot(r,rint)
%%
%预测结果带回验算  并绘图与原始数据对比
nf=1980:2004;
RLW=X*V;
figure;plot(nf,cc,'b:+');
hold on
plot(nf,RLW,'r-.');
grid on

%利用线性拟合出的数据来预测05、06的公路客运量和货运量
nh=[170.67	175.7
7.927	8.1949
2.1955	2.2818]';
nn=[zeros(2,1),nh];
ycc=nn*V;  %预测结果