算法 002. 线性回归方程的推导及其代码实现

该博客记录了2020.11.01 - 2020.11.02的学习内容，参考《概率论与数理统计教程》，涉及最小二乘法、线性回归方程。介绍了正态分布、假设函数、最大似然估计法，推导了二维线性回归方程，并给出例题求解，结果与书本相近，可修改代码调整小数位。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

1. 学习时间

2020.11.01 到 2020.11.02

2. 学习内容

参考

《概率论与数理统计教程》第四版（沈恒范） chapter 9.1、chapter 9.2

最小二乘法
线性回归方程
Python 编写线性回归方程

3. 学习产出

3.1 正态分布

为什么正态分布

中心极限定理说，在适当的条件下，大量相互独立随机变量的均值经适当标准化后依分布收敛于正态分布，误差的分布就应该是正态分布

参考： https://en.wikipedia.org/wiki/Normal_distribution

正态分布表达式：

${\sim} N(\mu, \sigma^{2})$

Y 的概率密度:

$\frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - {\mu}]^2 }$

在这里插入图片描述
\
\
\

记假设函数为
$μ(x)\mu(x)$

即下面 3.2 的假设函数
$y^=bTx+a\hat{y} = \boldsymbol{ {b}^{T} } \boldsymbol{x} + a$

当然其他情况下假设函数可以是其它形式，在线性回归方程这里就是
$μ(x)=y^=bTx+a\mu(x) = \hat{y} = \boldsymbol{ {b}^{T} } \boldsymbol{x} + a$

代入到 Y 的概率密度
$\frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - { \mu(x) }]^2 }$
$⇓\Downarrow$
$f(y)=(12πσ)ne−12σ2[y−y^]2 f(y) = ( { \frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - { \hat{y} }]^2 }$
$⇓\Downarrow$
$\frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - { \boldsymbol{ {b}^{T} } \boldsymbol{x} + a }]^2 }$

3.2 假设函数

n 维情况下
$y^=bTx+a \hat{y} = \boldsymbol{ {b}^{T} } \boldsymbol{x} + a$
其中，

$\boldsymbol{x} 表示 \{x_1,x_2, \dots,x_n\}^{T}，x_i 表示第 i 维向量； \boldsymbol{ {b}^{T} } 表示 \{ {b_1},{b_2},\dots,{b_n} \}^{T} 的转置，表示第 i 个维度的向量 x_i 前的系数$

二维平面情况
$y^=bTx+a \hat{y} = \boldsymbol{ {b}^{T} } \boldsymbol{x} + a$
其中，

$x表示{x}，二维平面向量只有一个维度（加上y^就二维）；bT表示{b}的转置，表示第1个维度的向量x前的系数 \boldsymbol{x} 表示 \{x\}，二维平面向量只有一个维度（加上 \hat{y}就二维）； \boldsymbol{ {b}^{T} } 表示 \{ b \} 的转置，表示第 1 个维度的向量 x 前的系数$

简单地写就是一条平面曲线
$y^=b^x+a\hat{y} = \hat{b}x+a$

三维情况下
$y^=bTx+a \hat{y} = \boldsymbol{ {b}^{T} } \boldsymbol{x} + a$
其中，

$\boldsymbol{x} 表示 \{x_1,x_2\}^{T}，x_i 表示第 i 维向量； \boldsymbol{ {b}^{T} } 表示 \{ {b_1},{b_2} \}^{T} 的转置，表示第 i 个维度的向量 x_i 前的系数$

简单地写就是
$y^=b1^x1+b2^x2+a\hat{y} = \hat{b_1} x_{1} + \hat{b_2} x_{2}+a$

换成常见的形式就是一条空间曲线
$z^=c^x+b^y+a\hat{z} = \hat{c}x + \hat{b}y + a$

3.3 最大似然估计法

因为在 n 次实验中得到观测值
$(x1,y1),(x2,y2),…,(xn,yn),(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n),$
所以使用最大似然估计法估计未知参数的估计值。似然函数如下：
$L(b)=L(b1,b2,…,bk)=∏i=1nf(yi)=(12πσ)ne−12σ2∑i=1n[yi−μ(xi)]2L(\boldsymbol{ b }) = L(b_1, b_2, \dots, b_k) = \prod_{i=1}^{n} f(y_i) = (\frac 1 {\sqrt{2 \pi} \sigma}) ^n e^{-{\frac 1 {2 {\sigma}^2} } \sum_{i=1}^n [y_i - {\mu(x_i)}] ^2 }$
要使似然函数取得最大值，则上式指数中的平方和
$\sum_{i=1}^n [y_i - \mu(x_i;b_1, b_2, \dots, b_k)]^2$
取最小值。

3. 意义：

$(x_i, y_i) (i = 1, 2, 3, \dots, n) 出现的可能性最大，应选择参数 b_1, b_2, \dots, b_k，使得观测值 y_i 与相应的函数值 \mu(x_i; b_1, b_2, \dots, b_k) 的偏差平方和最小。$

要计算 2 中 S 最小，则需求 S 的极小值。分别求 S 对其参数求导，得到方程组：
$(b_1, b_2, \dots, b_k) 求偏导：$
$∑i=1n[yi−μ(xi;b1,b2,…,bn)]∂∂a1μ(xi;b1,b2,…,bn)=0,\sum_{i=1}^n [y_i - \mu(x_i; b_1, b_2, \dots, b_n)] \frac {\partial} {\partial a_1} \mu(x_i; b_1, b_2, \dots, b_n) = 0,$
$∑i=1n[yi−μ(xi;b1,b2,…,bn)]∂∂a2μ(xi;b1,b2,…,bn)=0,\sum_{i=1}^n [y_i - \mu(x_i; b_1, b_2, \dots, b_n)] \frac {\partial} {\partial a_2} \mu(x_i; b_1, b_2, \dots, b_n) = 0,$
$\dots \dots \dots \dots \dots \dots$
$∑i=1n[yi−μ(xi;b1,b2,…,bn)]∂∂akμ(xi;b1,b2,…,bn)=0.\sum_{i=1}^n [y_i - \mu(x_i; b_1, b_2, \dots, b_n)] \frac {\partial} {\partial a_k} \mu(x_i; b_1, b_2, \dots, b_n) = 0.$

3.4 线性回归方程

只推导二维情况下（即平面下）的线性回归方程：
$y^=b^x+a^\hat{y} = \hat{b}x+\hat{a}$

Y 和 X 的关系
$Y∼N(a^+b^x,σ2),Y \sim N(\hat{a} + \hat{b} x, \sigma^2) ,$

按最小二乘法求解a, b:
$S=∑i=1n(yi−a^−b^xi)2. S = \sum_{i=1}^n (y_i - \hat{a} - \hat{b} x_i)^2.$

为了使 S 取最小值，分别求 S 对 a，b 求偏导，并令他们等于 0，得到方程组：
$∑i=1n(a^+b^xi−yi)=0\sum_{i=1}^n (\hat{a} + \hat{b} x_i - y_i) = 0$
$∑i=1n(a^+b^xi−yi)xi=0\sum_{i=1}^n (\hat{a} + \hat{b} x_i - y_i) x_i = 0$

整理得：
$na^+(∑i=1nxi)b^=∑i=1nyi n\hat{a} + (\sum_{i=1}^n x_i) \hat{b} = \sum_{i=1}^n y_i$
$(∑i=1nxi)a^+(∑i=1nxi2)b^=∑i=1nxiyi (\sum_{i=1}^n x_i) \hat{a} + (\sum_{i=1}^nx_i^2) \hat{b} = \sum_{i=1}^n x_i y_i$

解得：
$a^=yˉ−b^xˉ,\hat{a} = \bar{y} - \hat{b} \bar{x},$
$b^=lxylxx\hat{b} = \frac {l_{xy}} {l_{xx}}$

其中，
$xˉ=1n∑i=1nxi,\bar{x} = \frac {1} {n} \sum_{i=1}^n x_i,$
$\bar{y} = \frac {1} {n} \sum_{i=1}^n y_i,$
$l_{xy} = \sum_{i=1}^n (y_i-\bar{y}) (x_i-\bar{x}) = \sum_{i=1} ^ n (x_iy_i - n \bar{x} \bar{y}),$
$l_{xx} = \sum_{i=1}^n (x_i-\bar{x}) (x_i-\bar{x}) = n E[(X-E(X) )^2 ] = (n-1) s_x^2,$
$l_{yy} = \sum_{i=1}^n (y_i-\bar{y}) (y_i-\bar{y}) = n E[(Y-E(Y) )^2 ] = (n-1) s_y^2.$
$sx2,sy2分别是观测值(x1,x2,…,xn),(y1,y2,…,yn)的样本方差s_x^2, s_y^2 分别是观测值 (x_1, x_2, \dots, x_n), (y_1, y_2, \dots, y_n) 的样本方差$
$将a^,b^代入方程就可以得到平面上的线性回归方程:y^=a^+b^x 将 \hat{a}, \hat{b} 代入方程就可以得到平面上的线性回归方程 : \hat{y} = \hat{a} + \hat{b} x$

3.5 例题

题目（《概率论与数理统计教程》第四版（沈恒范））P242-P243

代码求解

准备数据文件（文件名字）: thickness.dat

0       30.0
1       29.1
2       28.4
3       28.1
4       28.0
5       27.7
6       27.5
7       27.2
8       27.0
9       26.8
10      26.5
11      26.3
12      26.1
13      25.7
14      25.3
15      24.8
16      24.0

"""
file: linear_regresion.py

"""

import numpy as np
from matplotlib import pyplot as plt


file_name = "thickness.dat"



def linear_regresion(x_datas, y_datas):

    # get the average for x, y
    x_bar = np.sum(x_datas) / len(x_datas)
    y_bar = np.sum(y_datas) / len(y_datas)

    # calculate lxx, lxy, lyy
    lxy = np.sum((x_datas - x_bar) * (y_datas - y_bar))
    lxx = np.sum((x_datas - x_bar) * (x_datas - x_bar))
    lyy = np.sum((y_datas - y_bar) * (y_datas - y_bar))

    # calculate coefficient a, b
    b = lxy / lxx
    a = y_bar - (b * x_bar)


    return (a, b)


# plotting lines to check results
# plot init line and fitted line
def paint_lines(x_datas, y_datas, a, b):

    X = x_datas
    Y_1 = y_datas
    Y_2 = b * X + a

    # plot init line
    plt.plot(X, Y_1, 'wheat', marker=".", ms=15, label="original data")
    # plot fitted line
    plt.plot(X, Y_2, 'r', label="fitted line", linewidth=3)

    plt.xlabel("time / h")
    plt.ylabel("thickness / cm")
    plt.legend(loc="upper right")
    plt.title("thickness-time")
    plt.axis([0, 17, 23, 32])

    plt.savefig("2.png")
    plt.show()



if __name__ == "__main__":


    # get all datas
    datas = np.loadtxt(file_name)
    x_col, y_col = (0, 1)

    # get datas for x label and y label
    x_datas = datas[:, x_col]
    y_datas = datas[:, y_col]

    a, b = linear_regresion(x_datas, y_datas)

    print("a:  " + str(a) + '\n' + "b:  " + str(b))
    print("formula:  y = {b}x  + {a} ".format(b=b, a=a))

    paint_lines(x_datas, y_datas, a, b)

结果：

In [54]: ! python linear_regression.py
a:  29.380392156862744
b:  -0.3012254901960784
formula:  y = -0.3012254901960784x  + 29.380392156862744

“2.png” 图如下

在这里插入图片描述

4. 可以看到，结果跟书本给出的一样，只是小数位不一样。这个可以通过修改第 71 行代码

print("formula:  y = {b:5.3f}x  + {a:5.2f} ".format(b=b, a=a))