什么是最小二乘估计
最小二乘方法是最常用的线性参数估计方法。这里的线性参数的意思并不是说我们估计出来的表达式的形式是线性的,即:只能估计线性函数的参数,如: y = a x + b y = ax+b y=ax+b,估计 a a a和 b b b,而是说我们应用一系列的观察值: ( x 0 , y 0 ) , ( x 1 , y 1 ) , . . . , ( x n − 1 , y n − 1 ) (x_0, y_0), (x_1, y_1), ..., (x_{n-1}, y_{n-1}) (x0,y0),(x1,y1),...,(xn−1,yn−1)来进行估计参数的过程中对观察值的操作是线性的操作,如: y = a 0 + a 1 x + . . . + a n − 1 x n − 1 y = a_0 + a_1x + ...+a_{n-1}x^{n-1} y=a0+a1x+...+an−1xn−1,通过对观察值的一系列线性组合,来估计参数 a 0 , a 1 , . . . , a n − 1 a_0, a_1, ..., a_{n-1} a0,a1,...,an−1。
有了这点概念之后,就来想怎么用这一系列点来估计参数呢?
最直观的一种就是,我们找到一组估计的参数,然后利用这一组估计的参数来得到估计的函数 f ^ ( x ) \hat{f}(x) f^(x),将 x x x代入得到一组估计的目标值,然后最理想的情况就是,我们估计出来的参数和原始真值相同,那么得到估计的目标值和真值之间的误差就是零,由于误差的存在,估计的目标值和真值之间存在一定的偏差,那么我们的目标就是选取一组估计的参数,使得其估计出来的目标值和真值之间的误差和最小。
引入目标函数
假设原真实的函数表达式为:
y
=
a
0
+
a
1
x
+
a
2
x
2
+
.
.
.
+
a
n
−
1
x
n
−
1
y = a_0+a_1x + a_2x^2 + ... + a_{n-1}x^{n-1}
y=a0+a1x+a2x2+...+an−1xn−1,此时我们观察到
m
m
m组值:
(
x
0
,
y
0
)
,
.
.
.
,
(
x
m
−
1
,
y
m
−
1
)
(x_0, y_0), ..., (x_{m-1}, y_{m-1})
(x0,y0),...,(xm−1,ym−1),此时我们就可以用一个矩阵方程来表示这些等式。
X
a
=
y
Xa = y
Xa=y,其中
X
=
(
1
x
0
x
0
2
.
.
.
x
0
n
−
1
1
x
1
x
1
2
.
.
.
x
1
n
−
1
.
.
.
1
x
m
−
1
x
m
−
1
2
.
.
.
x
m
−
1
n
−
1
)
X = \begin{pmatrix} 1 & x_0 & x_0^2 & ... & x0^{n-1} \\1 & x_1 & x_1^2 & ... & x_1^{n-1} \\ ... \\ 1 & x_{m-1} & x_{m-1}^2 & ... & x_{m-1}^{n-1} \end{pmatrix}
X=⎝
⎛11...1x0x1xm−1x02x12xm−12.........x0n−1x1n−1xm−1n−1⎠
⎞,
a
=
(
a
0
,
a
1
,
.
.
.
,
a
n
−
1
)
T
a = \begin{pmatrix} a_0, a_1, ..., a_{n-1}\end{pmatrix}^T
a=(a0,a1,...,an−1)T,
y
=
(
y
0
,
y
1
,
.
.
.
,
y
n
−
1
)
T
y = \begin{pmatrix} y_0, y_1, ..., y_{n-1}\end{pmatrix}^T
y=(y0,y1,...,yn−1)T其中需要估计的就是
a
a
a,一般情况下我们习惯将待估计的参数用
x
x
x代替,而用
A
A
A代表一个常数(已知)的矩阵,用
b
b
b表示方程的结果,于是上式就重写为:
A
x
=
b
Ax = b
Ax=b,其中
x
x
x就是我们代求的真值,而
A
A
A和
b
b
b是可以通过观察得到的已知量。
假设我们估计出一组值
x
^
=
(
a
0
^
,
a
1
^
,
.
.
.
,
a
n
−
1
^
)
T
\hat{x} = \begin{pmatrix} \hat{a_0}, \hat{a_1}, ..., \hat{a_{n-1}}\end{pmatrix}^T
x^=(a0^,a1^,...,an−1^)T,那么估计得到的值
b
^
=
A
x
^
\hat{b} = A\hat{x}
b^=Ax^,两者之间的差值就为
b
−
b
^
b-\hat{b}
b−b^,对其求二范数的平方可得:
(
b
−
b
^
)
T
×
(
b
−
b
^
)
(b-\hat{b})^T\times (b-\hat{b})
(b−b^)T×(b−b^),将其展开可得:
(
A
x
−
A
x
^
)
T
(
A
x
−
A
x
^
)
=
[
A
(
x
−
x
^
)
]
T
[
A
(
x
−
x
^
)
]
(Ax-A\hat{x})^T(Ax-A\hat{x}) = [A(x-\hat{x})]^T[A(x-\hat{x})]
(Ax−Ax^)T(Ax−Ax^)=[A(x−x^)]T[A(x−x^)]
将其进行展开,化简后为:
x
T
A
T
A
(
x
−
x
^
)
−
x
^
T
A
T
A
(
x
−
x
^
)
x^TA^TA(x-\hat{x}) - \hat{x}^TA^TA(x-\hat{x})
xTATA(x−x^)−x^TATA(x−x^)
由于我们需要的得到的是关于
x
^
\hat{x}
x^的极值点,所以将其对
x
^
\hat{x}
x^进行求导,并令其为0可得:
2
A
T
A
x
^
−
2
A
T
A
x
=
0
2A^TA\hat{x} - 2A^TAx = 0
2ATAx^−2ATAx=0
最终可得:
x
^
=
(
A
T
A
)
−
1
A
T
A
x
\hat{x} = (A^TA)^{-1}A^TAx
x^=(ATA)−1ATAx
由于
A
x
=
b
Ax = b
Ax=b,于是可得:
x
^
=
(
A
T
A
)
−
1
A
T
b
\hat{x} = (A^TA)^{-1}A^Tb
x^=(ATA)−1ATb
当
m
>
n
,
r
a
n
k
(
A
)
=
n
m>n, rank(A) = n
m>n,rank(A)=n,即
A
A
A为一个列满秩矩阵。
利用matlab进行仿真实现
说明:假设已知一个函数是一个三次函数,其表达式为: y = 1 + 1 x + 2 x 2 + 3 x 3 y = 1 + 1x + 2x^2 + 3x^3 y=1+1x+2x2+3x3,需要估计的参数就为: a 0 , a 1 , a 2 , a 3 a_0, a_1, a_2, a_3 a0,a1,a2,a3
x = 0:0.01:2;
x = x';
y_origin = 1 + 1*x + 2*x.^2 + 3*x.^3;
y = y_origin + 2*randn(length(x), 1);
A = [ones(length(y), 1) x x.^2 x.^3];
b = y;
param = inv((A'*A))*A'*b;
a0 = param(1);
a1 = param(2);
a2 = param(3);
a3 = param(4);
y_estimate = a0 + a1*x + a2*x.^2 + a3*x.^3;
figure;
plot(x, y);
hold on;
plot(x, y_origin);
hold on;
plot(x, y_estimate);
legend('观测数据', '标准数据', '估计数据');
得到的估计参数为:
大致的内容就是这样了,同时也是我暂时的理解,如果那里有问题,欢迎大家在评论中告知,如果有错,我会进行修改。