1 最小二乘拟合（方法一）
1.1 数学推导1.2 算例1.3 Python 代码 2.最小二乘拟合（方法二）
2.1 数学推导2.2 算例2.3 Python 代码 3 最小二乘法拟合（方法三）
3.1 数学推导3.2 算例3.3 Python 代码 4 利用sklearn.linear_model()
4.1 参考资料4.2 Python 代码
1 最小二乘拟合（方法一）本章介绍的是我在上研究生课程《数值分析》时所学的内容，也是最一般的解法，可以对任意函数进行拟合。
1.1 数学推导对一组数据
(
x
i
,
y
i
)
,
i
=
1
,
2
⋯
,
m
(x_i,y_i),i=1,2cdots,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)},nll m
Φ={ϕ0(x),ϕ1(x),⋯,ϕn(x)},n≪m中构造一个函数
ϕ
(
x
)
=
∑
k
=
0
n
a
k
ϕ
k
(
x
)
phi(x)=sum_{k=0}^{n}a_kphi_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_iright]^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_0phi_0(x_i)+a_1phi_1(x_i)+cdots+a_nphi_n(x_i) -y_iright]^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}&=2sum_{i=1}^mleft[ phi(x_i)-y_i right]frac{partial phi(x_i)}{partial a_k} \ &=2sum_{i=1}^m left[ 2sum_{j=0}^n a_jphi_j(x_i) - y_iright]phi_k(x_i)\ &=2left{ sum_{j=0}^n a_jleft[sum_{i=1}^m phi_j(x_i)phi_k(x_i)right] - sum_{i=1}^m y_iphi_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_iphi_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} = 2left[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_nend{bmatrix}
⎣⎢⎢⎢⎡a0a1⋮an⎦⎥⎥⎥⎤称为回归系数（Regression Coefficients）。用更一般的形式将上式再次改写
A
x
=
b
Ax=b
Ax=b
则可通过矩阵求逆的方式求出回归系数
x
=
A
−
1
b
x = A^{-1}b
x=A−1b
利用一阶多项式 Φ 1 = { 1 , x } Phi_1={1,x} Φ1={1,x}通过1.1节的方法对如下数据进行最小二乘拟合。
则法方程组为
[
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
拟合结果如下
#%% 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,2cdots,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=paralleldelta_iparallel ^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
Ain 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=FGquad (Fin mathbb{C} ^{m times r},Gin 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
Ain mathbb{C} ^{m times n},bin 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
本章算例与1.2节相同，以作对照。
利用一阶多项式
y
=
a
0
+
a
1
x
y=a_0+a_1x
y=a0+a1x对如下数据进行最小二乘拟合。
列写
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]
求解结果与第一章相同。
#%% 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,2cdots,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}=-2sum_{i=1}^{n}left[ y_i-(beta_0+beta_1x_i) right]&=0 \ Rightarrow qquad nbeta_0+beta_1sum_{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}=-2sum_{i=1}^{n}x_ileft[ y_i-(beta_0+beta_1x_i) right]&=0 \ Rightarrow qquad beta_0sum_{i=1}^nx_i+beta_1sum_{i=1}^nx_i^2 &=beta_1sum_{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)}{ nsum_{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 )}{ nsum_{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_1overline{x} \ hat{beta_1}&= frac{S_{xy}}{S_{xx}} end{aligned}
β0^β1^=y−β1x=SxxSxy
本章算例与1.2节相同，以作对照。
对如下数据进行线性回归。
计算系数
β
^
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_1overline{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
#%% 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函数的基本用法及其理解（傻瓜版）
#%% 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)
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）