function [ b, bint, r, rint, stats ] = Multiple_linear_regression( X,y,alpha )
%Multiple_linear_regression 多元线性回归
%对于输入数据(必须是去除量纲之后的数据)进行多元线性回归处理
% 输入变量
% X：由自变量组成的矩阵(第一列必须全为1)
% y：因变量（列向量）
% alpha：显著性水平（缺省时默认为0.05）
% 若X中两个原本全为变量，需第一列加上1，作为最后常数项的系数预测
% 输出变量
% b为线性回归模型Y = X*B回归系数（即斜率参数）第一个是常数项系数
% bint为b的置信区间(按行显示)
% r为残差（实际值与拟合值的差）
% rint为残差的置信区间
% stats包含如下四个统计量
% R^2统计量（越接近1，变量的线性相关性越强，说明模型有效）
% F统计量
% F(1,n-2)分布大于F的概率p值（显著性水平，判断原始假设是否正确的重要证据）
% 剩余方差s^2（主要用来比较模型是否有改进，越小则模型精度越高）
X = [ones(size(y)),X];
[b, bint, r, rint, stats] = regress(y,X,alpha);
if stats(3)<0.01
disp('F(1,n-2)分布大于F的概率p值为：');
disp( stats(3) );
disp('p值很小（P<0.001）,说明拟合模型有效');
else
disp('p值较大（P>0.001）,说明拟合模型效果不好');
end
%% 画出拟合效果图
figure
y_fitting = X*b;
t = 1:length(y);
plot(t,y_fitting,'r-',t,y,'b-',t,abs(y_fitting-y),'k-');
legend({'红——拟合值','蓝——实际值','黑——误差值'},'fontsize',12,'Location','NorthWest');
xlabel('点的序号','fontsize',12);
ylabel('y对应的值','fontsize',12);
set(gca, 'XGrid','on'); % X轴的网格
set(gca, 'YGrid','on'); % Y轴的网格
end
