基于MATALB语言实现线性回归分析

摘要

针对各学科领域中常遇到的一元或者多元线性回归问题，在简单介绍回归分析基本理论的基础上，结合具体实例，详细介绍了基于回归算法编写MATLAB程序、利用MATLAB预定义函数以及二者相结合解决多元线性回归问题的方法。再根据已得的实验结果以及以往的经验来建立统计模型，并研究变量之间的相关关系，建立起变量之间关系的近似表达式，并由此对相应的变量进行预测和控制。

关键字一元回归分析，多元回归分析，MATLAB

一、一元线性回归原理

1.1、数学模型

一元线性回归分析是在排除其他影响因素，分析某一个因素（自变量：X）是如何影响另外一个事物（因变量：Y）的过程，所进行的分析是比较理想化的对于一元线性回归来说，可以看成Y的值是随着X的值变化，每一个实际的X都会有一个实际的Y值，我们叫Y实际，那么我们就是要求出一条直线，每一个实际的X都会有一个直线预测的Y值，我们叫做Y预测，回归线使得每个Y的实际值与预测值之差的平方和最小，即达到一元线性回归的最终结果。

一般的，一元线性回归模型可由下表示：
$Y=\beta_0 +\beta_1\times X+\epsilon$Y=β0​+β1​×X+ϵ
固定的$\beta_0$β0​、$\beta_1$β1​称为回归系数，自变量X也成为回归变量，$Y=\beta_0 +\beta_1\times X+\epsilon$Y=β0​+β1​×X+ϵ,称为Y对X的回归直线方程，且$\epsilon$ϵ 的均值$E(\epsilon)=0$E(ϵ)=0，所以模型简化为$Y=\beta_0 +\beta_1\times X$Y=β0​+β1​×X。

对于实际问题，要建立回归方程，首先确定能否建立线性回归模型，其次确定如何对模型中未知参数$\beta_0$β0​、$\beta_1$β1​进行估计。所以首先对总体进行独立观测，在坐标系中画出Y与X的散点图，根据图判断Y与X直线关系是否符合一元线性回归模型。最后利用最小二乘法可以得到回归模型参数$\beta_0$β0​、$\beta_1$β1​的最小二乘估计，估计公式为：
$\begin{cases} \beta_{0}=\bar{y}-\bar{x} \beta_{1} \\ \beta_{1}=\frac{L_{x y}}{L_{x x}} \end{cases}${β0​=yˉ​−xˉβ1​β1​=Lxx​Lxy​​​
其中：$\bar{x}=\frac{1}{n} \sum_{i=1}^{n} X_{i}, \bar{y}=\frac{1}{n} \sum_{i=1}^{n} y_{i}$xˉ=n1​∑i=1n​Xi​,yˉ​=n1​∑i=1n​yi​，$L_{x x}=\sum_{i=1}^{n}\left(X_{i}-\bar{X}\right)^{2}, L_{x y}=\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)$Lxx​=∑i=1n​(Xi​−Xˉ)2,Lxy​=∑i=1n​(xi​−xˉ)(yi​−yˉ​)

于是我们建立公式模型：$y=\beta_{0}+\beta_{1} x$y=β0​+β1​x

接下来我们可以通过对一元线性回归分析总结为三点

1. 利用样本观测值对回归系数$\beta_0$β0​、$\beta_1$β1​做点估计
2. 对方程的线性关系做显著校验
3. 在X=$X_0$X0​处对Y做预测等

1.2、案例分析

分析总能耗与面积的关系，首先导入散点图，如下：

第一步:画出面积和总能耗之间的散点图，如下：

[Data,str]=xlsread('C:\Users\86188\Desktop\仿真数据\北京参数一百组新.xlsx','sheet1','A1:R101');%得到表格中所有数据
[Row,Col]=size(Data);                    % 得到数据的行和列宽
NewData=zeros(Row,2);                    % 建立两列行相等的数组
NewData(:,1)=Data(:,factor);             % 保存因素数据
NewData(:,2)=Data(:,Col);                % 保存自变量数据

subplot(311);                            % 分图1,肉眼看线性关系
plot(NewData(:,1)',NewData(:,2)','k+');  % 得出面积与总能耗之间的散点图，通过散点图可以明确发现存在线性关系
axis([0,inf,0,inf]);xlabel('面积');ylabel('总能耗');title('总能耗与面积散点图');legend('真实值');

可以观测处面积与总能耗之间存在线性关系，但是实际上不用标准的一元线性模型，而使用$Y=\beta_1\times X$Y=β1​×X,观测Y与X的具体关系，不加截距的方式。

第二步：求出面积与总能耗之间最大与最小斜率。

slop=zeros(Row,1);                        % 建立一个斜率数组
for i=1:Row
    slop(i)=NewData(i,2)/NewData(i,1);    % 得到每个数据对应的斜率
end
[MaxNumber,MaxIndex] = max(slop);         % 得到最大斜率和下标
[MinNumber,MinIndex] = min(slop);         % 得到最小效率和下标

第三步：求出每个斜率对应的误差值

最终得到在 K=63.88的时候可以得到最小误差RSS=5296700。

SlopNumber = (MaxNumber-MinNumber)/0.01;  % 得到以0.01为精度的斜率个数
Loss=zeros(SlopNumber,1);                 % 根据最大斜率和最小效率建立指定长度的损失数据       
count = 1;
for PreSlop = MinNumber:0.01:(MaxNumber-0.01)    % 斜率以0.01的增量去计算损失函数
    for i=1:Row
        Loss(count) = Loss(count)+(NewData(i,2)-PreSlop*NewData(i,1))^2/(2*100);
    end
    count=count+1;
end
Number=1:SlopNumber;                     % 得到损失函数列向量数量
subplot(312);                            % 分图2：看损失函数的最小值
plot(Number,Loss');ylabel('损失值');title('损失函数图');legend('损失函数');

第四步：得到最佳斜率并用模型进行预测总能耗

[MinLoss,MinLossindex] = min(Loss)      % 得到损失函数最小值及对应的下标
PreSlop = zeros(2,Row);                  % 建立预测曲线存取矩阵
for i=1:Row                              % 预测矩阵基于给定的数据，计算y=k*x的关系
    PreSlop(2,i)= (MinLossindex*0.01+MinNumber)*NewData(i,1);
    PreSlop(1,i)= NewData(i,1);
end
subplot(313);                            % 分图3：看实际数据与线性回归的关系
plot(NewData(:,1)',NewData(:,2)','k+',PreSlop(1,:)',PreSlop(2,:),'r');% 得出面积与总能耗之间的散点图，通过散点图可以明确发现存在线性关系
axis([0,inf,0,inf]);xlabel('面积');ylabel('总能耗');title('总能耗与面积的线性关系');
legend('真实值','预测值');grid on;hold on;

第五步：得到模型函数关系
$总能耗=63.88 \times 面积$总能耗=63.88×面积

二、多元线性回归原理

2.1、数学模型

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

如果一个因变量y与k个自变量$x_{1}, x_{2}, \dots, x_{k}$x1​,x2​,…,xk​存在线性相关关系，那么就可以用多元线性回归模型
$y=a_{0}+a_{1} x_{1}+a_{2} x_{2}+\ldots+a_{k} x_{k},\text{公式1}$y=a0​+a1​x1​+a2​x2​+…+ak​xk​,公式1
对其进行描述，其中未知常量$a_{0}, a_{1}, \dots, a_{k}$a0​,a1​,…,ak​称为回归模型系数，若n次抽样，第$i$i次抽样数据为$\left(y_{i}, x_{1 i}, x_{2 i}, \cdots, x_{k i}\right)$(yi​,x1i​,x2i​,⋯,xki​)那么就有
$\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​+a1​x11​+a2​x21​+…+ak​xk1​+ε1​y2​=a0​+a1​x12​+a2​x22​+…+ak​xk2​+ε2​⋮yn​=a0​+a1​x1i​+a2​x2i​+…+ai​xki​+εn​​公式2
其中$\varepsilon_{0}, \varepsilon_{1}, \ldots, \varepsilon_{n}$ε0​,ε1​,…,εn​为随机误差项，回归分析的主要任务就是以误差$\varepsilon_{0}, \varepsilon_{1}, \ldots, \varepsilon_{n}$ε0​,ε1​,…,εn​的平方和最小为原则，求多元回归模型的回归系数$a_{0}, a_{1}, \dots, a_{k}$a0​,a1​,…,ak​。

求解这个方程是要以$S=\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​+a1​x1i​+⋯+ak​xki​−yi​)2为最小原则，求$a_{0}, a_{1}, \dots, a_{k}$a0​,a1​,…,ak​要使得S最小，应该满足$\frac{\partial S}{\partial a_{j}}=0, j=0,1,\ldots,k$∂aj​∂S​=0,j=0,1,…,k
$即\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=n​2(a0​+a1​x1i​+a2​x2i​+⋯+ak​xki​−yi​)=0∑i=1i=n​2(a0​+a1​x1i​+a2​x2i​+⋯+ak​xki​−yi​)x1i​=0⋮∑i=1i=n​2(a0​+a1​x1i​+a2​x2i​+⋯+ak​xki​−yi​)xni​=0​,公式3

$有\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=n​x1i​a1​+⋯+∑i=1i=n​xki​ak​=∑i=1i=n​yi​∑i=1i=n​x1i​a0​+∑i=1i=n​x1i2​a1​+⋯+∑i=1i=n​x11​xki​ak​=∑i=1i=n​x1i​yi​⋮∑i=1i=n​xki​a0​+∑i=1i=n​x1i​xki​a1​+⋯+∑i=1i=n​xki2​ak​=∑i=1i=n​xki​yi​​,公式4

上式可以写成形式：$Y=X A$Y=XA

其中：
$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=n​x1i​⋮∑i=1i=n​xki​​∑i=1i=n​x1i​∑i=1i=n​x1i2​⋮∑i=1i=n​x11​xki​​…………​∑i=1i=n​xki​∑i=1n​x11​xki​⋮∑i=1i=n​xki2​​⎦⎥⎥⎥⎤​,Y=⎣⎢⎢⎢⎡​∑i=1i=n​yi​∑i=1i=n​x1i​yi​⋮∑i=1i=n​x2i​yi​​⎦⎥⎥⎥⎤​,A=⎣⎢⎢⎢⎡​a0​a1​⋮ak​​⎦⎥⎥⎥⎤​

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

其中：
$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=⎣⎢⎢⎢⎡​y1​y2​⋮yn​​⎦⎥⎥⎥⎤​,X=⎣⎢⎢⎢⎡​11⋮1​x1,1​x1,2​⋮x1,k​​…………​xk,1​xk,2​⋮x2,k​​​⎦⎥⎥⎥⎤​,A=⎣⎢⎢⎢⎡​a0​a1​⋮ak​​⎦⎥⎥⎥⎤​,E=⎣⎢⎢⎢⎡​ε1​ε2​⋮εn​​⎦⎥⎥⎥⎤​
那么我们的任务就是求解$A$A

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:count
    ConVariable(:,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:Row
    Loss = 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\times 体型系数 +32\times 面积 -101\times 人口密度 +431\times 内扰电耗$总能耗=33374+−53041×体型系数+32×面积−101×人口密度+431×内扰电耗

三、结论

对比一元线性回归和多元线性回归可以发现线性模型总体误差：

模型	RSS
一元	5296700
多元	1751100

得出多元线性回归在因变量受到多个因素可的影响下，可以大幅度减小真实值和预测值之间的误差。

参考文献：

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