線性迴歸的兩種演算法實現 梯度下降和解析法（西瓜書學習）
線性迴歸演算法實現
規定
d個屬性描述的示例
x
=
(
x
1
.
.
.
x
d
)
x= \begin{pmatrix} x_1\\ ...\\ x_d \end{pmatrix}
x=⎝⎛x1...xd⎠⎞
x
i
x_i
xi為第i個屬性的取值
即:
f
(
x
)
=
w
1
x
1
+
w
2
x
2
+
.
.
.
+
w
d
x
d
+
b
f(x)=w_1x_1+w_2x_2+...+w_dx_d+b
f(x)=w1x1+w2x2+...+wdxd+b
向量形式:
f
(
x
)
=
w
T
x
+
b
,
w
=
(
w
1
.
.
.
w
d
)
f(x)=w^Tx+b, w= \begin{pmatrix} w_1\\ ...\\ w_d \end{pmatrix}
將w和b一起寫為向量形式
即:
W
^
=
(
w
b
)
\hat{W}= \begin{pmatrix} w\\ b \end{pmatrix}
W^=(wb)
資料集表示為
D
=
{
(
x
1
,
y
1
)
,
.
.
.
,
(
x
m
,
y
m
)
}
,
x
i
=
(
x
i
1
.
.
.
x
i
d
)
,
y
i
∈
R
D=\{(x_1,y_1),...,(x_m,y_m)\}, x_i= \begin{pmatrix} x_{i1}\\ ...\\ x_{id} \end{pmatrix}, y_i\in R
資料集轉化為矩陣
X
=
(
x
11
x
12
.
.
.
x
1
d
1
x
21
x
22
.
.
.
x
2
d
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
x
m
1
x
m
2
.
.
.
x
m
d
1
)
X= \begin{pmatrix} x_{11} & x_{12} & ... & x_{1d} & 1\\ x_{21} & x_{22} & ... & x_{2d} & 1\\ ... & ... & ... & ... & ...\\ x_{m1} & x_{m2} & ... & x_{md} & 1 \end{pmatrix}
Y
=
(
y
1
.
.
.
y
m
)
Y= \begin{pmatrix} y_1\\ ...\\ y_m \end{pmatrix}
Y=⎝⎛y1...ym⎠⎞
代價函式
J
(
W
^
)
=
1
2
m
(
Y
−
X
W
^
)
T
(
Y
−
X
W
^
)
J(\hat{W})=\frac{1}{2m}(Y-X\hat{W})^T(Y-X\hat{W})
J(W^)=2m1(Y−XW^)T(Y−XW^)
引數估計和解釋
(
w
,
b
)
=
arg
min
(
w
,
b
)
J
(
W
^
)
(w,b)=\mathop{\arg\min}\limits_{(w,b)}J(\hat{W})
(w,b)=(w,b)argminJ(W^)
使得
f
(
x
)
=
X
W
^
≃
Y
f(x)=X\hat{W}\simeq Y
f(x)=XW^≃Y
均方誤差有非常好的幾何意義，它對應了常用的歐式距離。基於均方誤差最小化進行模型求解的方法加最小二乘法。線上性迴歸中，最小二乘法就是試圖找到一個條直線，使得所有樣本到直線上的距離最小
梯度下降
J
(
W
^
)
J(\hat{W})
J(W^)導數為
1
m
X
T
(
X
W
^
−
Y
)
\frac{1}{m}X^T(X\hat{W}-Y)
m1XT(XW^−Y)
W
^
−
=
J
′
(
W
^
)
\hat{W}-=J^{'}({\hat{W}})
W^−=J′(W^)
正規方程解法
公式為
W
^
=
(
X
T
X
)
−
1
X
T
Y
=
X
−
1
Y
\hat{W}=(X^TX)^{-1}X^TY =X^{-1}Y
W^=(XTX)−1XTY=X−1Y
匯入庫
import numpy as np
import matplotlib.pyplot as plt
準備資料
# 樣本數量
m=100
# 初始化資料
X=np.linspace(0,10,m)
Y=X+2+np.random.randn(m)
# 展示樣本資料
plt.scatter(X,Y)
plt.show()
# 將樣本資料矩陣化
X=X.reshape((m,1))
X1=np.hstack((np.ones_like(X),X))
Y=Y.reshape((m,1))
# 引數
omiga=np.zeros(2).reshape((2,1))
# 步長
alpha=0.01
precision=1e-3
代價函式
# 代價函式
def cost(omiga,X,Y):
return (Y-[email protected]).T @ (Y-[email protected])
梯度迭代
# 代價列表
costs=[]
while True:
costs.append(cost(omiga,X1,Y)[0][0])
omiga -=alpha /m * X1.T @ ( X1 @ omiga - Y )
if len(costs)>1:
if costs[-2]-costs[-1]<=precision:
break
# 代價函式曲線
plt.plot(np.arange(len(costs)),costs)
plt.show()
展示結果
plt.plot(X,X1 @ omiga)
plt.scatter(X,Y)
plt.show()
print(omiga)
[[1.80041713]
[1.05535637]]
正規方程
# omiga = np.linalg.pinv(X1.T @ X1) @ X1.T @ Y
# print(omiga)
omiga=np.linalg.pinv(X1) @ Y
print(omiga)
[[1.88910625]
[1.04202106]]
plt.plot(X,X1 @ omiga)
plt.scatter(X,Y)
plt.show()
