目录
1 最小二乘拟合(方法一)
本章介绍的是我在上研究生课程《数值分析》时所学的内容,也是最一般的解法,可以对任意函数进行拟合。
1.1 数学推导
对一组数据
(
x
i
,
y
i
)
,
i
=
1
,
2
⋯
,
m
(x_i,y_i),i=1,2\cdots,m
(xi,yi),i=1,2⋯,m,要在某个函数类
Φ
=
{
ϕ
0
(
x
)
,
ϕ
1
(
x
)
,
⋯
,
ϕ
n
(
x
)
}
,
n
≪
m
\Phi=\{\phi_0(x),\phi_1(x),\cdots,\phi_n(x)\},n\ll m
Φ={ϕ0(x),ϕ1(x),⋯,ϕn(x)},n≪m中构造一个函数
ϕ
(
x
)
=
∑
k
=
0
n
a
k
ϕ
k
(
x
)
\phi(x)=\sum_{k=0}^{n}a_k\phi_k(x)
ϕ(x)=∑k=0nakϕk(x),使得
I
=
∑
i
=
1
m
[
ϕ
(
x
i
)
−
y
i
]
2
I=\sum_{i=1}^{m}\left[ \phi(x_i)-y_i\right]^2
I=∑i=1m[ϕ(xi)−yi]2取得极小值。
此处的函数类
Φ
\Phi
Φ就是指很多种类函数的集合,例如幂函数、三角函数、指数函数、有理函数、多项式函数等。例如,常用的多项式函数类有
Φ
1
=
{
1
,
x
,
x
2
,
⋯
,
x
n
}
\Phi_1=\{1,x,x^2,\cdots,x^n\}
Φ1={1,x,x2,⋯,xn},
Φ
2
=
{
1
,
x
,
x
,
⋯
,
x
⏟
n
}
\Phi_2=\{1,\underbrace{x,x,\cdots,x}_{n}\}
Φ2={1,n
x,x,⋯,x}。其中
Φ
1
\Phi_1
Φ1为n阶多项式,
Φ
2
\Phi_2
Φ2为n元一阶多项式,分别对应着多项式拟合和多维线性回归。
残差函数
I
I
I定义为
I
(
a
0
,
a
1
,
⋯
,
a
n
)
=
∑
i
=
1
m
[
a
0
ϕ
0
(
x
i
)
+
a
1
ϕ
1
(
x
i
)
+
⋯
+
a
n
ϕ
n
(
x
i
)
−
y
i
]
2
I(a_0,a_1,\cdots,a_n)=\sum_{i=1}^{m}\left[a_0\phi_0(x_i)+a_1\phi_1(x_i)+\cdots+a_n\phi_n(x_i) -y_i\right]^2
I(a0,a1,⋯,an)=i=1∑m[a0ϕ0(xi)+a1ϕ1(xi)+⋯+anϕn(xi)−yi]2
令
I
I
I对系数
a
k
a_k
ak的偏导为0,即可求出使残差函数
I
I
I最小的系数
a
k
a_k
ak的值
∂
I
∂
a
k
=
2
∑
i
=
1
m
[
ϕ
(
x
i
)
−
y
i
]
∂
ϕ
(
x
i
)
∂
a
k
=
2
∑
i
=
1
m
[
2
∑
j
=
0
n
a
j
ϕ
j
(
x
i
)
−
y
i
]
ϕ
k
(
x
i
)
=
2
{
∑
j
=
0
n
a
j
[
∑
i
=
1
m
ϕ
j
(
x
i
)
ϕ
k
(
x
i
)
]
−
∑
i
=
1
m
y
i
ϕ
k
(
x
i
)
}
\begin{aligned} \frac{\partial I}{\partial a_k}&=2\sum_{i=1}^m\left[ \phi(x_i)-y_i \right]\frac{\partial \phi(x_i)}{\partial a_k} \\ &=2\sum_{i=1}^m \left[ 2\sum_{j=0}^n a_j\phi_j(x_i) - y_i\right]\phi_k(x_i)\\ &=2\left\{ \sum_{j=0}^n a_j\left[\sum_{i=1}^m \phi_j(x_i)\phi_k(x_i)\right] - \sum_{i=1}^m y_i\phi_k(x_i) \right\} \end{aligned}
∂ak∂I=2i=1∑m[ϕ(xi)−yi]∂ak∂ϕ(xi)=2i=1∑m[2j=0∑najϕj(xi)−yi]ϕk(xi)=2{j=0∑naj[i=1∑mϕj(xi)ϕk(xi)]−i=1∑myiϕk(xi)}
简便起见,记
(
ϕ
j
,
ϕ
k
)
=
∑
i
=
1
m
ϕ
j
(
x
i
)
ϕ
k
(
x
i
)
(
y
,
ϕ
k
)
=
∑
i
=
1
m
y
i
ϕ
k
(
x
i
)
\begin{aligned} (\phi_j,\phi_k)&=\sum_{i=1}^m \phi_j(x_i)\phi_k(x_i)\\ (y,\phi_k)&=\sum_{i=1}^m y_i\phi_k(x_i) \end{aligned}
(ϕj,ϕk)(y,ϕk)=i=1∑mϕj(xi)ϕk(xi)=i=1∑myiϕk(xi)
则
∂
I
∂
a
k
\frac{\partial I}{\partial a_k}
∂ak∂I可写成如下形式
∂
I
∂
a
k
=
2
[
∑
j
=
0
n
a
j
(
ϕ
j
,
ϕ
k
)
−
(
y
,
ϕ
k
)
]
=
0
\frac{\partial I}{\partial a_k} = 2\left[\sum_{j=0}^na_j(\phi_j,\phi_k) - (y,\phi_k)\right] = 0
∂ak∂I=2[j=0∑naj(ϕj,ϕk)−(y,ϕk)]=0
化简为
∑
j
=
0
n
a
j
(
ϕ
j
,
ϕ
k
)
=
(
y
,
ϕ
k
)
\sum_{j=0}^na_j(\phi_j,\phi_k) = (y,\phi_k)
j=0∑naj(ϕj,ϕk)=(y,ϕk)
注意到,公式中的系数
a
k
a_k
ak中
k
=
0
,
1
,
⋯
,
n
k=0,1,\cdots,n
k=0,1,⋯,n,则实际上的操作步骤为分别对每个系数求偏导,并令偏导为0,从而求出使残差函数
I
I
I取到极小值的所有的系数。写成矩阵的形式
[
(
ϕ
0
,
ϕ
0
)
(
ϕ
0
,
ϕ
1
)
⋯
(
ϕ
0
,
ϕ
n
)
(
ϕ
1
,
ϕ
0
)
(
ϕ
1
,
ϕ
1
)
⋯
(
ϕ
1
,
ϕ
n
)
⋮
⋮
⋮
(
ϕ
n
,
ϕ
0
)
(
ϕ
n
,
ϕ
1
)
⋯
(
ϕ
n
,
ϕ
n
)
]
[
a
0
a
1
⋮
a
n
]
=
[
(
y
,
ϕ
0
)
(
y
,
ϕ
1
)
⋮
(
y
,
ϕ
n
)
]
\begin{bmatrix} (\phi_0,\phi_0)& (\phi_0,\phi_1)& \cdots & (\phi_0,\phi_n) \\ (\phi_1,\phi_0)& (\phi_1,\phi_1)&\cdots&(\phi_1,\phi_n)\\ \vdots &\vdots&&\vdots\\ (\phi_n,\phi_0)& (\phi_n,\phi_1)&\cdots&(\phi_n,\phi_n) \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix}= \begin{bmatrix} (y,\phi_0)\\(y,\phi_1)\\ \vdots \\(y,\phi_n) \end{bmatrix}
⎣⎢⎢⎢⎡(ϕ0,ϕ0)(ϕ1,ϕ0)⋮(ϕn,ϕ0)(ϕ0,ϕ1)(ϕ1,ϕ1)⋮(ϕn,ϕ1)⋯⋯⋯(ϕ0,ϕn)(ϕ1,ϕn)⋮(ϕn,ϕn)⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡(y,ϕ0)(y,ϕ1)⋮(y,ϕn)⎦⎥⎥⎥⎤
上式方程组称为法方程组(Normal Equation),向量
[
a
0
a
1
⋮
a
n
]
\begin{bmatrix}a_0\\a_1\\ \vdots \\ a_n\end{bmatrix}
⎣⎢⎢⎢⎡a0a1⋮an⎦⎥⎥⎥⎤称为回归系数(Regression Coefficients)。用更一般的形式将上式再次改写
A
x
=
b
Ax=b
Ax=b
则可通过矩阵求逆的方式求出回归系数
x
=
A
−
1
b
x = A^{-1}b
x=A−1b
1.2 算例
利用一阶多项式 Φ 1 = { 1 , x } \Phi_1=\{1,x\} Φ1={1,x}通过1.1节的方法对如下数据进行最小二乘拟合。
Independent variable x | Dependent variable y |
---|---|
0 | 394.33 |
4 | 329.50 |
8 | 291.00 |
12 | 255.17 |
16 | 229.33 |
20 | 204.83 |
24 | 179.00 |
28 | 163.83 |
32 | 150.33 |
则法方程组为
[
9
∑
x
i
∑
x
i
∑
x
i
2
]
[
a
0
a
1
]
=
[
∑
y
i
∑
x
i
y
i
]
\begin{bmatrix} 9&\sum x_i\\\sum x_i&\sum x_i^2 \end{bmatrix} \begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} \sum y_i \\ \sum x_iy_i \end{bmatrix}
[9∑xi∑xi∑xi2][a0a1]=[∑yi∑xiyi]
带入数据
[
9
144
144
3264
]
[
a
0
a
1
]
=
[
2197.32
28167.72
]
\begin{bmatrix} 9&144\\144&3264 \end{bmatrix} \begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} 2197.32 \\ 28167.72 \end{bmatrix}
[91441443264][a0a1]=[2197.3228167.72]
求解得到
[
a
0
a
1
]
=
[
360.636667
−
7.280625
]
\begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} 360.636667\\-7.280625 \end{bmatrix}
[a0a1]=[360.636667−7.280625]
则可得到线性回归方程
y
=
−
7.280625
x
+
360.636667
y = -7.280625x + 360.636667
y=−7.280625x+360.636667
拟合结果如下
1.3 Python 代码
#%% import neccesary packages
import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset
x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%% calculate the normal equation
A11 = 9
A12 = np.sum(x)
A21 = A12
A22 = np.sum(np.power(x,2))
A = np.array([[A11,A12],[A21,A22]])
print(A)
b1 = np.sum(y)
b2 = np.sum(x*y)
b=np.array([b1, b2])
print(b)
#%% solve the polynamial coefficient a0 & a1
a = np.linalg.inv(A).dot(b)
a0 = a[0]
a1 = a[1]
print(a)
#%% plot the picture
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, a0 + a1 * x)
plt.show()
2.最小二乘拟合(方法二)
本章介绍我自己的推导思路,本文以多项式拟合为例说明思路,其他函数的拟合思路相同。
2.1 数学推导
对一组数据
(
x
i
,
y
i
)
,
i
=
1
,
2
⋯
,
m
(x_i,y_i),i=1,2\cdots,m
(xi,yi),i=1,2⋯,m进行多项式拟合,假设多项式为
f
(
x
)
=
a
n
x
n
+
a
n
−
1
x
n
−
1
+
⋯
+
a
1
x
+
a
0
f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0
f(x)=anxn+an−1xn−1+⋯+a1x+a0
对任意数据点
(
x
k
,
y
k
)
(x_k,y_k)
(xk,yk),它与多项式
f
(
x
)
f(x)
f(x)在
y
y
y轴方向上的距离称为残差
δ
k
\delta_k
δk,定义为
δ
k
=
f
(
x
k
)
−
y
k
\delta_k=f(x_k)-y_k
δk=f(xk)−yk
拟合的目标是寻找一组多项式系数
{
a
0
,
a
1
,
⋯
,
a
n
}
\{a_0,a_1,\cdots,a_n\}
{a0,a1,⋯,an}使得残差函数
I
I
I取最小值,本章将直接利用最小二乘法求解。将每个数据点的残差写成矩阵形式
[
1
x
1
x
1
2
⋯
x
1
n
1
x
2
x
2
2
⋯
x
2
n
⋮
⋮
⋮
⋮
1
x
m
x
m
2
⋯
x
m
n
]
[
a
0
a
1
⋮
a
n
]
−
[
y
1
y
2
⋮
y
m
]
=
[
δ
1
δ
2
⋮
δ
m
]
\begin{bmatrix} 1& x_1& x^2_1&\cdots & x^n_1 \\ 1&x_2& x^2_2&\cdots&x^n_2\\ \vdots &\vdots &\vdots&&\vdots\\ 1&x_m& x^2_m&\cdots&x^n_m \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix}- \begin{bmatrix} y_1\\y_2\\ \vdots \\y_m \end{bmatrix} = \begin{bmatrix} \delta_1\\\delta_2\\ \vdots \\\delta_m \end{bmatrix}
⎣⎢⎢⎢⎡11⋮1x1x2⋮xmx12x22⋮xm2⋯⋯⋯x1nx2n⋮xmn⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an⎦⎥⎥⎥⎤−⎣⎢⎢⎢⎡y1y2⋮ym⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡δ1δ2⋮δm⎦⎥⎥⎥⎤
定义残差函数
I
(
a
0
,
a
1
,
⋯
,
a
n
)
I(a_0,a_1,\cdots,a_n)
I(a0,a1,⋯,an)为
I
=
∥
δ
i
∥
2
2
=
∥
f
(
x
i
)
−
y
i
∥
2
2
I=\parallel\delta_i\parallel ^2_2=\parallel f(x_i) - y_i \parallel ^2_2
I=∥δi∥22=∥f(xi)−yi∥22
要使残差函数
I
I
I最小,本质上是求如下方程的最小二乘解
[
1
x
1
x
1
2
⋯
x
1
n
1
x
2
x
2
2
⋯
x
2
n
⋮
⋮
⋮
⋮
1
x
m
x
m
2
⋯
x
m
n
]
[
a
0
a
1
⋮
a
n
]
=
[
y
1
y
2
⋮
y
m
]
\begin{bmatrix} 1& x_1& x^2_1&\cdots & x^n_1 \\ 1&x_2& x^2_2&\cdots&x^n_2\\ \vdots &\vdots &\vdots&&\vdots\\ 1&x_m& x^2_m&\cdots&x^n_m \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots \\ a_n \end{bmatrix}= \begin{bmatrix} y_1\\y_2\\ \vdots \\y_m \end{bmatrix}
⎣⎢⎢⎢⎡11⋮1x1x2⋮xmx12x22⋮xm2⋯⋯⋯x1nx2n⋮xmn⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y1y2⋮ym⎦⎥⎥⎥⎤
上式可写成一般形式
A
x
=
b
Ax=b
Ax=b
(1)当
A
A
A为方阵时,其最小二乘解为
x
=
−
(
A
T
A
)
−
1
A
T
b
x=-(A^TA)^{-1}A^Tb
x=−(ATA)−1ATb
(2)在大多数情况下,
A
A
A不为方阵,可用Moore-Penrose逆
A
+
A^+
A+求解
定理1: 设
A
∈
C
m
×
n
A\in \mathbb{C} ^{m \times n}
A∈Cm×n是秩为
r
(
r
>
0
)
r(r>0)
r(r>0)的矩阵,且
A
A
A的满秩分解为
A
=
F
G
(
F
∈
C
m
×
r
,
G
∈
C
r
×
n
秩
都
为
r
)
A=FG\quad (F\in \mathbb{C} ^{m \times r},G\in \mathbb{C} ^{r \times n}秩都为r)
A=FG(F∈Cm×r,G∈Cr×n秩都为r)
则
A
+
=
G
H
(
G
G
H
)
−
1
(
F
H
F
)
−
1
F
H
A^+=G^H(GG^H)^{-1}(F^HF)^{-1}F^H
A+=GH(GGH)−1(FHF)−1FH
定理2: 设
A
∈
C
m
×
n
,
b
∈
C
m
A\in \mathbb{C} ^{m \times n},b\in \mathbb{C} ^{m}
A∈Cm×n,b∈Cm有解的充要条件是
A
A
+
b
=
b
AA^+b=b
AA+b=b
其通解为
x
=
A
+
b
+
(
I
−
A
+
A
)
y
x=A^+b+(I-A^+A)y
x=A+b+(I−A+A)y
其唯一极小范数最小二乘解为
x
=
A
+
b
x=A^+b
x=A+b
2.2 算例
本章算例与1.2节相同,以作对照。
利用一阶多项式
y
=
a
0
+
a
1
x
y=a_0+a_1x
y=a0+a1x对如下数据进行最小二乘拟合。
Independent variable x | Dependent variable y |
---|---|
0 | 394.33 |
4 | 329.50 |
8 | 291.00 |
12 | 255.17 |
16 | 229.33 |
20 | 204.83 |
24 | 179.00 |
28 | 163.83 |
32 | 150.33 |
列写
A
x
=
b
Ax=b
Ax=b的一般形式
[
1
0
1
4
1
8
1
12
1
16
1
20
1
24
1
28
1
32
]
[
a
0
a
1
]
=
[
394.33
329.50
291.00
255.17
229.33
204.83
179.00
163.83
150.33
]
\begin{bmatrix} 1& 0\\ 1& 4\\ 1& 8\\ 1& 12\\ 1& 16\\ 1& 20\\ 1& 24\\ 1& 28\\ 1& 32 \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \end{bmatrix}= \begin{bmatrix} 394.33\\ 329.50\\ 291.00\\ 255.17\\ 229.33\\ 204.83\\ 179.00\\ 163.83\\ 150.33 \end{bmatrix}
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡111111111048121620242832⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤[a0a1]=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡394.33329.50291.00255.17229.33204.83179.00163.83150.33⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
利用
P
y
t
h
o
n
Python
Python命令
n
u
m
p
y
.
l
i
n
a
l
g
.
p
i
n
v
(
)
numpy.linalg.pinv()
numpy.linalg.pinv()求解矩阵
A
A
A的
A
+
A^+
A+逆,并求其唯一极小范数最小二乘解
[
a
0
a
1
]
=
[
360.636667
−
7.280625
]
\begin{bmatrix} a_0\\a_1 \end{bmatrix}= \begin{bmatrix} 360.636667\\-7.280625 \end{bmatrix}
[a0a1]=[360.636667−7.280625]
求解结果与第一章相同。
2.3 Python 代码
#%% import neccesary packages
import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset
x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%% solve the polynamial coefficient a0 & a1
A = np.vstack((np.array([1]*9), x))
print(A)
pinv_A = np.linalg.pinv(A).T
print(pinv_A)
a = pinv_A.dot(y)
print(a)
a0 = a[0]
a1 = a[1]
#%% plot the picture
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, a0 + a1 * x)
plt.show()
3 最小二乘法拟合(方法三)
本章重点介绍线性回归算法。其实该算法本质上和第一章方法一样,只不过将其限定为一阶多项式,推导出更直观的表达方法。
3.1 数学推导
假设线性回归模型如下
y
i
=
β
0
+
β
1
x
i
+
ϵ
i
y_i = \beta_0+\beta_1x_i+\epsilon_i
yi=β0+β1xi+ϵi
式中
ϵ
i
\epsilon_i
ϵi为模型噪声服从高斯分布:
ϵ
∼
N
(
0
,
σ
2
)
\epsilon \sim N(0,\sigma^2)
ϵ∼N(0,σ2)。
假设需要对一组数据
(
x
i
,
y
i
)
,
i
=
1
,
2
⋯
,
m
(x_i,y_i),i=1,2\cdots,m
(xi,yi),i=1,2⋯,m进行线性回归,定义残差函数
Q
=
∑
i
=
1
n
[
y
i
−
(
β
0
+
β
1
x
i
)
]
2
Q = \sum_{i=1}^n \left[ y_i-(\beta_0+\beta_1x_i) \right]^2
Q=i=1∑n[yi−(β0+β1xi)]2
优化目标为最小化残差函数
m
i
n
β
0
,
β
1
:
Q
\mathop{min}\limits_{\beta_0,\beta_1}: Q
β0,β1min:Q
和第一章方法一样,令残差函数
Q
Q
Q对系数
β
0
、
β
1
\beta_0、\beta_1
β0、β1的偏导为0
∂
Q
∂
β
0
=
−
2
∑
i
=
1
n
[
y
i
−
(
β
0
+
β
1
x
i
)
]
=
0
⇒
n
β
0
+
β
1
∑
i
=
1
n
x
i
=
∑
i
=
1
n
y
i
(
1
)
\begin{aligned} \frac{\partial Q}{\partial \beta_0}=-2\sum_{i=1}^{n}\left[ y_i-(\beta_0+\beta_1x_i) \right]&=0 \\ \Rightarrow \qquad n\beta_0+\beta_1\sum_{i=1}^nx_i &=\sum_{i=1}^ny_i \qquad (1) \end{aligned}
∂β0∂Q=−2i=1∑n[yi−(β0+β1xi)]⇒nβ0+β1i=1∑nxi=0=i=1∑nyi(1)
∂
Q
∂
β
1
=
−
2
∑
i
=
1
n
x
i
[
y
i
−
(
β
0
+
β
1
x
i
)
]
=
0
⇒
β
0
∑
i
=
1
n
x
i
+
β
1
∑
i
=
1
n
x
i
2
=
β
1
∑
i
=
1
n
x
i
y
i
(
2
)
\begin{aligned} \frac{\partial Q}{\partial \beta_1}=-2\sum_{i=1}^{n}x_i\left[ y_i-(\beta_0+\beta_1x_i) \right]&=0 \\ \Rightarrow \qquad \beta_0\sum_{i=1}^nx_i+\beta_1\sum_{i=1}^nx_i^2 &=\beta_1\sum_{i=1}^nx_iy_i \qquad (2) \end{aligned}
∂β1∂Q=−2i=1∑nxi[yi−(β0+β1xi)]⇒β0i=1∑nxi+β1i=1∑nxi2=0=β1i=1∑nxiyi(2)
联立式(1)和式(2)解方程得到
β
0
、
β
1
\beta_0、\beta_1
β0、β1的代数表达式
β
^
0
=
(
∑
i
=
1
n
x
i
2
)
(
∑
i
=
1
n
y
i
)
−
(
∑
i
=
1
n
x
i
)
(
∑
i
=
1
n
x
i
y
i
)
n
∑
i
=
1
n
x
i
2
−
(
∑
i
=
1
n
x
i
)
2
\hat{\beta}_0=\frac{(\sum_{i=1}^nx_i^2)(\sum_{i=1}^ny_i)-(\sum_{i=1}^nx_i )(\sum_{i=1}^nx_iy_i)}{ n\sum_{i=1}^nx_i^2-(\sum_{i=1}^nx_i)^2}
β^0=n∑i=1nxi2−(∑i=1nxi)2(∑i=1nxi2)(∑i=1nyi)−(∑i=1nxi)(∑i=1nxiyi)
β
^
1
=
∑
i
=
1
n
x
i
y
i
−
(
∑
i
=
1
n
y
i
)
(
∑
i
=
1
n
x
i
)
n
∑
i
=
1
n
x
i
2
−
(
∑
i
=
1
n
x
i
)
2
\hat{\beta}_1=\frac{\sum_{i=1}^nx_iy_i-(\sum_{i=1}^ny_i)(\sum_{i=1}^nx_i )}{ n\sum_{i=1}^nx_i^2-(\sum_{i=1}^nx_i)^2}
β^1=n∑i=1nxi2−(∑i=1nxi)2∑i=1nxiyi−(∑i=1nyi)(∑i=1nxi)
为表达简便,记
S
x
y
=
∑
i
=
1
n
x
i
y
i
−
1
n
(
∑
i
=
1
n
x
i
)
(
∑
i
=
1
n
y
i
)
=
∑
i
=
1
n
(
x
i
−
x
‾
)
(
y
i
−
y
‾
)
S_{xy}=\sum_{i=1}^nx_iy_i-\frac{1}{n}(\sum_{i=1}^nx_i)(\sum_{i=1}^ny_i)=\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})
Sxy=i=1∑nxiyi−n1(i=1∑nxi)(i=1∑nyi)=i=1∑n(xi−x)(yi−y)
S
x
x
=
∑
i
=
1
n
x
i
2
−
1
n
(
∑
i
=
1
n
x
i
)
2
=
∑
i
=
1
n
(
x
i
−
x
‾
)
2
S_{xx}=\sum_{i=1}^nx_i^2-\frac{1}{n}(\sum_{i=1}^nx_i)^2=\sum_{i=1}^n(x_i-\overline{x})^2
Sxx=i=1∑nxi2−n1(i=1∑nxi)2=i=1∑n(xi−x)2
S
y
y
=
∑
i
=
1
n
y
i
2
−
1
n
(
∑
i
=
1
n
y
i
)
2
=
∑
i
=
1
n
(
y
i
−
y
‾
)
2
S_{yy}=\sum_{i=1}^ny_i^2-\frac{1}{n}(\sum_{i=1}^ny_i)^2=\sum_{i=1}^n(y_i-\overline{y})^2
Syy=i=1∑nyi2−n1(i=1∑nyi)2=i=1∑n(yi−y)2
其中,
x
‾
、
y
‾
\overline{x}、\overline{y}
x、y分别为拟合数据
x
i
、
y
i
x_i、y_i
xi、yi的平均值,
i
=
1
,
2
,
⋯
,
m
i=1,2,\cdots,m
i=1,2,⋯,m
x
‾
=
1
n
∑
i
=
1
n
x
i
y
‾
=
1
n
∑
i
=
1
n
y
i
\begin{aligned} \overline{x}&= \frac{1}{n}\sum_{i=1}^nx_i \\ \overline{y}&= \frac{1}{n}\sum_{i=1}^ny_i \end{aligned}
xy=n1i=1∑nxi=n1i=1∑nyi
则
β
0
、
β
1
\beta_0、\beta_1
β0、β1的表达式可改写为
β
0
^
=
y
‾
−
β
1
x
‾
β
1
^
=
S
x
y
S
x
x
\begin{aligned} \hat{\beta_0}&= \overline{y}-\beta_1\overline{x} \\ \hat{\beta_1}&= \frac{S_{xy}}{S_{xx}} \end{aligned}
β0^β1^=y−β1x=SxxSxy
3.2 算例
本章算例与1.2节相同,以作对照。
对如下数据进行线性回归。
Independent variable x | Dependent variable y |
---|---|
0 | 394.33 |
4 | 329.50 |
8 | 291.00 |
12 | 255.17 |
16 | 229.33 |
20 | 204.83 |
24 | 179.00 |
28 | 163.83 |
32 | 150.33 |
计算系数
β
^
0
、
β
^
1
\hat{\beta}_0、\hat{\beta}_1
β^0、β^1
x
‾
=
16
y
‾
=
244.15
S
x
y
=
−
6989.40
S
x
x
=
960
β
1
^
=
S
x
y
S
x
x
=
−
7.280625
β
0
^
=
y
‾
−
β
1
x
‾
=
360.636667
\begin{aligned} \overline{x}&=16 \\ \overline{y}&=244.15 \\ S_{xy}&=-6989.40 \\ S_{xx}&=960 \\ \hat{\beta_1}&= \frac{S_{xy}}{S_{xx}} = -7.280625 \\ \hat{\beta_0}&= \overline{y}-\beta_1\overline{x}=360.636667 \end{aligned}
xySxySxxβ1^β0^=16=244.15=−6989.40=960=SxxSxy=−7.280625=y−β1x=360.636667
计算结果和第一章、第二章结果一致。
y
=
360.636667
−
7.280625
x
y = 360.636667-7.280625x
y=360.636667−7.280625x
3.3 Python 代码
#%% import neccesary packages
import numpy as np
import matplotlib.pyplot as plt
#%% generate the dataset
x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%%
x_bar = np.sum(x)/9
y_bar = np.sum(y)/9
Sxy = np.sum((x-x_bar)*(y-y_bar))
Sxx = np.sum(np.power(x-x_bar,2))
beta_1 = Sxy/Sxx
beta_0 = y_bar - x_bar * beta_1
print(beta_0)
print(beta_1)
#%%
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x, beta_0 + beta_1 * x)
plt.show()
4 利用sklearn.linear_model()
4.1 参考资料
官网链接:Scikit-learn官网链接
教学1:python-sklearn学习笔记 第一节 linear_model
教学2:numpy中newaxis函数的基本用法及其理解(傻瓜版)
4.2 Python 代码
#%% import neccesary packages
import numpy as np
import matplotlib.pyplot as plt
import sklearn.linear_model as lm
#%% generate the dataset
x = np.arange(0,33,4)
y = np.array([394.33, 329.50, 291.00, 255.17, 229.33, 204.83, 179.00, 163.83, 150.33])
print(x)
#%%
x_lr = np.linspace(0,40,500)
lr = lm.LinearRegression()
lr.fit(x[:, np.newaxis],y)
y_lr = lr.predict(x_lr[:, np.newaxis])
#%%
fig, ax = plt.subplots(1, 1, figsize=(6, 3))
ax.scatter(x, y, color='black')
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
ax.plot(x_lr, y_lr)