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

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

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

1. 学习时间

2020.11.01 到 2020.11.02

2. 学习内容

  1. 参考

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

  1. 最小二乘法
  2. 线性回归方程
  3. Python 编写线性回归方程

3. 学习产出

3.1 正态分布

  1. 为什么正态分布

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

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

  1. 正态分布表达式:

Y∼N(μ,σ2) Y {\sim} N(\mu, \sigma^{2}) YN(μ,σ2)

Y 的概率密度:

f(y)=(12πσ)ne−12σ2[y−μ]2 f(y) = ( { \frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - {\mu}]^2 } f(y)=(2πσ1)ne2σ21[yμ]2

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

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

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

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

代入到 Y 的概率密度
f(y)=(12πσ)ne−12σ2[y−μ(x)]2 f(y) = ( { \frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - { \mu(x) }]^2 } f(y)=(2πσ1)ne2σ21[yμ(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 } f(y)=(2πσ1)ne2σ21[yy^]2
⇓\Downarrow
f(y)=(12πσ)ne−12σ2[y−bTx+a]2 f(y) = ( { \frac {1} {\sqrt{2\pi} {\sigma} } })^n e^{ -{ \frac{1}{2 {\sigma}^2} } [y - { \boldsymbol{ {b}^{T} } \boldsymbol{x} + a }]^2 } f(y)=(2πσ1)ne2σ21[ybTx+a]2







3.2 假设函数

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

x表示{x1,x2,…,xn}T,xi表示第i维向量;bT表示{b1,b2,…,bn}T的转置,表示第i个维度的向量xi前的系数 \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 前的系数 x{x1,x2,,xn}TxiibT{b1,b2,,bn}Tixi



\

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

x表示{x},二维平面向量只有一个维度(加上y^就二维);bT表示{b}的转置,表示第1个维度的向量x前的系数 \boldsymbol{x} 表示 \{x\},二维平面向量只有一个维度(加上 \hat{y}就二维); \boldsymbol{ {b}^{T} } 表示 \{ b \} 的转置,表示第 1 个维度的向量 x 前的系数 x{x}y^bT{b}1x

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

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

x表示{x1,x2}T,xi表示第i维向量;bT表示{b1,b2}T的转置,表示第i个维度的向量xi前的系数 \boldsymbol{x} 表示 \{x_1,x_2\}^{T},x_i 表示 第 i 维向量; \boldsymbol{ {b}^{T} } 表示 \{ {b_1},{b_2} \}^{T} 的转置,表示第 i 个维度的向量 x_i 前的系数 x{x1,x2}TxiibT{b1,b2}Tixi

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

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



3.3 最大似然估计法

  1. 因为在 n 次实验中得到观测值
    (x1,y1),(x2,y2),…,(xn,yn),(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n),(x1,y1),(x2,y2),,(xn,yn),
    所以使用最大似然估计法估计未知参数的估计值。似然函数如下:
    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 }L(b)=L(b1,b2,,bk)=i=1nf(yi)=(2πσ1)ne2σ21i=1n[yiμ(xi)]2

  2. 要使似然函数取得最大值,则上式指数中的平方和
    S=∑i=1n[yi−μ(xi;b1,b2,…,bk)]2S = \sum_{i=1}^n [y_i - \mu(x_i;b_1, b_2, \dots, b_k)]^2 S=i=1n[yiμ(xi;b1,b2,,bk)]2
    取最小值。

3. 意义:

为了使观测值(xi,yi)(i=1,2,3,…,n)出现的可能性最大,应选择参数b1,b2,…,bk,使得观测值yi与相应的函数值μ(xi;b1,b2,…,bk)的偏差平方和最小。 为了使观测值 (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) 的偏差平方和最小。 使(xi,yi)(i=1,2,3,,n)b1,b2,,bk使yiμ(xi;b1,b2,,bk)

  1. 要计算 2 中 S 最小,则需求 S 的极小值。分别求 S 对 其参数求导,得到方程组:
    分别对(b1,b2,…,bk)求偏导:分别对 (b_1, b_2, \dots, b_k) 求偏导:(b1,b2,,bk)
    ∑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)]a1μ(xi;b1,b2,,bn)=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,i=1n[yiμ(xi;b1,b2,,bn)]a2μ(xi;b1,b2,,bn)=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.i=1n[yiμ(xi;b1,b2,,bn)]akμ(xi;b1,b2,,bn)=0.



3.4 线性回归方程

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

Y 和 X 的关系
Y∼N(a^+b^x,σ2),Y \sim N(\hat{a} + \hat{b} x, \sigma^2) ,YN(a^+b^x,σ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=i=1n(yia^b^xi)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) = 0i=1n(a^+b^xiyi)=0
∑i=1n(a^+b^xi−yi)xi=0\sum_{i=1}^n (\hat{a} + \hat{b} x_i - y_i) x_i = 0i=1n(a^+b^xiyi)xi=0


整理得:
na^+(∑i=1nxi)b^=∑i=1nyi n\hat{a} + (\sum_{i=1}^n x_i) \hat{b} = \sum_{i=1}^n y_i na^+(i=1nxi)b^=i=1nyi
(∑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 (i=1nxi)a^+(i=1nxi2)b^=i=1nxiyi


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


其中,
xˉ=1n∑i=1nxi,\bar{x} = \frac {1} {n} \sum_{i=1}^n x_i, xˉ=n1i=1nxi,
yˉ=1n∑i=1nyi, \bar{y} = \frac {1} {n} \sum_{i=1}^n y_i, yˉ=n1i=1nyi,
lxy=∑i=1n(yi−yˉ)(xi−xˉ)=∑i=1n(xiyi−nxˉyˉ), 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}), lxy=i=1n(yiyˉ)(xixˉ)=i=1n(xiyinxˉyˉ),
lxx=∑i=1n(xi−xˉ)(xi−xˉ)=nE[(X−E(X))2]=(n−1)sx2, 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, lxx=i=1n(xixˉ)(xixˉ)=nE[(XE(X))2]=(n1)sx2,
lyy=∑i=1n(yi−yˉ)(yi−yˉ)=nE[(Y−E(Y))2]=(n−1)sy2. 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. lyy=i=1n(yiyˉ)(yiyˉ)=nE[(YE(Y))2]=(n1)sy2.
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) 的样本方差sx2,sy2(x1,x2,,xn),(y1,y2,,yn)
将a^,b^代入方程就可以得到平面上的线性回归方程:y^=a^+b^x 将 \hat{a}, \hat{b} 代入方程 就可以得到平面上的线性回归方程 : \hat{y} = \hat{a} + \hat{b} x a^,b^线:y^=a^+b^x






3.5 例题

  1. 题目 (《概率论与数理统计教程》 第四版 (沈恒范))P242-P243
    在这里插入图片描述
    在这里插入图片描述





  1. 代码求解

准备数据文件(文件名字): 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)



  1. 结果:
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))
实现
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值