数据拟合

来源：互联网发布：淘宝照片大小是多少编辑：程序博客网时间：2024/04/30 08:41

假设对观测数据进行拟合，得到的拟合曲线为 $\hat y%$ 。将观测数据 ${x_i}%$ 代入 $\hat y%$ ，得到 $\hat y\left( {{x_i}} \right)%$ ，其和 ${y_i}%$ 的偏差定义为

${d_i} = g\left( {{y_i} - \hat y\left( {{x_i}} \right)} \right)%$ （1）

评价 $\hat y%$ 拟合结果好坏的函数称为指标函数

$I\left( {\hat y} \right) = h\left( {{d_1},{d_2}, \cdots ,{d_n}} \right)%$ （2）

拟合函数 $\hat y%$ 在观测数据上总的偏差越小，说明拟合的越好，因此 $I\left( {\hat y} \right)%$ 可以具体写成

$I\left( {\hat y} \right) = h\left( {{d_1},{d_2}, \cdots ,{d_n}} \right) = \sum\limits_{i = 1}^n {{\omega _i}{d_i}}%$ （3）

其中 $\sum\limits_{i = 1}^n {{\omega _i}} = 1%$ ， ${\omega _i}%$ 是加权系数，表示每个观测数据的重要程度。一般情况下，我们认为每个观测数据是等重要的。因此，上式可以简化为

$I\left( {\hat y} \right) = \sum\limits_{i = 1}^n {{d_i}}%$ （4）

其中偏差可以定义为

$g\left( \cdot \right) = {\left| \cdot \right|_p}%$ （5）

当p=0时， ${d_i} = g\left( {{y_i} - \hat y\left( {{x_i}} \right)} \right) = {\left| {{y_i} - \hat y\left( {{x_i}} \right)} \right|^0}%$

当p=1时， ${d_i} = g\left( {{y_i} - \hat y\left( {{x_i}} \right)} \right) = {\left| {{y_i} - \hat y\left( {{x_i}} \right)} \right|}%$

当p=2时， ${d_i} = g\left( {{y_i} - \hat y\left( {{x_i}} \right)} \right) = {\left| {{y_i} - \hat y\left( {{x_i}} \right)} \right|^2}%$

。。。

评价函数 $I\left( {\hat y} \right)%$ 一般写成p范数的p次方

$I\left( {\hat y} \right) = \sum\limits_{i = 1}^n {{d_i}} = \sum\limits_{i = 1}^n {g\left( {{y_i} - \hat y\left( {{x_i}} \right)} \right)} = \sum\limits_{i = 1}^n {{{\left| {{y_i} - \hat y\left( {{x_i}} \right)} \right|}_p}} = \left\| {y - \hat y\left( x \right)} \right\|_p^p%$ （6）

求最优拟合函数 $\hat y%$ 的过程是在 $\hat y%$ 和 $I\left( {\hat y} \right)%$ 构成的空间上寻优的过程

$\mathop {\min }\limits_{\hat y} I\left( {\hat y} \right) = \mathop {\min }\limits_{\hat y} \left\| {y - \hat y\left( x \right)} \right\|_p^p%$ （7）

对应的解是

$\hat y^*} = \mathop {\arg \min }\limits_{\hat y} I\left( {\hat y} \right) = \mathop {\arg \min }\limits_{\hat y} \left\| {y - \hat y\left( x \right)} \right\|_p^p%$ （8）

指标函数是一个以函数为自变量的函数。至此，这个问题是变分问题（至于说如何使用变分来求解，暂且压一压，以后再细说）。

假设拟合函数为 $\hat y = X\beta%$ 。 $Y=\begin{bmatrix}y_{1}\\ y_{2}\\ \vdots\\ y_{n}\end{bmatrix}%$ ， $X=\begin{bmatrix} 1 \, x_{1}\\ 1 \, x_{2}\\ \vdots \, \vdots\\ 1 \, x_{n} \end{bmatrix}%$ ， $\beta=\begin{bmatrix} \beta_{1}\\ \beta_{2}\\ \end{bmatrix}%$ 。

如何求出 $\beta%$ 呢？，从p=2开始推导。

$\min I\left( \beta \right) = \min \left\| {y - \hat y\left( x \right)} \right\|_2^2 = \min \left\| {Y - X\beta } \right\|_2^2%$

评价函数对 $\beta%$ 求导（因为在 $\beta%$ 的取值范围是连续的，且向量的 $\left\| \cdot \right\|_2^2%$ 是一个凸函数，至少是不凹的），

$\frac{d}{{d\beta }}I\left( \beta \right) = \frac{d}{{d\beta }}\left\| {Y - X\beta } \right\|_2^2 = - 2{\left( {Y - X\beta } \right)^T}X%$ （9）

$\frac{d}{{d{\beta ^T}}}\left( {\frac{d}{{d\beta }}I\left( \beta \right)} \right) = \frac{d}{{d{\beta ^T}}}\left( { - 2{{\left( {Y - X\beta } \right)}^T}X} \right) = 2{X^T}X%$ （10）

当 $\frac{d}{{d\beta }}I\left( \beta \right) = 0%$ 时， $\frac{d}{{d{\beta ^T}}}\left( {\frac{d}{{d\beta }}I\left( \beta \right)} \right) \ge 0%$ ，因此极值点为极小值点。

$\frac{d}{{d\beta }}I\left( \beta \right) = 0%$ 等价于

${\left( {Y - X\beta } \right)^T}X = 0%$ （11）

上式两边进行转置，得

${X^T}Y = {X^T}X\beta%$ （12）

X为列满秩矩阵， ${X^T}X%$ 的逆矩阵存在。因此

$\hat \beta = {\left( {{X^T}X} \right)^{ - 1}}{X^T}Y%$ （13）

$\hat Y = X\hat \beta = X{\left( {{X^T}X} \right)^{ - 1}}{X^T}Y%$ （14）

当p=0时， $\min I\left( \beta \right) = \min {\left\| {y - \hat y\left( x \right)} \right\|_0} = \min {\left\| {Y - X\beta } \right\|_0}%$ 。该式表示，拟合直线穿过观测数据点越多越好。该方法对数据的分布以及其中掺杂的噪声比较敏感，解不稳定。该方法可能不能拟合出真实数据的曲线，而是拟合了噪声数据。

当p=1时， $\min I\left( \beta \right) = \min {\left\| {y - \hat y\left( x \right)} \right\|_1} = \min {\left\| {Y - X\beta } \right\|_1}%$ 。评价函数导数不连续，求解不像二范数求解那么方便，因此很少使用，但在有些情况下，1范数拟合的曲线更好。

下面用一段代码，简单说明一下p范数对拟合结果的影响

[plain] view plaincopy
<span style="font-size:14px;">lens=50;  
b=50;  
   
x=1:lens;  
y=2*x+b*randn(1,lens);  
   
x=[1 2 3 4 5];  
y=[1 2 2 3 5];  
   
for k=1:200  
    for b=1:200  
        yk=(k-100)/10*x+(b-100)/10;  
        d=y-yk;  
        sign=ones(1,lens);  
        sign(find(abs(d)<0.5))=0;  
        err0(k,b)=sum(sign);  
        err1(k,b)=sum(abs(d));  
        err2(k,b)=sum(d.*d);  
    end  
end  
figure(1)  
mesh(err0)  
figure(2)  
mesh(err1)  
figure(3)  
mesh(err2)  
   
figure(4)  
plot(x,y,'ro')  
hold on  
   
[a1 a2]=min(err0);  
[b1 b2]=min(a1);  
ka0=(a2(b2)-100)/10;  
ba0=(b2-100)/10;  
yk0=ka0*x+ba0;  
plot(x,yk0,'md')  
   
   
[a1 a2]=min(err1);  
[b1 b2]=min(a1);  
ka1=(a2(b2)-100)/10;  
ba1=(b2-100)/10;  
yk1=ka1*x+ba1;  
plot(x,yk1,'g+')  
   
   
[a1 a2]=min(err2);  
[b1 b2]=min(a1);  
ka2=(a2(b2)-100)/10;  
ba2=(b2-100)/10;  
yk2=ka2*x+ba2;  
plot(x,yk2,'bs')  
   
   
e10=sum(abs(y-yk0))  
e11=sum(abs(y-yk1))  
e12=sum(abs(y-yk2))  
   
e20=sum((y-yk0).*(y-yk0))  
e21=sum((y-yk1).*(y-yk1))  
e22=sum((y-yk2).*(y-yk2))  
   
legend('data','L0','L1','L2')  
hold off</span>  

这段代码首先生成一组数据（x，y），然后分别使用0、1和2范数进行拟合求解。搜索范围k=[-10,10]，b=[-10,10]。在搜索空间中找到拟合误差最小的最小p乘解。

从评价函数的2范数出发， $\beta%$ 的解的形式中包含着复杂的矩阵运算关系，这其中应该蕴含着什么。让我们首先从一个简单的例子入手吧。

有a，b两个向量，b在a上的投影p可以写成投影长度x和a方向上单位向量的乘积

$p = \frac{a}{{\sqrt {{a^T}a} }}x = \frac{a}{{\sqrt {{a^T}a} }}\frac{{{a^T}b}}{{\sqrt {{a^T}a} }} = \frac{{a{a^T}}}{{{a^T}a}}b%$ （15）

往a方向上做投影的投影变换矩阵为

$P=\frac{{a{a^T}}}{{{a^T}a}}%$ （16）

$p=Pb%$ （17）

向量b到向量a的距离等于e=b-p的模，垂直方向上的投影矩阵为

${P_ \bot } = I - P%$ （18）

$e = b - p = b - Pb = \left( {I - P} \right)b = {P_ \bot }b%$ （19）

在数据拟合中，数据的个数远远多于未知数（待求解参数）的个数，因此这个方程不能得到精确解。我们需要找到距离Y最近的一个空间，将Y投影到该空间中。X的列向量构成的空间叫做列空间。该空间内的任意向量Xv都是X的列向量的线性组合得到，v就是组合系数。Y的近似解 $\hat Y%$ 形如

$\hat Y = X\hat \beta%$ （20）

近似解 $\hat Y%$ 与Y之间的向量为

$e = Y - \hat Y%$ （21）

该向量垂直于X的列空间中的任意向量Xv，有

${\left( {Xv} \right)^T}e = 0%$ （22）

${v^T}{X^T}\left( {Y - \hat Y} \right) = 0%$ （23）

因为v为任意向量，因此

${X^T}\left( {Y - \hat Y} \right) = 0%$ （24）

（20）式两边同时左乘矩阵X的转置，将（24）式代入得到

${X^T}Y = {X^T}X\hat \beta%$ （25）

解出

$\hat \beta = {\left( {{X^T}X} \right)^{ - 1}}{X^T}Y%$ （26）

近似解

$\hat Y = X\hat \beta = X{\left( {{X^T}X} \right)^{ - 1}}{X^T}Y%$ （27）

Y在投影变换矩阵 $P=X{\left({{X^T}X} \right)^{-1}}{X^T}%$ 的作用下，投影到X的列空间上，投影后得到的向量为 $\hat Y=PY%$ ，误差向量 $e = Y - \hat Y = \left( {I - P} \right)Y%$ ，误差的大小为 ${\left\| {Y - \hat Y} \right\|_p}%$ 。

投影矩阵从另外一个角度解释了最小二乘法，同时也是最小p（p>2）乘法的解释。我们发现，最小二乘法给出的解是近似解。误差的来源方方面面，比如系统偏差，观测误差，记录误差等等。这些因素之间是一个什么样的关系，线性的还是非线性的，一时半会儿说不清楚，我们偷个懒，将拟合问题重新形式化如下

$Y = X\hat \beta + \varepsilon%$ （28）

我们要求的问题等价于

$\hat \beta = \mathop {\arg \max }\limits_\beta p\left( {\left. Y \right|X,\beta } \right)%$ （29）

搞工程的搞来搞去不经意的发现，当随机不可知的因素很多，独立随机试验的次数很大时，由这些随机因素造成的随机误差服从高斯分布。学术界的人按耐不住了，不能让搞工程的压下去，整出了一个中心极限定理，各种分布都会渐进服从高斯分布。

我们俗气一把，还是从高斯分布入手，假设误差服从零均值高斯分布，即 $\varepsilon \sim N\left( {0,{\sigma ^2}} \right)%$ 。为什么是零均值，因为我们不想估计出个没用的有偏的分布出来。这种零均值误差还有个很带感的名字，白噪声。对其进行傅里叶变换后，各个频率都有响应，也就是说这种噪声是由不同频率的噪声合成的。我们常见的由各个频率合成的事物就是白色的光。所以，按照国际惯例，这种噪声叫做白噪声。

另外一个假设是独立同分布，也就是每次实验都是独立的，但是服从的分布相同。

$p\left( \varepsilon \right) = \frac{1}{{\sqrt {2\pi {\sigma ^2}} }}{e^{ - \frac{1}{2}{{\left( {\frac{\varepsilon }{\sigma }} \right)}^2}}}%$ （30）

$p\left( {\left. Y \right|X,\beta } \right) = \frac{1}{{\sqrt {2\pi {\sigma ^2}} }}{e^{ - \frac{1}{2}{{\left( {\frac{{Y - X\beta }}{\sigma }} \right)}^2}}}%$ （31）

令上式对 $\beta%$ 的导数为零，又见 ${X^T}Y = {X^T}X\hat \beta%$ 。也就是说，最小二乘法的概率解释是，拟合误差服从零均值高斯分布，拟合直线通过均值点。

假设噪声服从形如 $p\left( \varepsilon \right) = {e^{ - \left\| \varepsilon \right\|_p^p}}%$ ，那么其解对应的就是最小p范数解。这种分布貌似就是指数族分布，当p=1时叫做拉普拉斯分布，p=2时就是大名鼎鼎的高斯分布。

通过上面的分析，我们可以察觉到，噪声的模型影响了最终的拟合结果。如果噪声是形如1范数的，那么用2范数的最小二乘法拟合出来的直线就存在偏差。如果噪声中有粗大误差，那么如果不能事先去除，拟合的结果很有可能拟合了噪声而没能拟合真实数据。从这个角度看，其实我们能够解决的问题还比较有限，因为我们对于噪声的认识还不够。

从p范数的角度出发，均值是误差形如2范数的解，中位数是误差形如1范数的解，众数是误差形如0范数的解。猜想，p阶矩就是误差形如p范数的解。

到这里，我丧心病狂的把凸优化，线性代数和概率论貌似完美的在数据拟合的框架下联系在了一起。从其反方向看，这些理论本来就是为了解决数据拟合而被提出来的。只是数学教学时，人为的把本来应该在一起的拆开了。

图中的数据在y=2x上加入了标准差为10的零均值高斯白噪声 $\varepsilon \sim N\left(0,100 \right)%$ 。使用上面推导的公式解出拟合曲线为

拟合误差的样本均值为0.abc%#%$^#e-13，样本标准差为9.44。同时，该拟合曲线通过样本的均值点。

[plain] view plaincopy
<span style="font-size:14px;">a=4;  
b=10;  
lens=100;  
   
x=1:lens;  
y=zeros(1,lens);  
y=2*x+b*randn(1,lens);  
   
X=[x;y]';  
   
Y0=zeros(lens,2);  
Y1=zeros(lens,2);  
Y2=zeros(lens,2);  
Y3=zeros(lens,2);  
pY=zeros(lens+10,2);  
   
A=zeros(lens,4);  
V=zeros(lens,4);  
D=zeros(lens,4);  
pY1=zeros(lens+10,2);  
   
figure  
for i=1:lens  
    Y0(i,:)=X(i,:);  
         
    if i>3 && a>=4  
       Y3(i,:)=Y0(i,:)-3*Y0(i-1,:)+3*Y0(i-2,:)-Y0(i-3,:);  
    end  
    if i>2 && a>=3  
        Y2(i,:)=Y0(i,:)-2*Y0(i-1,:)+Y0(i-2,:);  
    end  
    if i>1 && a>=2  
        Y1(i,:)=Y0(i,:)-Y0(i-1,:);  
    end  
     
    pY(i+1,:)=Y0(i,:)+Y1(i,:)+0.5*Y2(i,:)+1/6*Y3(i,:);  
     
    if i<lens  
        t=X(i+1,:)'*pinv(X(i,:)');  
        pY1(i+2,:)=(t*X(i+1,:)')';  
        [v d]=eig(t(1:2,1:2));  
        A(i,:)=reshape(t,1,4);  
        V(i,:)=reshape(v,1,4);  
        D(i,:)=reshape(d,1,4);  
    end  
end  
   
X1=[ones(1,lens);x]';  
beta=inv(X1'*X1)*X1'*y';  
pY2=X1*beta;  
   
plot(x,y,'go')  
hold on  
plot(pY(a:lens,1),pY(a:lens,2),'r+');  
plot(pY1(a:lens,1),pY1(a:lens,2),'bs');  
plot(x,pY2,'mv');</span>  

上面的代码中pY是基于泰勒级数展开的近似估计，pY1是基于局部最小二乘估计，pY2是全局最小二乘估计。

说了这么多，我们原地踏步在一阶线性估计上，高阶怎么求解。改改X和 $\beta%$ 就可以了。求解方法还是原来的配方，还是原来的味道。这里就不在多啰嗦了。

$X=\begin{bmatrix} 1&{x_1^1}& \cdots &{x_1^m}\\ 1&{x_2^1}& \cdots &{x_2^m}\\ \vdots & \vdots & \ddots & \vdots \\ 1&{x_n^1}& \cdots &{x_n^m} \end{bmatrix}%$ （32）

$\beta=\begin{bmatrix} {\beta_0}\\ {\beta_1}\\ \vdots\\ {\beta_m} \end{bmatrix}%$ （33）

给大家再拜上一记大杀器

$X=\begin{bmatrix} {h_0}\left( 1 \right)&{h_1}\left( {{x_1}} \right)& \cdots &{h_m}\left( {{x_1}} \right)\\ {h_0}\left( 1 \right)&{h_1}\left( {{x_2}} \right)& \cdots &{h_m}\left( {{x_2}} \right)\\ \vdots & \vdots & \ddots & \vdots \\ {h_0}\left( 1 \right)&{h_1}\left( {{x_n}} \right)& \cdots &{h_m}\left( {{x_n}} \right) \end{bmatrix}%$ （34）

小伙伴们不要害怕，看 $h\left( \cdot \right)%$ 一副卖萌的样子就知道这货是传说中的核函数。加权函数一块来的时候，更加凶残的公式也就来了

$\hat \beta = {\left( {{X^T}WX} \right)^{ - 1}}{X^T}WY%$ （35）

从数据出发，我们总能对它们进行合适的解释，发现合适的模型，拟合出合适的曲线，并给出这种解释的好坏程度。到目前为止，我们对于数据拟合能够给出的最简洁的表达是

${\hat y^*} = \mathop {\arg \min }\limits_{\hat y} I\left( {\hat y} \right) = \mathop {\arg \min }\limits_{\hat y} \sum\limits_{i = 1}^n {{\omega _i}\left\| {{y_i} - {x_i}\beta } \right\|_p^p}%$ （36）

貌似我们可以满足的洗洗睡了，但是一个关键的问题是X的具体形式能不能由数据自己说出来，也就是说从无到有可不可能。此外，还有两个问题没有解决，如何在保证上式成立的同时将参数的空间最小化，也就说能够用线性拟合的不用二阶多项式拟合。第二个问题就是迭代拟合。现实情况中，我们接收到的数据往往是以时间序列的形式呈现的，那t时刻和t-1时刻的模型出现不一致时，怎么处理。

除了最小p乘，还有PCA

如图1所示，最小p乘法求得是 $\min \sum\limits_{i = 1}^n {\left\| {{d_i}} \right\|_p^p}%$ ，而真实值到拟合曲线的距离为 ${d'_i}%$ 。那么， $\min \sum\limits_{i = 1}^n {\left\| {{d'_i}} \right\|_p^p}%$ 对应的是什么样的数据分析呢？

图1 最小p乘法的使用的误差是 ${d_i}%$ 。真实值到拟合曲线的距离为 ${d'_i}%$

假如存在拟合曲线，设直线方程为 $y = kx + b%$ 。真实值 ${X_i} = {\left( {{x_i},{y_i}} \right)^T}%$ 到该曲线的投影点为 ${X'_i} = {\left( {{x'_i},{y'_i}} \right)^T}%$ 。p=2时，则两点之间的距离为

$\left\| {{{d'}_i}} \right\|_2^2 = \left\| {{X_i} - {{X'}_i}} \right\|_2^2%$ （37）

$I = \sum\limits_{i = 1}^n {\left\| {{X_i} - {{X'}_i}} \right\|_2^2}%$ （38）

点 ${X'_i} = {\left( {{x'_i},{y'_i}} \right)^T}%$ 在直线上，同时 $\left( {{X_i} - {{X'}_i}} \right) \bot \left( {1,k} \right)%$ 。这两个条件构成如下方程组

$\left\{\begin{matrix} {\left( {{x_i} - {{x'}_i},{y_i} - {{y'}_i}} \right) \cdot \left( {1,k} \right) = 0}\\ {y' = kx' + b} \end{matrix}\right.%$ （39）

联立上述方程组求得

$\left\| {{{d'}_i}} \right\|_2^2 = \frac{1}{{1 + {k^2}}}{\left( {{y_i} - k{x_i} - b} \right)^2}%$ （40）

代入式（37）（38）得

$I = \sum\limits_{i = 1}^n {\frac{1}{{1 + {k^2}}}{{\left( {{y_i} - k{x_i} - b} \right)}^2}}%$ （41）

上式两边对b求偏导，令偏导数为零得

$\frac{{\partial I}}{{\partial b}} = \sum\limits_{i = 1}^n {\frac{{ - 2}}{{1 + {k^2}}}\left( {{y_i} - k{x_i} - b} \right)} = 0%$ （42）

化简为

$\sum\limits_{i = 1}^n {\left( {{y_i} - k{x_i} - b} \right)} = 0%$ （43）

$\sum\limits_{i = 1}^n {{y_i}} - k\sum\limits_{i = 1}^n {{x_i}} - nb = 0%$ （44）

令

${m_x} = \frac{1}{n}\sum\limits_{i = 1}^n {{x_i}}%$ （45）

${m_y} = \frac{1}{n}\sum\limits_{i = 1}^n {{y_i}}%$ （46）

将式（45）（46）代入式（44）得

$n{m_y} - nk{m_x} - nb = 0%$ （47）

$b = {m_y} - k{m_x}%$ （48）

因此，该直线通过均值点 $m = {\left( {{m_x},{m_y}} \right)^T}%$ ，投影点可以改写为

${X'_i} = m + {e^T}({X_i} - m)e%$ （49）

其中e是直线方向的单位向量。将式（49）代入式（38）得

$I = \sum\limits_{i = 1}^n {\left\| {m + {e^T}({X_i} - m)e - {X_i}} \right\|_2^2}%$

$I = \sum\limits_{i = 1}^n {\left\| {{e^T}({X_i} - m)e} \right\|_2^2} + \sum\limits_{i = 1}^n { - 2{e^T}({X_i} - m){{({X_i} - m)}^T}e} + \sum\limits_{i = 1}^n {\left\| {({X_i} - m)} \right\|_2^2}%$

$I = - \sum\limits_{i = 1}^n {{e^T}({X_i} - m){{({X_i} - m)}^T}e} + \sum\limits_{i = 1}^n {\left\| {({X_i} - m)} \right\|_2^2}%$

化简为

$I = - {e^T}Se + \sum\limits_{i = 1}^n {\left\| {({X_i} - m)} \right\|_2^2}%$ （50）

其中 $S = ({X_i} - m){({X_i} - m)^T}%$ 。上式中等号右边的第二项是个常量，不影响I取得极值时对应的e，可以去掉。同时，我们假设e是单位向量，则 $\left\| e \right\|_2^2 = 1%$ 。重写I如下

$u = {e^T}Se + \lambda \left( {1 - {e^T}e} \right)%$ （51）

上式两边对e求导得

$\frac{{\partial u}}{{\partial e}} = 2Se - 2\lambda e = 0%$ （52）

化简得

$Se = \lambda e%$ （53）

上式成立时，u取得最大值，I取得最小值。对上两边同时除以（n-1），得到数据矩阵的协方差矩阵。 $\lambda%$ /(n-1)是协方差矩阵的特征值，e是对应的特征向量。上述推导过程可以较为简单的推广到m维空间。对特征值按降序排列， ${\lambda _1} \ge {\lambda _2} \ge ... \ge {\lambda _m}%$ ，其中m为数据变量的维度。 ${\lambda_1}%$ 对应着数据的主方向。经过特征向量矩阵的映射，将协方差矩阵投影为对角阵，变量之间的相关性被消除，而数据方差最大的方向就是主方向。

当计算出数据协方差矩阵的特征向量后，我们计算贡献率

$pareto = \frac{{{\lambda _1} + {\lambda _2} + ... + {\lambda _k}}}{{{\lambda _1} + {\lambda _2} + ... + {\lambda _m}}} \ge 1 - \alpha%$ （54）

求出

${k^*} = \mathop {\arg pareto}\limits_k \left( k \right) \ge 1 - \alpha%$ （55）

使用前 ${k^*}%$ 个特征值和特征向量压缩原来的数据的表达空间，同时还能保证压缩后的数据矩阵损失最小。上述方法就是我们熟悉的PCA。

主方向线通过数据矩阵的均值点，这个点对应的是使用PCA做人脸识别时求出的均值脸。

总结一下PCA的推导过程，

1、去数据变量样本间的均值，并将该均值从数据矩阵中减去，得到零均值矩阵。

2、求零均值数据矩阵的协方差矩阵。

3、求协方差矩阵的特征向量和特征值。

4、按照一定的比例选择特征值和特征向量，实现降维。

上面推导的是线性关系的PCA，对于非线性的数据上面的方法可能会失效。解决方法，使用核函数将数据映射到高维再进行上述分析，这貌似就是kernel PCA。

同时，PCA分析的主方向通过数据的均值。而数据的均值很采样数据紧相关，如果数据中存在粗大误差，那么此时的均值不能反应真实的数据均值。如果不进行预处理，后续的PCA分析很可能会是错误的，达不到预期的效果。因此，对数据进行预处理是很必要的，剔除粗大误差后再进行PCA分析，貌似就是robust PCA。

在PCA推导的过程中，我们可以较为清晰地看到，如果将数据标签揉到推导中，修改优化的目标函数，我们应该能推导出SVM。因此，不同的误差定义，不同的优化目标函数推导出了不同的数据分析方法。无论这些怎么变换花样，其依托的数学思想都是一致的。

说到这里，我们差不多吧数据拟合相关的数据分析方法说了遍，但是说来说去关键问题还是没有触及，我们最关心最希望自动化的东西没有设计，那就是数据的模式，线性的还是非线性的，一阶的还是二阶的等等问题。因为，我们明明可以看到数据在图像上或者几何上呈现出了某种分布，但是却不能通过数学推导自动化的把它从数据中挖出来。

到底能不能是一个终极问题，这个系列的文章只能做个抛砖引玉，希望能够激发出大家的进一步迭代思考，也许这种模式对应的数学公式就在不远处。

Kalman滤波

数据拟合能够估计出数据变化的趋势，另外一个同等重要的应用是如何利用这一趋势，预测下一时刻数据可能的值。通俗点儿说，你观察苍蝇（蚊子，蜜蜂）飞了几秒，你也许会想“它下一个时刻可能在哪儿”，“呈现出什么样的状态”诸如此类的问题。预知未来这档子事儿对我们有一种不可抗拒的吸引力。别看我们预测的未来很近，但这对于实际应用有很大的帮助。比如减小解空间的范围，便于搜索。对于搜索问题，预测可以看成是对从当前状态到目标状态的启发评价函数。好吧，我承认我陷得太深了，都是复习人工智能搞得。扯得有点儿远了，继续说我们的主题，预测。

古人每遇到重大活动，都会卜上一卦。念几句咒语，抽个签，看看签释，心里大概对所问之事有了个谱儿。再比如，这几天你的左眼皮一直在跳，你想知道这是为什么，意味着什么。你跑去算了一卦。抽签的时候，你心里默念着是不是要捡到钱了等等，结果抽了一个上上签，说你要遇到好事儿。“这几天眼皮跳”是你的观察数据。“你想知道未来会发生什么”是我们想要预测的东西。抽签的时候你心里默念的话，签儿，签上的符号和某些事件的对应关系，这些都是预测的算法。虽然占卜的过程包含了观察，有预测算法，有预测结果，同时也有结果的方差范围等等。但是我们说这种预测是不科学的，因为预测算法不科学，因果关系不见得成立等等。那有没有科学的预测呢，让我们进入今天的话题，Kalman滤波。

假设这样一个场景，A先生使用遥控器控制一架四轴直升飞机F在一个空旷的场地上飞行。直升飞机F上有一个GPS模块，通过无线发射模块实时的将直升飞机F的位置发给计算终端C。B先生在终端C上运行一个“打”直升飞机F的程序D。程序D根据终端C接收的GPS数据，指导一个虚拟的导弹E去跟踪直升飞机F，并试图将F“击落”。

A先生控制的直升飞机F飞行轨迹多变，很难被跟踪。同时，终端接收的GPS数据中还有噪声。B先生引导的导弹E燃料有限，因此不能长时间、频繁地机动。因此，B先生希望程序D要尽可能准的估计出F的位置，尽可能少机动，跟踪F并将其击落。

假如，A先生控制F急速攀升，如图1a。D得到的数据如图1b中的红色点。如果不进行预测，直接根据GPS数据控制E机动，E的运动轨迹如图1c，绿色的轨迹线。E很可能因燃料不足提前爆炸，而没有击中F。B先生很希望D能够根据GPS数据计算出如图1d所示的轨迹（橙色的轨迹线），来引导E去追踪F。

a b c d

图1 a、被观测对象实际的运动轨迹；b、我们观测到的被观测对象的运动轨迹；

c、如果不滤波的话，预测的轨迹；d、滤波后的预测轨迹。

假如A先生很狡猾，他控制F飞行，其飞行轨迹如图2所示。B先生深感压力巨大，如何才能有效的跟踪F，并将其击落呢。

图2 A先生控制的直升飞机F可能飞出的轨迹（终端C得到的GPS数据）。

颜色越深，获取数据的时间越早；反之，颜色越浅，获取时间越晚。

程序D根据终端C提供的GPS数据，估计F的位置（x，y），时间t的采样记为 ${s_t} = \left( {{x_t},{y_t}} \right)%$ 。使用前t个时刻的采样 ${s_1},{s_2},...,{s_t}%$ ，估计 $\hat f = g\left( {{s_1},{s_2},...,{s_t}} \right)%$ ，使得该估计满足，

$\min \sum\limits_{i = 1}^t {\left\| {\hat f\left( i \right) - {s_i}} \right\|} _p^p%$ (1)

$\min \left\| {\hat f\left( {i + 1} \right) - {s_{i + 1}}} \right\|_p^p%$ (2)

其中，为 ${\left\| \cdot \right\|_p}%$ 范数。公式(1)是对历史数据的平滑（smooth，filter），公式(2)是对未来数据的预测（predict）。求解 $\hat f = g\left( {{s_1},{s_2},...,{s_t}} \right)%$ 的过程是数据模型的更新。正如绪论中讨论的那样，数据模型可以形式化为，

$\hat f = g\left( {s,s',s'',...} \right)%$

其中s为观测直接得到的数据， $s'%$ 为观测数据的一阶微分或者偏微分， $s''%$ 为二阶微分或者偏微分，省略的部分为更高阶的微分或者偏微分。假如模型的复杂度函数h和模型涉及的数据的阶数相关，阶数越小复杂度越小，阶数越高复杂度越高。估计需要满足的第三个公式是

$\min h\left( {\hat f} \right)%$ (3)

模型的复杂度控制就是正则化。

B将D指导E跟踪F轨迹这一问题抽象为这样一个模型，其中涉及观察变量和状态变量。观察变量 $Z_{t}=\left [x\; y \right ]_{t}^{T}%$ 是终端C得到的GPS数据，状态变量 $X_{t}=\left [x\; y\; \dot{x}\; \dot{y} \right ]_{t}^{T}%$ 是程序D用于估计F位置的。Ft-1是状态转移矩阵，描述F的运动模型。Gt-1是控制矩阵， $U_{t-1}=\left [\ddot{x}\; \ddot{y} \right ]^{T}%$ 是外控制变量。Ht是观测矩阵，描述观测和状态之间的关系。Wt-1和Vt是高斯白噪声，covariance分别是Q和R，假设其不随状态变量变化。

${X_t} = {F_{t - 1}}{X_{t - 1}} + {G_{t - 1}}{U_{t - 1}} + {W_{t - 1}}%$ (4)

${Z_t} = {H_t}{X_t} + {V_t}%$ (5)

Ft和Ht的如何确定的呢？我们首先插入一段广告。对于一个具有n阶导数的函数f，其在x处的泰勒展开为

$f\left( x \right) = f\left( {{x_0}} \right) + \left( {x - {x_0}} \right)f'\left( {{x_0}} \right) + \frac{1}{2}{\left( {x - {x_0}} \right)^2}f''\left( {{x_0}} \right) + O\left( {{{\left( {x - {x_0}} \right)}^3}} \right)%$ (6)

忽略2阶以上的项，取x=t，x0=t-1，则上式可以写成

$f\left( t \right) \buildrel\textstyle.\over= f\left( {t - 1} \right) + f'\left( {t - 1} \right) + \frac{1}{2}f''\left( {t - 1} \right)%$ (7)

对上式分别求1阶和2阶导数有

$f'\left( t \right) \buildrel\textstyle.\over= f'\left( {t - 1} \right) + f''\left( {t - 1} \right)%$ (8)

$f''\left( t \right) \buildrel\textstyle.\over= f''\left( {t - 1} \right)%$ (9)

用矩阵的形式重写公式（7）如下

$\begin{pmatrix}f\\ f'\\ f''\\\end{pmatrix}_{t} \doteq \begin{pmatrix} 1& 1 & \frac{1}{2}\\ 0& 1 &1 \\ 0& 0 &1 \end{pmatrix}\begin{pmatrix}f\\ f'\\ f''\\\end{pmatrix}_{t-1}%$ (10)

对于离散模型的，微分用差分近似表示，式（8）（9）改写为

$f'\left( t \right) = f\left( t \right) - f\left( {t - 1} \right)%$ (11)

$f''\left( t \right) = f'\left( t \right) - f'\left( {t - 1} \right)\\ \;\;= \left( {f\left( t \right) - f\left( {t - 1} \right)} \right) - \left( {f\left( {t - 1} \right) - f\left( {t - 2} \right)} \right)\\ \;\;= f\left( t \right) - 2f\left( {t - 1} \right) + f\left( {t - 2} \right)%$ (12)

$\begin{pmatrix}f\\ f'\\ f''\\\end{pmatrix}_{t} \doteq \begin{pmatrix} 1& 0 &0 \\ 1& -1 &0 \\ 1& -2 &1 \end{pmatrix}\begin{pmatrix}f_{t}\\ f_{t-1}\\ f_{t-2}\\\end{pmatrix}_{t-1}%$ (13)

式（13）给出了要估计的函数、导数与观测数据之间的关系。

广告时间结束，言归正传Ft描述的是状态之间的关系。该关系受到运动学的基本关系式的约束。牛顿运动学定律可以使用式（10）表示。

如果，我们想实时更新状态变量的值，式（13）告诉我们，观测数据是如何影响状态变量的。如果不想实时更新，就可以仅用式（10）。

根据式（10），Ft和Ht的具体形式为

$F_{t}=\begin{pmatrix}1 &0 &1 &0 \\ 0 &1 &0 &1 \\ 0 &0 &1 &0 \\ 0 &0 &0 &1 \end{pmatrix}%$ (14)

$G_{t-1}=\begin{pmatrix}\frac{1}{2} &0 \\ 0 &\frac{1}{2} \\ 1 &0 \\ 0 &1 \end{pmatrix}%$ (15)

$H_{t}=\begin{pmatrix}1 &0 &0 &0 \\ 0 &1 &0 &0 \end{pmatrix}%$ (16)

在t时刻，根据式（4）预测F的当前位置 ${\hat X_t}%$ ，根据式（5）得到终端C得到的GPS数据的预测值 ${\hat Z_t}%$ 。使用t-1时刻的最优状态估计 $X\left( {\left. {t - 1} \right|t - 1} \right)%$ ，代入式（4）得

$\hat X\left( {\left. t \right|t - 1} \right) = {F_{t - 1}}\hat X\left( {\left. {t - 1} \right|t - 1} \right) + {G_{t - 1}}{U_{t - 1}}%$ （17）

$X\left( {\left. {t} \right|t - 1} \right)%$ 的covariance更新如下，其中covariance用P表示

$P\left( {\left. t \right|t - 1} \right) = {F_{t - 1}}P\left( {\left. {t - 1} \right|t - 1} \right){F_{t - 1}}^T + {Q_{t - 1}}%$ （18）

$P\left( {\left. {t} \right|t - 1} \right)%$ 是 $\hat X\left( {\left. {t} \right|t - 1} \right)%$ 的covariance， $P\left( {\left. {t - 1} \right|t - 1} \right)%$ 是 $\hat X\left( {\left. {t - 1} \right|t - 1} \right)%$ 的covariance。式（17）、（18）完成了预测，如何结合新的观测求解最优估计呢，继续往后看。t时刻的观测变量的预测

$\hat Z\left( {\left. t \right|t - 1} \right) = {H_t}\hat X\left( {\left. t \right|t - 1} \right)%$ （19）

观测变量的covariance

${S_t} = {H_t}P\left( {\left. t \right|t - 1} \right){H_t}^T + {R_t}%$ （20）

Kalman增益

${K_t} = P\left( {\left. t \right|t - 1} \right)H_t^TS_t^{ - 1}%$ （21）

t时刻观测变量的真实值与预测值之间的残差

${v_t} = {Z_t} - \hat Z\left( {\left. t \right|t - 1} \right)%$ （22）

t时刻观测变量的最优估计

$\hat X\left( {\left. t \right|t} \right) = \hat X\left( {\left. t \right|t - 1} \right) + {K_t} \cdot {v_t}%$ （23）

其covariance的最优估计是

$P\left( {\left. t \right|t} \right) = \left( {I - {K_t}{H_t}} \right)P\left( {\left. t \right|t - 1} \right)%$ （24）

公式（17）-（24）可以使用图3解释（图3使用观测变量来表示，而没有具体描述状态变量的预测和寻优的过程）。公式（17）-（19）对应的是图3b，根据图3a的t-1时刻的最优估计预测t时刻的观测变量。公式（20）-（24）对应的是图3d，根据图3c的新观测变量计算最优状态变量和观测变量。

$P\left( {\left. t \right|t - 1} \right) = E\left( {\left( {{X_t} - \hat X\left( {\left. t \right|t - 1} \right)} \right){{\left( {{X_t} - \hat X\left( {\left. t \right|t - 1} \right)} \right)}^T}} \right)%$ （25）

$P\left( {\left. t \right|t} \right) = E\left( {\left( {{X_t} - \hat X\left( {\left. t \right|t} \right)} \right){{\left( {{X_t} - \hat X\left( {\left. t \right|t} \right)} \right)}^T}} \right)%$ （26）

首先根据状态转移模型计算状态值的预测，求得观测变量的预测值。然后获得新的观测变量。再结合观测变量和观测变量的预测值，求出状态和观测变量的最优估计值。下面给出的是t时刻最优估计的模型，依然是高斯的。

$E\left( {{X_t}} \right) = \hat X\left( {\left. t \right|t} \right)%$

$E\left( {\left( {{X_t} - \hat X\left( {\left. t \right|t} \right)} \right){{\left( {{X_t} - \hat X\left( {\left. t \right|t} \right)} \right)}^T}} \right) = P\left( {\left. t \right|t} \right)%$

$p\left( {\left. {{X_t}} \right|{Z_t}} \right) \sim N\left( {\hat X\left( {\left. t \right|t} \right),P\left( {\left. t \right|t} \right)} \right)%$

a b c d

图3 F位置估计

Kalman Filter的matlab代码

[plain] view plaincopy
% the data to estimate  
lens=100;  
a=2;  
b=50;  
x=1:lens;  
y=a*x+b*randn(1,lens);  
D=[x;y];  
   
% the number of the stateparamters  
StateParamNum=4;  
% the number of thecontrol parameters  
ContrParamNum=2;  
% the number of theobservation parameters  
ObsevParamNum=2;  
% the motion transitionmatrix  
F=[1 0 1 0;0 1 0 1;0 0 1 0;0 00 1];  
% the control matrix  
G=[0.5 0;0 0.5;1 0;0 1];  
% the observation matrix  
H=[1 0 0 0;0 1 0 0];  
% the state vector  
X=zeros(StateParamNum,1);  
X=[D(1,1);D(2,1);0.001;0.001];  
% the control vector  
U=0*randn(2,1);  
% the observation vector  
Z=zeros(ObsevParamNum,1);  
% the covariance of thestate  
P=eye(StateParamNum,StateParamNum);  
P(1,1)=10;  
P(2,2)=10;  
P(3,3)=10;  
P(4,4)=10;  
% the covariance of thestate noise  
q=eye(StateParamNum,StateParamNum);  
q(1,1)=0.1;  
q(2,2)=0.1;  
q(3,3)=0.01;  
q(4,4)=0.01;  
% the covariance of theobserve noise  
r=eye(ObsevParamNum,ObsevParamNum);  
r(1,1)=10;  
r(2,2)=10;  
% the optimal estimationof the the state  
Xf=zeros(StateParamNum,lens);  
% the optimal estimationof the the observation  
Zf=zeros(ObsevParamNum,lens);  
V=zeros(ObsevParamNum,lens);  
Pf=zeros(StateParamNum,lens);  
   
for i=1:lens  
    % theestimation of the state in time t  
    Xest=F*X+G*U;  
    % thecovariance of the estimated state  
    Pest=F*P*F'+q;  
    % theestimation of the observation in time t  
    Zest=H*Xest;  
    % thecovariance of the estimated observation  
    Sest=H*Pest*H'+r;  
    % theKalman Gain  
    K=Pest*H'*inv(Sest);  
    % thedifference between estimation and observation  
    v=D(:,i)-Zest;  
    % theoptimal estimation of the state in time t  
    X=Xest+K*v;  
    % thecovariance of the optimal state  
   P=(eye(StateParamNum,StateParamNum)-K*H)*Pest;  
    % theoptimal estimation of the observation in time t  
    Z=H*X;  
     
    Xf(:,i)=X;  
    Zf(:,i)=Z;  
    V(:,i)=v;  
    Pf(:,i)=diag(P);  
end  
   
figure(1)  
hold on  
colormax=lens+1;  
c=(1:lens-1)/colormax;  
c1=repmat(c,3,1);  
c2=[ones(1,lens-1);repmat(c,2,1)];  
   
for i=1:lens-1  
    plot(D(1,i:i+1),D(2,i:i+1),...  
        'LineWidth',3,'Color',c1(:,i)');  
    plot(Zf(1,i:i+1),Zf(2,i:i+1),'-+',...  
        'LineWidth',1,'Color',[1 0 0]);  
end  
hold off  

a b

图4 a、原始数据和滤波数据；b、拟合误差

上面这段代码中去掉了控制部分，调整q和r可以改善滤波的效果。

在上面的推导过程中，对Kalman滤波有以下几点认识

1、模型是线性的，体现在公式（4）中；

2、模型是高斯的，体现在公式（4）、（5）中；

3、式（4）中的F和G矩阵没有更新；

4、状态变量是根据设计者的知识给出的。

综合以上几点，Kalman滤波是一个预设的跟踪器，物体的运动模型，运动之间的关系都是给定的。

我们可以默认这些预设都是正确的，直接来用。但是，人作为第一发现者，是如何从数据中抽象出这些状态，如何从状态到状态的转移求得运动模型，这些都没有解决。如果，这个问题没有解决，我们将无法进入下一个螺旋上升的阶段。说到这里，我不得不怀疑这样一点——我们获得的所谓的运动模型，所谓的状态变量的内容，是不是某些非人类的文明灌输给我们的。如果是这样的，那么我们不可能发现螺旋上升的途径；如果不是这样，这些是我们自己发现的，那么我们就有办法重新发现“发现知识”的过程，指导我们进入下一个螺旋。

1 0