卡尔曼滤波器小结

      断断续续，大概弄了将近一个月的时间，今天总算是把卡尔曼滤波器用于流量矩阵估计的问题给解决了。实验结果很理想，和作者的实验结果差不多。

      回顾整个过程：

      1，觉得卡尔曼有现成的预测和更新公式，只要确定参数C、Q、R就可以使用，很简单，便用最小二乘从直接测量的一天的数据算出C，再根据样本方差确定Q，由于Y是直接通过恒定的Y计算得来，所以假设R=0。但是，通过这些参数，使用卡尔曼进行估计时，发现估计的值是发散的，相对误差几乎成指数增长，令人费解。

      2，反省，觉得还是由于对卡尔曼理解不够，所以导致不知道问题出在哪。于是开始着手对卡尔曼滤波的两个方程的理解，包括公式的推导。主要是其中的增益矩阵K的推导，soul的那篇论文里并没有给出推导公式，我又看了其参考文献的关于卡尔曼的介绍，也没有给出推导，后来反而是在网上瞎找时，偶然在维基百科上找到了关于卡尔曼滤波的完整介绍，包括各个公式所依据的理论和推导过程，知道原来卡尔曼的理论基础原来是马尔可夫过程和贝叶斯统计。通过这上面的介绍也消除了我以前在看卡尔曼的一个疑问，就是关于更新公式中的估计方差的公式不一致的问题，有的论文（流量估计方面的论文中关于卡尔曼滤波）提出Pe = （I-KA）Pp（I - KA）'   （1）而有的论文（专门单独介绍卡尔曼滤波的论文）为Pe = （I-KA）Pp   （2）。其实两种写法都对，当K取最优增益时（使估计误差最小），(1)可以通过一系列推导简写为（2） ，当然在实际应用时，K都是取的最优增益，所以便取（2）。

       3，通过2对卡尔曼滤波的理解，我觉得可能是由于参数的确定出了问题，因为soul提出用最大似然估计的EM算法确定C和Q，而我用的是最小二乘，而我们知道卡尔曼滤波在使用时是假设X（t）=CX（t-1）+W（t），即X为马尔可夫链式的高斯分布，所以明显地用最大似然估计得到的参数更适用于卡尔曼估计。今天的实验结论也证明了这一点。但是问题又来了，soul本人的论文中关于怎么样求参数这一问题却写的相当简略而且有点自相矛盾，而实际上，卡尔曼滤波的参数确定需要根据实际的问题来决定，这也是卡尔曼滤波应用的难点和关键点。soul的混乱与矛盾在于，论文的第一段说通过24小时直接测量到的X 算出C，第二段又笔锋一转：只测量X0，根据观测数据Y，通过EM算法确定参数C、Q的最大似然估计。这二段文字我读了很多遍，联系上下文揣摩了很多遍，觉得还是矛盾，后来我又在其参考文献里专门琢磨其EM算法，由于引入了SWITCH，导致公式符号复杂化，看起来很吃力，后来终于搞定。为了弄清楚EM算法，我又专门看了些专门介绍EM算法的论文，茹正亮的《EM算法在不完全参数估计中的应用》写的比较好，至少容易看懂，至此我明白了EM算法的原理，知道了soul论文中的矛盾不矛盾的唯一可能性：如果EM算法的的输入数据包括24小时直接测量到的X和Y，那么这个矛盾就不矛盾了。但我知道这是不可能的，因为EM算法在E步和M步递推时，必须要借助于隐藏数据X的概率密度，这样才能从上一下的参数确定出下一步的参数。

      4，在这一矛盾下，我选择了后者（只测量到X的一个初始态和24小时的Y，假设其余的X为未知的马尔可夫高斯分布），因为作者对于后者的介绍的多一些，于是我用代码直接编写的参数估计的EM最大似然估计算法，其具体的EM公式实在是很难看懂，运行后发现结果很不理想，因为中间根据上一步参数求出下一步参数时，需要用到卡尔曼平滑器（根据24小时的Y对X和P进行平滑处理），这样首先又需要应用卡尔曼滤波求出预测的估计的所有X，算法代码在运行时往往会发散。我昨天去图书馆借了几本关于卡尔曼的书，了解到原来是当P在小型机上计算时中间有一个矩阵减法，这样容易导致其不正定，进而导致迭代计算时会发散。由于初始参数的选取我是采用了过程1中提到的最小二乘算出的C和Q，所以便出现了1中出现的发散现象。总之EM算法对参数的初始值很敏感，这样，初值的选定成了问题。这一部分现在我还没有搞清楚。

      5，我决定走自己的路，选择矛盾的前者，即通过直接测量的24小时的X，通过最大似然估计求出C才Q（假设R=0，即Y的测量无误差），这一条路其实中间我尝试过，似然函数为矩阵C和Q的函数，由于无法用EM算法求解，所以只好直接求偏导算极值，但是遇到了函数矩阵的问题，即把C和Q这样的矩阵看作变量，如何求导？查阅矩阵分析的书，都是介绍矩阵函数的（矩阵的元素为函数），昨天在图书馆看到一本特殊矩阵的书，里面才有提到函数矩阵的概念，和一些性质（这时候我才知道这个叫函数矩阵）。至此我已经觉得彻底无路可走了，因为函数矩阵的问题我是无法搞定了。

      6，转念一想，既然写成矩阵无法求解，那我就再将矩阵的各个元素分别列出来，独立求解，昨天晚上推导了一会，发现了可行性，便编写了代码验证，由于算法比较粗糙，用MATLAB运行了1个多小时，未出结果，因为矩阵过大，MATLAB存不下。今天上午重新对算法改进，一行一行的求解C，这样终出求了出来。根据C再求出Q，将C和Q代入之前的卡尔曼，经过对一些错误的修正，终于取了满意的实验结果。

 

　　整个过程中，一方面碰到矛盾、问题时好多次都绝望了，想过放弃，改做别的方法，中间我还做了一些实验，关于相关系数平稳性的实验，发现任何周期的相关系数都极不稳定，没有任何规律可言，而OD流的周期性通过我采用的一种独特的平滑技术才表现出天的周期性。这让我觉得卡尔曼在用于TM估计时是无力的，因为２４小时估计出来的时空关系C（为常数）是粗糙的。我还据此产生的新的想法，略去卡尔曼的传输过程，直接通过Y求X，运用平滑技术的周期性，但广义逆的问题不好解决，这样做出来的东西是个什么东西，没有理论支持，就暂且搁下了。

　　今天的实验结果让我看到了卡尔曼的强大，它的优秀自适应调整能力。到目前为止，我还是很惊异于它居然有这么好的实验结果。我的关于卡尔曼的改进还有没有可提升的空间？这一问题变得更容易被否定。目前周静静提出了一种UD分解的平方根卡尔曼滤波，对我的研究冲击很大。我也是因此而昨天去的图书馆，找到了卡尔曼发散的原因呵。

　　先写到这里。路还很长，而且分岔口很多。