基于matlab的多元线性回归分析

二、多元线性回归原理

2.1、数学模型

在社会生活及生产实践中会经常遇到一种问题，即我们非常关注一个量的变化，而这个量受到另一个或是多个因素的影响，我们想要了解这些因素是如何影响我们最为关注的这个量的以及这些因素对我们最为关注的这个量的影响权重分别有多大，知道了这些，我们就可以对该量变化所反映的相关问题做出分析和评价，并对其未来发展趋势进行预测和控制，这里就要用到数理统计中一个非常重要而普遍的分析方法，即回归分析法。

如果一个因变量y与k个自变量x1,x2,…,xkx_{1}, x_{2}, \dots, x_{k}x1,x2,…,xk存在线性相关关系，那么就可以用多元线性回归模型
y=a0+a1x1+a2x2+…+akxk,公式1y=a_{0}+a_{1} x_{1}+a_{2} x_{2}+\ldots+a_{k} x_{k},\text{公式1} y=a0+a1x1+a2x2+…+akxk,公式1
对其进行描述，其中未知常量a0,a1,…,aka_{0}, a_{1}, \dots, a_{k}a0,a1,…,ak称为回归模型系数，若n次抽样，第iii次抽样数据为(yi,x1i,x2i,⋯,xki)\left(y_{i}, x_{1 i}, x_{2 i}, \cdots, x_{k i}\right)(yi,x1i,x2i,⋯,xki)那么就有
{y1=a0+a1x11+a2x21+…+akxk1+ε1y2=a0+a1x12+a2x22+…+akxk2+ε2⋮yn=a0+a1x1i+a2x2i+…+aixki+εn公式2\left\{\begin{array}{l}y_{1}=a_{0}+a_{1} x_{11}+a_{2} x_{21}+\ldots+a_{k} x_{k1}+\varepsilon_{1} \\ y_{2}=a_{0}+a_{1} x_{12}+a_{2} x_{22}+\ldots+a_{k} x_{k2}+\varepsilon_{2} \\ \vdots \\ y_{n}=a_{0}+a_{1} x_{1 i}+a_{2} x_{2 i}+\ldots+a_{i} x_{ki}+\varepsilon_{n}\end{array}\right.\text{公式2} ⎩⎪⎪⎪⎨⎪⎪⎪⎧y1=a0+a1x11+a2x21+…+akxk1+ε1y2=a0+a1x12+a2x22+…+akxk2+ε2⋮yn=a0+a1x1i+a2x2i+…+aixki+εn公式2
其中ε0,ε1,…,εn\varepsilon_{0}, \varepsilon_{1}, \ldots, \varepsilon_{n}ε0,ε1,…,εn为随机误差项，回归分析的主要任务就是以误差ε0,ε1,…,εn\varepsilon_{0}, \varepsilon_{1}, \ldots, \varepsilon_{n}ε0,ε1,…,εn的平方和最小为原则，求多元回归模型的回归系数a0,a1,…,aka_{0}, a_{1}, \dots, a_{k}a0,a1,…,ak。

求解这个方程是要以S=∑i=1i=nεi2=∑i=1i=n(a0+a1x1i+⋯+akxki−yi)2S=\sum_{i=1}^{i=n} \varepsilon_{i}^{2}=\sum_{i=1}^{i=n}\left(a_{0}+a_{1} x_{1 i}+\dots+a_{k} x_{k i}-y_{i}\right)^{2}S=∑i=1i=nεi2=∑i=1i=n(a0+a1x1i+⋯+akxki−yi)2为最小原则，求a0,a1,…,aka_{0}, a_{1}, \dots, a_{k}a0,a1,…,ak要使得S最小，应该满足∂S∂aj=0,j=0,1,…,k\frac{\partial S}{\partial a_{j}}=0, j=0,1,\ldots,k∂aj∂S=0,j=0,1,…,k
即{∑i=1i=n2(a0+a1x1i+a2x2i+⋯+akxki−yi)=0∑i=1i=n2(a0+a1x1i+a2x2i+⋯+akxki−yi)x1i=0⋮∑i=1i=n2(a0+a1x1i+a2x2i+⋯+akxki−yi)xni=0,公式3即\left\{\begin{array}{l}\sum_{i=1}^{i=n} 2\left(a_{0}+a_{1} x_{1 i}+a_{2} x_{2 i}+\dots+a_{k} x_{k i}-y_{i}\right)=0 \\ \sum_{i=1}^{i=n} 2\left(a_{0}+a_{1} x_{1 i}+a_{2} x_{2 i}+\dots+a_{k} x_{k i}-y_{i}\right) x_{1 i}=0 \\ \vdots \\ \sum_{i=1}^{i=n} 2\left(a_{0}+a_{1} x_{1 i}+a_{2} x_{2 i}+\dots+a_{k} x_{k i}-y_{i}\right) x_{n i}=0\end{array}\right.,\text{公式3} 即⎩⎪⎪⎪⎨⎪⎪⎪⎧∑i=1i=n2(a0+a1x1i+a2x2i+⋯+akxki−yi)=0∑i=1i=n2(a0+a1x1i+a2x2i+⋯+akxki−yi)x1i=0⋮∑i=1i=n2(a0+a1x1i+a2x2i+⋯+akxki−yi)xni=0,公式3

有{na0+∑i=1i=nx1ia1+⋯+∑i=1i=nxkiak=∑i=1i=nyi∑i=1i=nx1ia0+∑i=1i=nx1i2a1+⋯+∑i=1i=nx11xkiak=∑i=1i=nx1iyi⋮∑i=1i=nxkia0+∑i=1i=nx1ixkia1+⋯+∑i=1i=nxki2ak=∑i=1i=nxkiyi,公式4有\left\{\begin{array}{l}n a_{0}+\sum_{i=1}^{i=n} x_{1 i} a_{1}+\dots+\sum_{i=1}^{i=n} x_{k i} a_{k}=\sum_{i=1}^{i=n} y_{i} \\ \sum_{i=1}^{i=n} x_{1 i} a_{0}+\sum_{i=1}^{i=n} x_{1 i}^{2} a_{1}+\dots+\sum_{i=1}^{i=n} x_{11} x_{k i} a_{k}=\sum_{i=1}^{i=n} x_{1 i} y_{i} \\ \vdots \\\sum_{i=1}^{i=n} x_{k i} a_{0}+\sum_{i=1}^{i=n} x_{1 i} x_{k i} a_{1}+\dots+\sum_{i=1}^{i=n} x_{k i}^{2} a_{k}=\sum_{i=1}^{i=n} x_{k i} y_{i}\end{array}\right.,\text{公式4} 有⎩⎪⎪⎪⎨⎪⎪⎪⎧na0+∑i=1i=nx1ia1+⋯+∑i=1i=nxkiak=∑i=1i=nyi∑i=1i=nx1ia0+∑i=1i=nx1i2a1+⋯+∑i=1i=nx11xkiak=∑i=1i=nx1iyi⋮∑i=1i=nxkia0+∑i=1i=nx1ixkia1+⋯+∑i=1i=nxki2ak=∑i=1i=nxkiyi,公式4

上式可以写成形式：Y=XAY=X AY=XA

其中：
X=[n∑i=1i=nx1i…∑i=1i=nxki∑i=1i=nx1i∑i=1i=nx1i2…∑i=1nx11xki⋮⋮…⋮∑i=1i=nxki∑i=1i=nx11xki…∑i=1i=nxki2],Y=[∑i=1i=nyi∑i=1i=nx1iyi⋮∑i=1i=nx2iyi],A=[a0a1⋮ak]X=\left[\begin{array}{ccc}n & \sum_{i=1}^{i=n} x_{1 i} & \ldots & \sum_{i=1}^{i=n} x_{k i} \\ \sum_{i=1}^{i=n} x_{1 i} & \sum_{i=1}^{i=n} x_{1 i}^{2} & \ldots& \sum_{i=1}^{n} x_{11} x_{k i} \\ \vdots& \vdots & \ldots & \vdots \\ \sum_{i=1}^{i=n} x_{k i} & \sum_{i=1}^{i=n} x_{11} x_{k i} & \ldots & \sum_{i=1}^{i=n} x_{k i}^{2}\end{array}\right], Y=\left[\begin{array}{c}\sum_{i=1}^{i=n} y_{i} \\ \sum_{i=1}^{i=n} x_{1 i} y_{i} \\ \vdots \\\sum_{i=1}^{i=n} x_{2 i} y_{i}\end{array}\right], A=\left[\begin{array}{c}a_{0} \\ a_{1} \\ \vdots\\ a_{k}\end{array}\right] X=⎣⎢⎢⎢⎡n∑i=1i=nx1i⋮∑i=1i=nxki∑i=1i=nx1i∑i=1i=nx1i2⋮∑i=1i=nx11xki…………∑i=1i=nxki∑i=1nx11xki⋮∑i=1i=nxki2⎦⎥⎥⎥⎤,Y=⎣⎢⎢⎢⎡∑i=1i=nyi∑i=1i=nx1iyi⋮∑i=1i=nx2iyi⎦⎥⎥⎥⎤,A=⎣⎢⎢⎢⎡a0a1⋮ak⎦⎥⎥⎥⎤

公式2也可以直接写成矩阵表达式：Y=XA+EY=X A+EY=XA+E

其中：
Y=[y1y2⋮yn],X=[1x1,1…xk,11x1,2…xk,2⋮⋮…⋮1x1,k…x2,k],A=[a0a1⋮ak],E=[ε1ε2⋮εn]Y=\left[\begin{array}{l}y_{1} \\ y_{2} \\ \vdots \\ y_{n}\end{array}\right], X=\left[\begin{array}{ccc}1 & x_{1,1} &\ldots & x_{k,1} \\ 1 & x_{1,2} &\ldots& x_{k,2} \\ \vdots & \vdots &\ldots &\vdots & \\ 1 & x_{1,k}&\ldots & x_{2,k}\end{array}\right], A=\left[\begin{array}{l}a_{0} \\ a_{1} \\ \vdots \\ a_{k}\end{array}\right], E=\left[\begin{array}{c}\varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{n}\end{array}\right] Y=⎣⎢⎢⎢⎡y1y2⋮yn⎦⎥⎥⎥⎤,X=⎣⎢⎢⎢⎡11⋮1x1,1x1,2⋮x1,k…………xk,1xk,2⋮x2,k⎦⎥⎥⎥⎤,A=⎣⎢⎢⎢⎡a0a1⋮ak⎦⎥⎥⎥⎤,E=⎣⎢⎢⎢⎡ε1ε2⋮εn⎦⎥⎥⎥⎤
那么我们的任务就是求解AAA

2.2、案例分析

分析总能耗与其他变量的关系，首先分析其相关性，做相似矩阵分析，如下：

第一步：得出总能耗与那些变量之间的关系，求出相关系数

[Data,str]=xlsread('C:\Users\86188\Desktop\仿真数据\北京参数一百组新.xlsx','sheet1','A1:R101',0.4);% 得到表格中所有数据
Resemblance=corrcoef(Data);                                 % 得到系数相关矩阵
[Row,Col]=size(Resemblance);                                % 得到Resemblance矩阵的行和列
site=[];count=0;                                            %
for row=1:(Row-1)                                           % 得到总能源与哪些因素有关if(abs(Resemblance(row,Col))>0.4)disp(['总能源与第' num2str(row) '列' str(row) '相关度较高，相关系数为 ' num2str(Resemblance(Col,row)) ]);count=count+1;                                      % 保存相关变量的数量site(count)=row;                                    % 保存相关变量的列地址end
end

第二步：得到回归系数和置信区间

[Row,Col]=size(Data);                                       % 得到数据矩阵的行和列的大小
ConVariable=zeros(Row,count);                               % 创建一个行相等列指定的矩阵
for row=1:countConVariable(:,row)=Data(:,site(row));
end
TotalEnergy=Data(:,Col);                                   % 得到总能源原始数据
IndeVariable=[ones(Row,1),ConVariable];                    % 创建相关变量数组
[b,bint,r,rint,stats]=regress(TotalEnergy,IndeVariable);   % 求回归系数的点估计和区间估计、并检验回归模型

第三步：做残差分析，得出那些点偏差太大

得出四个点偏差太大，可以考虑去除这四个点

subplot(211);                                                           % 画残差图
rcoplot(r,rint);                                                        % 画残差图

第四步：得出预测模型，用预测模型去和真实值去对比

NewTotalEnergy=zeros(Row,1);                                            % 预测数据矩阵
for row = 1:Row                                                         % 预测赋值NewTotalEnergy(row)=b(1);for PaRow=2:(count+1)NewTotalEnergy(row)=NewTotalEnergy(row)+b(PaRow)*ConVariable(row,PaRow-1);end
end
Loss = 0;                                   % 根据最大斜率和最小效率建立指定长度的损失数据
for i=1:RowLoss = Loss+(TotalEnergy(i)-NewTotalEnergy(i))^2/(2*Row);
end
subplot(212);                                                         % 画预测和真实图
number=[1:1:Row];
plot(number,TotalEnergy','r',number,NewTotalEnergy','b');
xlabel('数量序列');ylabel('总能耗');title('总能耗与相关参数散点图');legend('真实值','预测值')grid on;

第五步：得出模型函数关系
总能耗=33374+−53041×体型系数+32×面积−101×人口密度+431×内扰电耗总能耗=33374 +-53041\times 体型系数 +32\times 面积 -101\times 人口密度 +431\times 内扰电耗总能耗=33374+−53041×体型系数+32×面积−101×人口密度+431×内扰电耗

参考文献

数学建模与数学试验
多元线性回归MATLAB实现
一元线性回归模型及其假设条件
[MATLAB]逐步回归详解(stepwise使用指南)
基于Matlab的数据多元回归分析的研究