张悦++王新梅1++伍恒
摘 要：马尔可夫过程是一个具备了马尔可夫性质的随机过程。具备离散状态的马尔可夫过程，称为马尔可夫链。马尔科夫链通常被用来建模排队理论与统计学中的建模，同时还能够运用到信号模型的熵编码技术中。在文中则主要是将其运用到院子轨迹的预测算法之中。本文描述的是基于马尔可夫链围绕预测化学原子运动轨迹设计的算法。这种预测算法旨在节省分子动力学中高昂的运算成本。它通过学习回溯并纠正的机制对状态转移概率矩阵进行更新，并进行原子运动轨迹的预测。最后用实验验证了该算法的有效性。
关键词：马尔可夫链 马尔可夫模型 原子运动轨迹 分子动力学
中图分类号：O211.62 文献标识码：A 文章编号：1674-098X（2016）7（b）-0000-00
Abstract：The Markov process meant a random process with Markov properties. A Markov process with discrete state was called Markov chain. This paper was to forecast the algorithm of chemical atomic trajectory design on the basis of Markov chain. This kind of forecast method aimed to save the expensive operation cost in molecular dynamics. A state transition probability matrix， which tends to be stable， was obtained by learning backtracking and correction mechanism to forecast the atomic trajectory. Finally， the experiment proved the effectiveness of the algorithm.
Key words： markov chain； markov model； atomic trajectory； molecular dynamics
轨迹分析作为对物体的行为进行学习或描述的基本问题现今已经成为了一个十分热门的话题。本文设计的算法用于预测原子的运动轨迹，在这一领域，还未有使用马尔可夫链进行预测的研究成果。但在普遍的运动轨迹预测领域，对马尔可夫链的使用是十分广泛的。
分子动力学依靠计算机来模拟分子、原子的运动。它是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中的能量、速度、加速度等。但是其高昂的计算成本和不断增长的用户需求之间的矛盾逾演逾烈。因此分子动力学仿真的优化技术研究一直是数学界、计算化学界以及计算机界的一个研究热点[1]。在分子动力学中，通常，分子、原子的轨迹是通过数值求解牛顿运动方程得到的。正因为在分子动力学中假定牛顿运动方程决定原子的运动，这也意味着原子的运动可能会产生特定的轨道，因此我们可以考虑能否利用历史数据对原子运动轨迹进行预测，来节省大量运算产生的开支。
本文提出了一个基于齐次马尔可夫链C-K方程，通过回溯纠正的机制来更新马尔可夫链状态转移概率矩阵的算法。通过纠正更新来使状态转移概率矩阵适应数据源，并趋于稳定。以稳定的状态转移概率矩阵来预测数据，尽可能减少分子动力学仿真中的运算量。
本文第二章将对马尔可夫链的相关概念和公式进行介绍；第三章将对整个算法进行详细描述和分析；第四章将对算法进行数据的测试并且分析；第五章对整个研究进行总结。
1 马尔可夫链
1.1 马尔可夫链与马尔可夫过程
在一个事件的发展过程中，若事件的每次状态的转移都只和前一时刻的状态有关，和过去的状态无关，或者说事件的状态转移过程是无后效性的，则这样的状态转移过程就称为马尔可夫过程。举例来说，设某个随机过程的状态X可取到一个离散集合中的值，这个值随时间t变化，我们将该值表示为X（t）。此时，时间变量离散或是连续并不影响结果。假设我们将“过去的时间”集合表示为（…，p2，p1），任何“当前时间”n，以及任何“未来时间”t，这些所有时间都在X的取值范围内，有公式（2.1）所示如下。
… < p2 < p1 < n < t（2.1）
则该过程为马尔可夫过程。如果对于所有的取值（…，x（p2），x（p1），x（n），x（t）），满足公式（2.2）所示如下。
P{X（t）=x（t）|X（n）=x（n），X（p1）=x（p1），X（p2）=x（p2），...}
= P{X（t）=x（t）|X（n）=x（n）}（2.2）
则称{X（n），n=0，1，2，…}为马尔可夫链。在本文中我们讨论的是连续时间马尔可夫链，其状态发生变化的时刻是任意时刻，是连续值。
1.2 齐次马尔可夫链
马尔可夫链在n时刻的k步转移概率记为，Pij（n，n+k）[2]。若转移概率Pij（n，n+k）只与出发、到达状态i，j以及转移步数k有关，不依赖于时刻n，则将此马尔可夫链称为齐次马尔可夫链。
1.3 k步状态转移矩阵
当k为1时，Pij（1）称为一步转移概率，记为Pij。所有的k步转移概率为Pij（k），其组成的矩阵P（n）就称为马尔可夫链的k步状态转移概率矩阵。根据C-K方程[3]可得到公式（2.3）如下。
P（k）=PP（k-1）=Pk（2.3）
2 算法详解及思路
2.1 数据预处理与状态划分
经过计算得出的数据的精确度是十分高的，但在实际研究过程中我们得到精确度在小数点后五位的数据便可。我们的算法并不是将一个位置作为一个状态，而是将两个位置之间的相对位移作为一个状态。因此在我们实验的过程中，状态是有正负之分的。位移区间的基数为0.000005，因此状态区间以0.000005分段，四舍五入。如[0.000005，0.000015]之间的相对位移省略为0.00001。此时的状态划分为1状态。在数据处理并进行状态划分之后，为了避免给无用的状态分配空间，状态区间需要分为异边、同边的情况。在这里我们使用一个灵活分配空间的二维数组bound来储存状态，bound的一维固定长度为2。并且储存状态的0和1数组其首位值都用于储存状态的个数。如果状态区间异边，则0代表的一维数组储存正状态，1代表的一维数组储存负状态；如果状态区间同边且为正，状态（包括零状态）都储存在0代表的一维数组里，1代表的一维数组的首位值为0（无负状态）；如果状态区间同边且为负，状态（包括零状态）都储存在1代表的一维数组里，0代表的一维数组首位值为0（无正状态）。
2.2 原子运动轨迹模型的建立
原子运动的轨迹由位置决定，而原子运动的位置序列是按时间顺序记录下来的。基于马尔可夫链的预测主要是根据状态转移概率矩阵来计算状态之间直接转换的概率。
2.2.1 建立类马尔可夫链
根据马尔可夫链的定义，将每个原子运动的过程用一个类马尔可夫链来描述。一个原子的类马尔可夫链预测模型我们用一个三元组表示
P = （pij） = （3.1）
2.2.2 确定状态转移概率矩阵
假设我们的原子位置之间的位移既马尔可夫链中的状态点为m个。状态转移概率矩阵P为一个m×m的矩阵。pij表示从i相对位移到j相对位移的概率。我们通过统计来得到pij。假设从i相对位移到j相对位移的次数用Nij表示。如公式（3.2）所示如下。
Pij= （i≥1，j≤m）（3.2）
这是一步状态转移概率矩阵，根据公式（2.3）可得k步转移概率矩阵如公式（3.3）所示如下。
P（k）=Pk（3.3）
2.3 预测算法的思路
2.3.1 算法思路
对于原子下一步位置的预测是基于一个数据源的学习来进行的。我们的思路是假设首先取20步数据进行学习，利用这20步数据可以得到一个初始概率以及状态转移概率矩阵。当预测第21步的相对位移时，由初始概率以及状态转移概率矩阵决定2个可能的相对位移取值，第一次返回的预测值如果超出容错范围，则回溯返回第二个预测值，如果预测仍然不准确，则第21步继续学习，初始概率以及状态转移概率矩阵随之改变。如果预测值在容错范围以内则可将此预测值加入现有的数据源中，更新初始概率以及状态转移矩阵。继续预测下一步的值。
2.3.2 主要算法描述及分析
整个算法与普通的马尔可夫链预测算法比较改进的重点在于对状态转移概率矩阵的更新。因为我们的目的在于节省计算成本，那么我们可以在较少学习数据的基础上进行预测，预测3-5步，重新进行数据学习。这样既可以保证预测的准确度，亦可以保证控制程序运行的时间。
在普遍的运用齐次马尔可夫链进行预测的应用中，因为要保证状态转移矩阵必须具有一定的稳定性。因此必须以足够多的统计数据来保证预测的精确。这也就决定了马尔可夫预测模型必须具有足够统计数据的局限性。我们改进的分子动力学仿真预测算法，在基于马尔可夫链C-K方程的基础上，在数据的预测中，采用纠正并回溯的方法更新状态转移矩阵，来缩小所需的数据源，节省需要的计算开支。
3 实验分析
我们进行实验来验证改进算法的有效性。硬件环境：Intel Core i5-3210M CPU 2.49GHz，内存为12GB。软件环境：Micros-oft windows 8.1 Professional，Visual Studio 2013。
本文的实验数据是CO2分子使用Hessian-based Predictor-Correct算法计算得到，计算环境为NWChem-5.1.1[5]。
实验数据分为三份，三份分别是三次计算的结果。每一份有300组数据，每一组将有9个位置数据，分别为碳原子的三维数据和两个氧原子的三维数据。
实验中，其中20步数据作为学习数据。为了保障状态转移概率矩阵的有效性，我们每一100步重新学习一次。因此每一份300组数据都成3次进行实验。第一次学习1-20步，分别预测至50，80，100步；第二次学习101-120步，分别预测至150，180，200步；第三次学习201-220步，分别预测至250，280，300步。在实验中，预测准确的标准是：与实际数据相差小于0.00005，每一步的9维数据预测准确的个数大于等于7判定为这一步预测准确。常规的马尔可夫链预测方法没有更新这一环节，在表中显示了常规方法与改进方法之间预测效果的对比。因三份数据的情况基本类似，在此列出一份数据的测试结果。数据的预测情况如表4.1所示如下。
通过实验可以看出，预测的次数在60次540个数据时，次数准确率与步数准确率整体较高。预测的次数在80次750个数据时，次数准确率仍可达到65%以上，而步数预测准确率可达到80%以上。可见基于改进的马尔可夫链预测出的数据是可以考虑使用的。
尽管影响原子运动轨迹的因素我们无法考虑完全，并且数据的预测有一定随机性。但在超过80%的步数预测准确率的基础上，对计算的数据进行一定程度的预测来节省高昂的运算成本是可行的。
4 结语
本文提出了一种基于马尔可夫链，改进更新状态转移概率矩阵的预测算法，并详细阐述了预测的方法和步骤。这种方法为得到原子的轨迹提出了一种新的思路，利用历史原子轨迹基于改进的马尔可夫链进行预测，在学习数据较少的基础上得到一个可靠的预测结果，为节省计算的开支提供了一种可行方式。
但原子的轨迹的影响因素众多，在单一的由C-K方程计算来得到下一步的状态转移概率矩阵，并由此来计算出下一步的状态概率的这种方法，对于外界因素的考虑是不够的。对于状态的划分使用简单的数据预处理也一定程度影响了预测效果，后续的研究可以考虑引入聚类划分状态。我们仍需要进一步的研究，改进预测的方法来得到更加可靠的预测结果。
参考文献
[1] 伍恒.分子动力学仿真化研究[D].湖南：中南大学.2015：10-11.
[2] 彭曲，丁治明，郭黎敏.基于马尔可夫链的轨迹预测[J].计算机科学.2010.37（8）：189-191.
[3] 汪荣鑫.随机过程[M].西安：西安交通大学出版社，2006.
[4] 何丽.基于Markov链的用户浏览行为预测方法[J].计算机工程.2008.34（22）：32-33.
[5] Wu H， Lu S， Zhu N， et al. A high order predictor–corrector integration algorithm for first principle chemical dynamics simulations[J].Journal of Theoretical and Computational Chemistry.2016.15（01）：1650003.
[6] 夏莉，黄正洪.马尔可夫链在股票价格预测中的应用[J].商业研究.2003（10）：62-65.
[7] Gilks W R. Markov chain monte carlo[M]. New York： John Wiley & Sons，2005.
[8] Lourderaj U， Song K， Windus T L， et al. Direct dynamics simulations using Hessian-based predictor-corrector integration algorithms[J]. The Journal of chemical physics.2007.126（4）.
[9] Doob J L.Stochastic Processes[M].New York： John Wiley and Sons，1953.
