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
- 正态分布表达式:
Y∼N(μ,σ2) Y {\sim} N(\mu, \sigma^{2}) Y∼N(μ,σ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)ne−2σ21[y−μ]2
\
\
\
- 记假设函数为
μ(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)ne−2σ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)ne−2σ21[y−y^]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)ne−2σ21[y−bTx+a]2
3.2 假设函数
- 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}T,xi表示第i维向量;bT表示{b1,b2,…,bn}T的转置,表示第i个维度的向量xi前的系数
\
- 二维平面情况
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}的转置,表示第1个维度的向量x前的系数
简单地写就是一条平面曲线
y^=b^x+a\hat{y} = \hat{b}x+ay^=b^x+a
- 三 维情况下
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}T,xi表示第i维向量;bT表示{b1,b2}T的转置,表示第i个维度的向量xi前的系数
简单地写就是
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 最大似然估计法
-
因为在 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=1∏nf(yi)=(2πσ1)ne−2σ21∑i=1n[yi−μ(xi)]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=1∑n[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)的偏差平方和最小。
- 要计算 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=1∑n[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=1∑n[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=1∑n[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) ,Y∼N(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=1∑n(yi−a^−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=1∑n(a^+b^xi−yi)=0
∑i=1n(a^+b^xi−yi)xi=0\sum_{i=1}^n (\hat{a} + \hat{b} x_i - y_i) x_i = 0i=1∑n(a^+b^xi−yi)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=1∑nxi)b^=i=1∑nyi
(∑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=1∑nxi)a^+(i=1∑nxi2)b^=i=1∑nxiyi
解得:
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=1∑nxi,
yˉ=1n∑i=1nyi, \bar{y} = \frac {1} {n} \sum_{i=1}^n y_i, yˉ=n1i=1∑nyi,
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=1∑n(yi−yˉ)(xi−xˉ)=i=1∑n(xiyi−nxˉ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=1∑n(xi−xˉ)(xi−xˉ)=nE[(X−E(X))2]=(n−1)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=1∑n(yi−yˉ)(yi−yˉ)=nE[(Y−E(Y))2]=(n−1)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 例题
- 题目 (《概率论与数理统计教程》 第四版 (沈恒范))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))