线性回归算法实现
规定
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}
f(x)=wTx+b,w=⎝⎛w1...wd⎠⎞
将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
D={(x1,y1),...,(xm,ym)},xi=⎝⎛xi1...xid⎠⎞,yi∈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}
X=⎝⎜⎜⎛x11x21...xm1x12x22...xm2............x1dx2d...xmd11...1⎠⎟⎟⎞
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-X@omiga).T @ (Y-X@omiga)
梯度迭代
# 代价列表
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()
1210
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
