线性回归是利用数理统计中回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，在线性回归分析中，只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归分析称为一元线性回归分析。如果回归分析中包括两个或两个以上的自变量，且因变量和自变量之间是线性关系，则称为多元线性回归分析。

一元线性回归

定义

一般地，称由y=β_0+β_1x+ε确定的模型称为一元线性回归模型，记作：

\left\{
\begin{aligned}
&y=β_0+β_1x+ε\\
&E(ε)=0,D(ε)=δ^2
\end{aligned}\right.

固定的未知参数β_0、β_1称为回归系数，自变量x也称为回归变量。y=β_0+β_1x称为y对x的回归直线方程。

主要任务

• 用试验值（样本值）对β_0、β_1和δ作点估计；
• 对回归系数β_0、β_1作假设检验；
• 在x=x_0处对y作预测，对y作区间估计。

普通最小二乘法（Ordinary Least Square, OLS）

给定一组样本观测值（Xi, Yi），i=1,2,…n，假如模型参数估计量已经求得，并且是最合理的参数估计量，那么样本回归函数应该能够最好地拟合样本数据，即样本回归线上的点与真实观测点的“总体误差”应该
尽可能地小。
普通最小二乘法给出的判断标准是：二者之差的平方和最小，即：

Q=\sum^n_{i=1}(Y_i-\hat{Y_i})^2=\sum^n_{i=1}(Y_i-(\hat{β_0}+\hat{β_1}X_i))^2→最小

换句话说，也就是要使得所有样本点到样本回归线的数值距离平方和最小。

由于Q=\sum^n_{i=1}(Y_i-\hat{Y_i})^2=\sum^n_{i=1}(Y_i-(\hat{β_0}+\hat{β_1}X_i))^2是\hat{β_0}、\hat{β_1}的二次函数，并且非负。所以其极小值总是存在的。根据极值存在的条件，当Q对β_0、β_1的一阶偏导数为0时，Q达到最小，即：

\left\{
\begin{aligned}
&\frac{∂Q}{∂\hat{β_0}}=0\\
&\frac{∂Q}{∂\hat{β_1}}=0
\end{aligned}\right.\ →\ 
\left\{
\begin{aligned}
&\sum(Y_i-\hat{β_0}-\hat{β_1}X_i)=0\\
&\sum(Y_i-\hat{β_0}-\hat{β_1}X_i)X_i=0
\end{aligned}\right.\ →\ 
\left\{
\begin{aligned}
&\sum Y_i=n\hat{β_0}+\hat{β_1}\sum X_i\\
&\sum Y_iX_i=\hat{β_0}\sum X_i+\hat{β_1}\sum X_i^2
\end{aligned}\right.

解得：

\left\{
\begin{aligned}
&\hat{β_1}=\frac{n\sum Y_iX_i-\sum Y_i \sum X_i}{n\sum X_i^2-(\sum X_i)^2}\\
&\hat{β_0}=\overline Y-\hat{β_1}\overline X
\end{aligned}\right.

由于β_0、β_1的估计结果是从最小二乘原理得到的，故称为最小二乘估量(least-squares estimators)。

普通最小二乘参数估计量的离差形式（deviation form）

记：

\left\{
\begin{aligned}
&\overline X=\frac{1}{n}\sum X_i\\
&\overline Y=\frac{1}{n}\sum Y_i\\
&x_i=X_i-\overline X\\
&y_i=Y_i-\overline Y
\end{aligned}\right.

则参数估计量可以写成：

\left\{
\begin{aligned}
&\hat{β_1}=\frac{\sum x_iy_i}{\sum x_i^2}\\
&\hat{β_0}=\overline Y-\hat{β_1}\overline X
\end{aligned}\right.

注：在计量经济学中，往往以大写字母表示原始数据（观测值），而以小写字母表示对均值的离差（deviation)。

随机误差项方差的估计量

记e_i+Y_i=\hat{Y_i}为第i个样本观察点的残差，即被解释变量的估计值与观测值之差，则随机误差项方差的估计量为：

\hat{δ}^2_μ=\frac{\sum e^2_i}{n-2}

回归方程的显著性检验

对回归方程y=β_0+β_1x的显著性检验，归结为对假设：

H_0:β_1=0;H_1:β_1\not=0

进行检验。假设H_0:β_1=0被拒绝，则回归显著，认为y与x存在线性关系，所求的线性回归方程也有意义；否则认为回归不显著，y与x的关系不能用一元线性回归模型来描述，所求的线性回归方程也无意义。

F检验法

当H_0成立时，

F=\frac{U}{Q_e/(n-2)} \sim F(1,n-2)

其中回归平方和U：

U=\sum^n_{i=1}(\hat{y_i}-\overline y)^2

若F>F_{1-α}(1,n-2)，拒绝H_0，否则就接受H_0。

T检验法

当H_0成立时，

T=\frac{\sqrt{L_{xx}}\hat{β_1}}{\hat{δ_e}}\sim t(n-2),其中L_{xx}=\sum^n_{i=1}(x_i-\overline x)^2=\sum^n_{i=1}x^2_i-n\overline x^2

若|T|>t_{1-\frac{α}{2}}(n-2)，拒绝H_0，否则就接受H_0。

r检验法

记

r=\frac{\sum^n_{i=1}(x_i-\overline x)(y_i-\overline y)}{\sqrt{\sum^n_{i=1}(x_i-\overline x)^2\sum^n_{i=1}(y_i-\overline y)^2}}

若|r|>r_{1-α}，拒绝H_0，否则就接受H_0。其中：

r_{1-α}=\sqrt{\frac{1}{1+(n-2)/F_{1-α}(1,n-2)}}

回归系数的置信区间

β_0置信水平为1-α的置信区间为：

\left[\hat{β_0}-t_{1-\frac{α}{2}}(n-2)\hat{δ_e}\sqrt{\frac{1}{n}+\frac{\overline x^2}{L_{xx}}},\hat{β_0}+t_{1-\frac{α}{2}}(n-2)\hat{δ_e}\sqrt{\frac{1}{n}+\frac{\overline x^2}{L_{xx}}}\right]

β_1置信水平为1-α的置信区间为：

\left[\hat{β_1}-\frac{t_{1-\frac{α}{2}}(n-2)\hat{δ_e}}{\sqrt{L_{xx}}},\hat{β_1}+\frac{t_{1-\frac{α}{2}}(n-2)\hat{δ_e}}{\sqrt{L_{xx}}}\right]

δ^2置信水平为1-α的置信区间为：

\left[\frac{Q_e}{x^2_{1-\frac{α}{2}}(n-2)},\frac{Q_e}{x^2_{\frac{α}{2}}(n-2)}\right]

回归方程的预测与控制

预测

用y_0的回归值\hat{y_0}=\hat{β_0}+\hat{β_1}x_0作为y_0的预测值。

y_0的置信水平为1-α的预测区间为：

[\hat{y_0}-δ(x_0),\hat{y_0}+δ(x_0)]

其中：

δ(x_0)=\hat{δ_e}t_{1-\frac{α}{2}}(n-2)\sqrt{1+\frac{1}{n}+\frac{(x_0-\overline x)^2}{L_{xx}}}

特别地，当n很大且x_0在\overline x附近取值时，y的置信水平为1-α的预测区间近似为：

[\hat y-\hat{δ_e}u_{1-\frac{α}{2}},\hat y+\hat{δ_e}u_{1-\frac{α}{2}}]

控制

要求;、y=β_0+β_1x+ε的值以1-α的概率落在指定区间(y’,y”)，只要控制x满足以下两个不等式：

\hat y-δ(x)\ge y',\hat y+δ(x)\le y''

要求y''-y'\ge 2δ(x)，若ョx’,x”满足：

\hat y-δ(x')= y',\hat y+δ(x'')= y''

则(x’,x”)就是所求的x的控制区间。

多元线性回归

定义

一般称：

\left\{
\begin{aligned}
&Y=Xβ+ε\\
&E(ε)=0,COV(ε,ε)=δ^2I_n
\end{aligned}\right.

为高斯-马尔可夫线性模型（k元线性回归模型），并且简记为(Y,Xβ,δ^2I_n)。

Y=\left[
\begin{aligned}
&y_1\\
&y_2\\
&...\\
&y_n
\end{aligned}\right],
X=\left[
\begin{aligned}
&1&\ &x_{11}&\ &...&\ &x_{1k}&\\
&1&\ &x_{21}&\ &...&\ &x_{2k}&\\
&...&\ &...&\ &...&\ &...\\
&1&\ &x_{n1}&\ &...&\ &x_{nk}&
\end{aligned}\right],
β=\left[
\begin{aligned}
&β_0\\
&β_1\\
&...\\
&β_k
\end{aligned}\right],
ε=\left[
\begin{aligned}
&ε_1\\
&ε_2\\
&...\\
&ε_n
\end{aligned}\right]

y=β_0+β_1x_1+...+β_kx_k称为回归平面方程。

主要任务

• 用试验值（样本值）对β和δ作点估计与假设检验；
• 在x_1=x_{01},...,x_k=x_{0k}处对y作预测与控制，对y作区间估计。

模型参数估计

用最小二乘法求β的估计量：作离差平方和

Q=\sum^n_{i=1}(y_i-β_0-β_1x_{i1}-...-β_kx_{ik})^2

选择合适的β值使得Q达到最小，解得估计值：

\hat β=(X^TX)^{-1}X^TY

其中(X^TX)^{-1}X^T称为X的伪逆。得到的\hat{β_i}带入回归平面方程得：

y=\hat{β_0}+\hat{β_1}x_1+...+\hat{β_k}x_k

称为经验回归平面方程，\hat{β_i}称为经验回归系数。

线性模型和回归系数的检验

假设H_0:β_0=β_1=...=β_k=0。

F检验法

当H_0成立时，

F=\frac{U/k}{Q_e/(n-k-1)}\sim F(k,n-k-1)

如果F>F_{1-α}(k,n-k-1)，则拒绝H_0，认为有显著的线性关系，否则就接受H_0，认为线性关系不显著。其中：

回归平方和U=\sum^n_{i=1}(\hat y_i-\overline y)^2,
残差平方和Q_e=\sum^n_{i=1}(y_i-\hat y_i)^2

R检验法

定义：

R=\sqrt{\frac{U}{L_{yy}}}=\sqrt{\frac{U}{U+Q_e}}

为y与x的多元相关系数或复相关系数。

由于F=\frac{n-k-1}{k}\frac{R^2}{1-R^2}，故用F和用R检验是等效的。

逐步回归分析

“最优”的回归方程就是包含所有对Y有影响的变量,而不包含对Y影响不显著的变量回归方程。选择“最优”的回归方程有以下几种方法：

• 从所有可能的因子（变量）组合的回归方程中选择最优者
• 从包含全部变量的回归方程中逐次剔除不显著因子
• 从一个变量开始，把变量逐个引入方程
• “有进有出”的逐步回归分析

以第四种方法，即逐步回归分析法在筛选变量方面较为理想。

逐步回归分析法的思想：

• 从一个自变量开始，视自变量Y对作用的显著程度，从大到小地依次逐个引入回归方程
• 当引入的自变量由于后面变量的引入而变得不显著时，要将其剔除掉
• 引入一个自变量或从回归方程中剔除一个自变量，为逐步回归的一步
• 对于每一步都要进行Y值检验，以确保每次引入新的显著性变量前回归方程中只包含对Y作用显著的变量
• 这个过程反复进行，直至既无不显著的变量从回归方程中剔除，又无显著变量可引入回归方程时为止

线性回归的python实现

python中不像R中，默认的函数可以做回归分析lm，可以做方差分析aov，python中进行统计分析需要载入外在的包，这里经常用到的是statsmodels和sklearn包，statsmodels风格还是和R语言类似，sklearn则偏向机器学习了，是机器学习的入门包。

从数学建模角度出发，本文侧重于使用Statsmodels库。

在后文的所有代码中，导入的库和初始化设置如下：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import chi2_contingency

mpl.rcParams['font.family']='SimHei'
plt.rcParams['axes.unicode_minus']=False
pd.set_option('display.unicode.ambiguous_as_wide',True)
pd.set_option('display.unicode.east_asian_width',True)

Statsmodels库

Statsmodels库是Python中一个强大的统计分析库，包含假设检验、回归分析、时间序列分析等功能，能够很好的和Numpy和Pandas等库结合起来，提高工作效率。

一元线性回归示例

使用数据：

代码及输出结果：

file='E:\DEV\python\math\data\线性回归.xlsx'
data=pd.read_excel(file,sheet_name='一元线性回归')
sns.lmplot(x=data.columns[0],y=data.columns[1],data=data,truncate=True)
plt.show()
fit=smf.ols(f'{data.columns[1]}~{data.columns[0]}',data=data).fit()
print(f'y={"%f"%fit.params[0]}+{"%f"%fit.params[1]}x')
print(fit.summary())

y=7.909167+0.217500x
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   增大容积   R-squared:                       0.677
Model:                            OLS   Adj. R-squared:                  0.652
Method:                 Least Squares   F-statistic:                     27.27
Date:                Fri, 27 Jan 2023   Prob (F-statistic):           0.000165
Time:                        18:16:57   Log-Likelihood:                -14.794
No. Observations:                  15   AIC:                             33.59
Df Residuals:                      13   BIC:                             35.00
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      7.9092      0.416     19.023      0.000       7.011       8.807
使用次数           0.2175      0.042      5.222      0.000       0.128       0.307
==============================================================================
Omnibus:                       13.833   Durbin-Watson:                   0.711
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               10.255
Skew:                          -1.589   Prob(JB):                      0.00593
Kurtosis:                       5.512   Cond. No.                         23.3
==============================================================================

参数解释

Summary中值的解释如下：

• Dep. Variable
  • Which variable is the response in the model
  • 输出Y变量的名称
• Model
  • What model you are using in the fit
  • 使用的参数确定的模型
• Method
  • How the parameters of the model were calculated
  • 确定参数的方法
• Date
  • 日期
• Time
  • 时间
• No. Observations
  • The number of observations (examples)
  • 样本数目
• DF Residuals
  • Degrees of freedom of the residuals. Number of observations – number of parameters
  • 残差的自由度（等于观测数(No.Observations)-参数数目(Df Model+1(常量参数))）
  • 残差:指实际观察值与估计值（拟合值）之间的差
• Df Model
  • Number of parameters in the model (not including the constant term if present)
  • 模型参数个数（不包含常量参数），对应于coef中的行数
• R-squared 该数据反应回归数据的拟合程度
  • The coefficient of determination. A statistical measure of how well the regression line approximates the real data points
  • 可决系数的平方，说明估计的准确性。“可决系数”是通过数据的变化来表征一个拟合的好坏。由上面的表达式可以知道“可决系数”的正常取值范围为[0,1]，越接近1，表明方程的变量对y的解释能力越强，这个模型对数据拟合的也较好
• Adj. R-squared
  • The above value adjusted based on the number of observations and the degrees-of-freedom of the residuals
  • 根据观测值的数量和残差的自由度进行调整后产生的可决系数的平方
• F-statistic
  • A measure how significant the fit is. The mean squared error of the model divided by the mean squared error of the residuals
  • 衡量拟合度有多大。模型的平均平方误差除以残差的平均平方误差
• Prob (F-statistic) 该数据为通过F检验法获得的显著性指标。
  • The probability that you would get the above statistic, given the null hypothesis that they are unrelated
  • 在无效假设下，你会得到上述统计数字的概率，即它们是不相关的。
  • 通常在实际的应用中将概率值p与0.05做比较，小于0.05则拒绝原假设，否则接受原假设。
  • 对于F检验来说，如果无法拒绝原假设，则认为模型是无效的，通常解决办法是增加数据量、改变自变量或选择其他模型；
• Log-likelihood
  • The log of the likelihood function.
  • 概率函数的对数。
• AIC
  • The Akaike Information Criterion. Adjusts the log-likelihood based on the number of observations and the complexity of the model.
  • Akaike信息准则。根据观察值的数量和模型的复杂性来调整对数可能性。
• BIC
  • The Bayesian Information Criterion. Similar to the AIC, but has a higher penalty for models with more parameters.
  • 贝叶斯信息准则。与AIC相似，但对有更多参数的模型有更高的惩罚。
• coef
  • The estimated value of the coefficient
  • 系数的估计值
• std err
  • The basic standard error of the estimate of the coefficient. More sophisticated errors are also available.
  • 系数估计的基本标准误差。也可提供更复杂的误差。
• t
  • The t-statistic value. This is a measure of how statistically significant the coefficient is.
  • t统计值。这是一个衡量系数的统计学意义的指标。
• P > |t| 该数据为通过T检验法获得的显著性指标。
  • P-value that the null-hypothesis that the coefficient = 0 is true. If it is less than the confidence level, often 0.05, it indicates that there is a statistically significant relationship between the term and the response.
  • 证明系数=0的无效假设是真的P值。如果它小于置信度，通常是0.05，则表明该术语和反应之间存在统计上的显著关系。
  • 对于t检验，如果无法拒绝原假设，则认为对应的自变量与因变量之间不存在线性关系，通常的解决办法是剔除该变量或修正该变量（如选择对于的数学转换函数对其修正处理）
• [95.0% Conf. Interval]
  • The lower and upper values of the 95% confidence interval
  • 95%置信区间的下限和上限值
• Skewness
  • A measure of the symmetry of the data about the mean. Normally-distributed errors should be symmetrically distributed about the mean (equal amounts above and below the line).
  • 衡量数据对平均值的对称性。正常分布的误差应该是围绕平均值对称分布的（线上和线下的数量相等）。
• Kurtosis
  • A measure of the shape of the distribution. Compares the amount of data close to the mean with those far away from the mean (in the tails).
  • 对分布形状的衡量。将接近平均值的数据量与远离平均值的数据量（在尾部）进行比较。
• Omnibus
  • D’Angostino’s test. It provides a combined statistical test for the presence of skewness and kurtosis.
  • D’Angostino’s测试。它为偏度和峰度的存在提供了一个综合的统计测试。
• Prob(Omnibus)
  • The above statistic turned into a probability
  • 将Omnibus统计数字变成了一个概率
• Jarque-Bera
  • A different test of the skewness and kurtosis
  • 对偏度和峰度的另一种测试
• Prob (JB)
  • The above statistic turned into a probability
  • 将Jarque-Bera统计数字变成了一个概率
• Durbin-Watson
  • A test for the presence of autocorrelation (that the errors are not independent.) Often important in time-series analysis
  • 检验是否存在自相关（即误差不独立），这在时间序列分析中往往很重要
• Cond. No
  • A test for multicollinearity (if in a fit with multiple parameters, the parameters are related with each other).
  • 多重共线性的测试（如果在有多个参数的拟合中，参数之间有关联）。

使用Statsmodels库构造的线性回归类

# 线性回归
class LinearRegression:     
    def get_independent_var(self):
        if type(self.independent_var)!=type(None):
            return self.independent_var
        columns=self.vars.columns.tolist()
        if self.dependent_var in columns:
            columns.remove(self.dependent_var)
        else:
            raise ValueError('因变量不在引入的数据列表中，请检查！')
        if type(self.remove_var)!=type(None):
            for elem in self.remove_var:
                columns.remove(elem)
        return columns
    def independent_var_toString(self):
        string=''
        for i_var in self.independent_var:
            string+=f'{i_var}+'
        return string.strip('+')
    def __init__(self,data:pd.DataFrame,vars:pd.DataFrame,dependent_var:str,independent_var:list=None,remove_var:list=None):
        """
        介绍
        ----------
        线性回归是利用数理统计中回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，在线性回归分析中，\n
        只包括一个自变量和一个因变量，且二者的关系可用一条直线近似表示，这种回归分析称为一元线性回归分析。如果回归分析中包\n
        括两个或两个以上的自变量，且因变量和自变量之间是线性关系，则称为多元线性回归分析。

        参数
        ----------
        data : pd.DataFrame
            完整的原始数据表格
        vars : pd.DataFrame
            仅包含数字变量的数据表格
        dependent_var : str
            因变量
        independent_var : list, optional
            自变量，留空时默认取var中包含的所有除因变量外的其他变量
        remove_var : list, optional
            从var中去除的自变量，仅在independent_var留空时此参数有效，否则自变量以前者为准
            
        属性
        ----------
        data : pd.DataFrame
            完整的原始数据表格，由同名参数传入
        vars : pd.DataFrame
            仅包含数字变量的数据表格，由同名参数传入
        dependent_var : str
            因变量，由同名参数传入
        independent_var : list, optional
            自变量，由同名参数传入
        remove_var : list, optional
            从var中去除的自变量，由同名参数传入
        
        方法
        ----------
        summary :
            展示回归分析结果表
        show_lmplot :
            展示因变量与某个自变量之间的拟合效果图
        forecast :
            根据回归结果进行预测
        show :
            展示线性回归结果分析表和拟合效果图
        """        
        # 传入参数
        self.data=data                      # 原始数据
        self.vars=vars                      # 作为变量的数据
        self.dependent_var=dependent_var    # 因变量
        self.remove_var=remove_var          # 从原始数据中移除自变量
        # 统计自变量
        self.independent_var=independent_var
        self.independent_var=self.get_independent_var() # 自变量列表
        self.i_var_str=self.independent_var_toString()  # 处理自变量成字符串
        # 回归操作
        self.fit=smf.ols(f'{self.dependent_var}~{self.i_var_str}',data=self.vars).fit()
        self.Summary=self.fit.summary()     # 简报
    def summary(self):
        """
        展示回归分析结果表
        """        
        print(self.Summary)
    def show_lmplot(self,var=None,hue=None,col=None,row=None,col_wrap=None,order=1):
        """
        展示因变量与某个自变量之间的拟合效果图

        Parameters
        ----------
        var : str or list, optional
            选用于做图的自变量，留空时默认为全部自变量
        hue : _type_, optional
            子图变量1
        col : _type_, optional
            子图变量2
        row : _type_, optional
            子图变量3
        col_wrap : _type_, optional
            每行显示的子图数量
        order : int, optional
            多项式拟合次数，默认为1
        """        
        if type(var)==type(None):
            var=self.independent_var
        if type(var)==str:
            if var in self.independent_var: 
                sns.lmplot(x=var,y=self.dependent_var,data=self.data,hue=hue,col=col,row=row,col_wrap=col_wrap,order=order)
            else:
                raise ValueError('var不在自变量列表中，请检查！')
        elif type(var)==list:
            for elem in var:
                if elem in self.independent_var: 
                    sns.lmplot(x=elem,y=self.dependent_var,data=self.data,hue=hue,col=col,row=row,col_wrap=col_wrap,order=order)
                else:
                    raise ValueError(f'var中元素"{elem}"不在自变量列表中，请检查！')
        else:
            raise TypeError('var的类型只接受str或list')
    def forecast(self,vals):
        """
        根据回归结果进行预测

        参数
        ----------
        vals : array like
            用于预测的各自变量系数

        返回值
        -------
        预测的结果
        """        
        res=self.fit.params['Intercept']
        for val,var in zip(vals,self.independent_var):
            res+=self.fit.params[var]*val
        return res
    def show(self):
        """
        展示线性回归结果分析表和拟合效果图
        """        
        print('线性回归分析结果表：')
        self.summary()
        print('拟合效果图：')
        plt.figure(figsize=(20,8),dpi=80)
        plt.plot(self.data.index,self.data[self.dependent_var],marker='o')
        forecast_list=[]
        for i in range(0,len(self.data.index)):
            sub_list=[]
            for var in self.independent_var:
                sub_list.append(self.data[var][i])
            forecast_list.append(self.forecast(sub_list))
        plt.plot(self.data.index,forecast_list,marker='o')
        plt.show()