最小二乘法(Least Square，LS)
思想
-
最小二乘法目标：残差(观测点和估计点的差值)的平方和最小
-
用于在线性回归模型中估计未知参数，非线性需要迭代的方法
-
方程的数目要多于未知的参数数目
推导
? 变量y与x可控变量x1 x2 …xt之间关系：
y
=
β
1
x
+
β
2
x
+
.
.
.
+
β
t
x
+
ε
y=β_1x+β_2x+...+β_tx+ε
y=β1?x+β2?x+...+βt?x+ε
? 转成矩阵格式
V = A x ? B V=Ax-B V=Ax?B
? 式中,
KaTeX parse error: Undefined control sequence: \matrix at position 11: A=\left[ \?m?a?t?r?i?x?{ β_1^1 & β_2…
KaTeX parse error: Undefined control sequence: \matrix at position 11: B=\left[ \?m?a?t?r?i?x?{ y^1-ε^1\\ …
? 转化成最小二乘形式
m
i
n
∣
∣
V
∣
∣
min||V||
min∣∣V∣∣
? 求欧几里得空间的2范数作为距离
∣
∣
V
∣
∣
2
2
=
∣
∣
A
x
?
B
∣
∣
2
2
=
(
A
x
?
B
)
T
(
A
x
?
B
)
=
x
T
A
T
A
x
?
B
T
A
x
?
x
T
A
T
B
+
B
T
B
||V||_2^2= ||Ax-B||_2^2\\ =(Ax-B)^T(Ax-B)\\ =x^TA^TAx-B^TAx-x^TA^TB+B^TB
∣∣V∣∣22?=∣∣Ax?B∣∣22?=(Ax?B)T(Ax?B)=xTATAx?BTAx?xTATB+BTB
? 导数为0，求最小值
?
∣
∣
A
x
?
B
∣
∣
2
2
?
x
=
2
A
T
A
x
?
2
A
T
B
=
0
\frac{?||Ax-B||_2^2}{?x}=2A^TAx-2A^TB=0
?x?∣∣Ax?B∣∣22??=2ATAx?2ATB=0
? 最小二乘的解为：
x
=
(
A
T
A
)
?
1
(
A
T
B
)
x=(A^TA)^{-1}(A^TB)
x=(ATA)?1(ATB)
? 加权最小二乘(WLS)
x
估
计
=
(
A
T
P
A
)
?
1
(
A
T
P
B
)
x_{估计}=(A^TPA)^{-1}(A^TPB)
x估计?=(ATPA)?1(ATPB)
? 模型误差为(n为条件数,t为未知数,分别为矩阵A的行数的列数，n>t)
σ
=
(
∣
∣
A
x
估
计
?
B
∣
∣
2
2
n
?
t
)
σ=\sqrt(\frac{||Ax_{估计}-B||_2^2}{n-t})
σ=(?n?t∣∣Ax估计??B∣∣22??)
应用
直线拟合
- 仿真一组参数
X = range(0, 100, 10)
Y = [5 * _+ random.randint(0, 100) for _ in X]
- LS 拟合直线
y=kx+b
def cal_model(X, Y):
# y = k*x + b 最小二乘法
Bval, Lval = [], []
for _x, _y in zip(X, Y):
Bval.append([_x, 1])
Lval.append([_y])
Bmat = np.matrix(Bval)
Lmat = np.matrix(Lval)
gRes = (Bmat.T * Bmat).I * (Bmat.T * Lmat)
Vmat = Bmat * gRes - Lmat
sigma = np.sqrt(Vmat.T * Vmat / (len(X) - 2))
return gRes.tolist(), sigma
- 直线拟合误差
def cal_model_error(para, x, y):
y_gj = para[0][0] * x + para[1][0]
return y_gj - y
- 仿真结果可视化
x1 = 0
y1 = para[0][0] * x1 + para[1][0]
x2 = 90
y2 = para[0][0] * x2 + para[1][0]
fig = plt.figure()
plt.scatter(X, Y)
plt.plot([x1, x2], [y1, y2], color='r')
plt.rcParams['font.sans-serif'] = 'SimHei'
plt.title('LS直线拟合')
plt.show()
