基本思想
存在一组观察值
(
x
i
,
y
i
)
(x_i, y_i)
(xi,yi)，其中
y
i
y_i
yi和
x
i
x_i
xi之间满足一定的线性关系，如
y
=
a
0
f
0
(
x
)
+
a
1
f
1
(
x
)
+
.
.
.
+
a
m
−
1
f
m
−
1
(
x
)
y = a_0 f_0(x) + a_1 f_1(x) + ... + a_{m-1} f_{m-1}(x)
y=a0f0(x)+a1f1(x)+...+am−1fm−1(x)
其中
f
i
(
x
)
f_i(x)
fi(x)是已知的，
a
i
a_i
ai未知。
由于测量过程中存在一定的误差，所以导致测量得到的
y
i
y_i
yi具有一定的误差，即
y
i
=
a
0
f
0
(
x
i
)
+
a
1
f
1
(
x
i
)
+
.
.
.
+
a
m
−
1
f
m
−
1
(
x
i
)
+
e
i
y_i = a_0 f_0(x_i) + a_1 f_1(x_i) + ... + a_{m-1} f_{m-1}(x_i) + e_i
yi=a0f0(xi)+a1f1(xi)+...+am−1fm−1(xi)+ei
假设，我们由
n
n
n组观察值，将其写为矩阵的形式
y
⃗
=
H
a
⃗
+
e
⃗
\vec{y} = H\vec{a} + \vec{e}
y=Ha+e
y
⃗
\vec{y}
y和
e
⃗
\vec{e}
e为一个
n
×
1
n\times1
n×1的向量，
a
⃗
\vec{a}
a为
m
×
1
m\times 1
m×1的向量，
H
H
H为一个已知的
n
×
m
n\times m
n×m的矩阵
(
f
0
(
x
0
)
f
1
(
x
0
)
.
.
.
f
m
−
1
(
x
0
)
f
0
(
x
1
)
f
1
(
x
1
)
.
.
.
f
m
−
1
(
x
1
)
.
.
.
.
.
.
.
.
.
.
.
.
f
0
(
x
n
−
1
)
f
1
(
x
n
−
1
)
.
.
.
f
m
−
1
(
x
n
−
1
)
)
\begin{pmatrix} f_0(x_0) & f_1(x_0) & ... & f_{m-1}(x_0) \\ f_0(x_1) & f_1(x_1) & ... & f_{m-1}(x_1) \\ ... & ... & ... & ... \\ f_0(x_{n-1}) & f_1(x_{n-1}) & ... & f_{m-1}(x_{n-1})\end{pmatrix}
f0(x0)f0(x1)...f0(xn−1)f1(x0)f1(x1)...f1(xn−1)............fm−1(x0)fm−1(x1)...fm−1(xn−1)
由于测量误差满足高斯分布，由最大似然估计可得，选取使得
∣
y
⃗
−
H
a
⃗
∣
2
|\vec{y} - H\vec{a}|^2
∣y−Ha∣2达到最小值的
a
^
\hat{a}
a^就是最真实
a
⃗
\vec{a}
a的最优估计。
最小二乘估计表达式推导过程
令代价函数
L
(
a
⃗
)
=
∣
y
⃗
−
H
a
⃗
∣
2
L(\vec{a}) = |\vec{y} - H\vec{a}|^2
L(a)=∣y−Ha∣2，将
L
(
a
⃗
)
L(\vec{a})
L(a)对
a
⃗
\vec{a}
a进行求导可得
∂
L
(
a
⃗
)
∂
a
⃗
=
−
2
H
T
y
⃗
+
2
H
T
H
a
⃗
\frac{\partial{L(\vec{a})}}{\partial{\vec{a}}} = -2H^T\vec{y} + 2H^TH\vec{a}
∂a∂L(a)=−2HTy+2HTHa
令其为零可得
a
^
=
(
H
T
H
)
−
1
H
T
y
⃗
\hat{a} = (H^TH)^{-1}H^T\vec{y}
a^=(HTH)−1HTy
举例说明
利用matlab生成一组数据，满足如下的表达式
y
=
a
0
f
0
(
x
)
+
a
1
f
1
(
x
)
+
a
2
f
2
(
x
)
+
e
(
x
)
y = a_0 f_0(x) + a_1 f_1(x) + a_2 f_2(x)+e(x)
y=a0f0(x)+a1f1(x)+a2f2(x)+e(x)
其中
a
0
=
1
a_0 = 1
a0=1 ,
a
1
=
2
a_1 = 2
a1=2,
a
2
=
10
a_2 = 10
a2=10,
f
0
(
x
)
=
x
f_0(x) = x
f0(x)=x,
f
1
(
x
)
=
e
x
f_1(x) = e^x
f1(x)=ex,
f
2
(
x
)
=
x
2
f_2(x) = x^2
f2(x)=x2，
e
(
x
)
e(x)
e(x)服从标准正态分布
% 测试LSM算法
% f0 = 1; f1 = e^x; f2 = x^2
% a0 = 1; a1 = 2; a2 = 10
t = 0 : 0.01 : 3;
f0 = ones(1, length(t));
f1 = exp(t);
f2 = t.^2;
% f1 = t;
% f2 = t.^2;
a0 = 1;
a1 = 2;
a2 = 10;
H = [f0' f1' f2'];
y = a0*f0 + a1*f1 + a2*f2 + randn(1, length(t));
y = y';
est_a = inv(H'*H)*H'*y;
plot(y);
hold on
plot(est_a' * [f0;f1;f2]);
legend('带有噪声的观测量', '估计得到的观测量')
得到的结果如下：
得到的
a
^
=
[
0.9807
,
1.9565
,
10.0977
]
\hat{a} = [0.9807, 1.9565, 10.0977]
a^=[0.9807,1.9565,10.0977]，可以看到能够较好地估计原理论值
a
⃗
=
[
1
,
2
,
10
]
\vec{a} = [1, 2, 10]
a=[1,2,10]
