泊松、非时齐泊松、离散马氏链仿真.docx

随机模拟报告 2011211173 航院 陈金秀在本次所有的模拟过程中，下标为1作为开始时刻，如在泊松过程中，是用表示第0次跳跃的时间，取为0。表示第i－1次跳跃的时间，表示的时刻，也是采取类似的取法。1.泊松过程的模拟1）对于参数为时齐泊松过程的模拟。 时齐泊松过程的模拟是基于，每次增加的时间间隔，即每次跳跃的时间间隔是独立且服从参数为的指数分布。所以时齐泊松过程的模拟本质上是模拟每次跳跃的时间间隔。方法为生成一个在0到1上的随机数X，然后取，则Y就是参数为的指数分布。对于3000次的跳跃的模拟过程，就是用此方法生成3000个，则第N次跳跃的时刻。表示表示初始时间取为0，在编程的时候，由于不指定如何用向量表示，都是采用表示第i－1次的跳跃时间间隔，其余相关符号也采用类似方法,详情可以见代码中的注释。 以S为横坐标，纵坐标表示跳跃后的值，就可以得到跳跃过程，如图所示：局部放大后如下图所示：2）对于强度函数的/非时齐泊松过程根据定理2.6.2（2），对于非时齐泊松过程，先由计算函数，对于的时齐泊松过程，令，则就是强度函数为的泊松过程。所以采用的方法是先模拟一个参数为1的时齐泊松过程，然后对于每次跳跃时间，代入的反函数，可以求得非时齐泊松过程的各次跳跃时间，本例中取，10次跳跃，900条轨道，对首达时间进行频数统计，得到直方图如下图所示：2.计算积分计算积分，我们可以将此积分看成的数学期望，其中。若为独立同分布，则可以取作为的估计，由大数定律当时，，这表明可以用随机模拟的方法计算积分。本例中计算标准正态分布概率密度函数在上的积分，其参考解为。首先生成的随机数，则上的随机数为，则积分式=。每次蒙特卡洛模拟使用10000个随机数，并进行100次循环取每次结果的平均值作为最终计算结果。下图为某次计算结果：其中y即为最终计算结果，与参考值完全符合。为验证算法的稳定性，运算十次发现与参考值的误差均在之内。3.离散马氏链的模拟对于离散马氏链的模拟，其关键是已知第i次取值X，判断下次跳跃后所在位置，可以生成一个0到1的随机数a，根据随机变量的值和一步转移矩阵第X行的值，来确定跳跃后，即第i＋1次的取值。当a大于第X行前j-1项和，小于或等于前j项和时，第i＋1跳跃的取值为j。初值取法类似，只是改由初始分布π(0)来确定。对于状态空间，初始分布，一步概率转移矩阵生成的一条具有19000个数据点的轨道的图像如下（由于数据点较多轨道路径稠密，所以只取局部放大便于观看）：横坐标表示跳跃次数，纵坐标表示跳跃后的值4. Metropolis-Hastings算法实现及应用实现步骤为：任取马尔可夫链的初始状态为；由转移核产生一个尝试移动；生成随机数，若，则令，否则保持当前状态不变，即；重复上述步骤，依次生成。应用举例：利用Metropolis-Hastings算法对复杂积分进行模拟计算。其基本原理是对于，若为某一概率密度函数，则，即可以通过从中产生独立同分布的随机数，根据大数定律依概率收敛到，所以问题的关键在于如何得到的独立同分布的样本点。而Metropolis-Hastings方法可以得到收敛到的马尔可夫链，例如模拟产生10000条长度为10000的马氏链，那么取每条链的最后一个值就可以视为得到了一个来自于极限分布的有10000个独立样本点的样本，然后就可以进行随机模拟。本例为计算标准正态分布的高阶原点矩，当k为奇数时易知其为0，当k为偶数时，本例中仅以4阶原点矩为例。采用独立抽样， 。为使独立抽样能有较好的结果，选取的应与接近，本例中选用类指数分布作为转移核。下图为抽样分布与目标分布的比较，可以看到符合程度很好。具体程序参见源代码，下图为某一次的计算结果：其中，M4即为四阶原点矩计算值。为检验算法的准确性，查阅相关文献可知标准正态分布的n阶原点矩为 故四阶原点矩的精确值为3，计算值与其十分接近，另外为检验算法的稳定性，多次重复计算得到10次计算结果如下：2.9767 3.0709 3.1102 3.0173 2.9221 3.0458 3.2083 2.8625 2.9100 2.8717 可见算法的稳定性也较好，单次运算值即呈现出良好的模拟性能。