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

随机模拟报告 2011211173航院陈金秀 在本次所有的模拟过程中，下标为 1作为开始时刻，如在泊松过程中，是用 S(1)表示 第0次跳跃的时间，取为 0。 S(i)表示第i — 1次跳跃的时间，S(3)表示N(t) 2的时刻, X(i)也是采取类似的取法。 1?泊松过程的模拟 对于参数为 时齐泊松过程的模拟。 时齐泊松过程的模拟是基于， N(t)每次增加的时间间隔， 即每次跳跃的时间间隔是独立 且服从参数为 的指数分布。所以时齐泊松过程的模拟本质上是模拟每次跳跃的时间间隔。 方法为生成一个在 0到1上的随机数X，然后取丫 In(X)，则丫就是参数为 的指数 分布。对于3000次的跳跃的模拟过程，就是用此方法生成 3000个Y(i)，则第N次跳跃的 N 时刻SN Y(i)。Y(0)表示表示初始时间取为 0，在编程的时候，由于不指定如何用向 i 0 量表示Y(0),都是采用Y(i)表示第i— 1次的跳跃时间间隔，其余相关符号也采用类似方法 详情可以见代码中的注释。 以S为横坐标，纵坐标表示跳跃后的值，就可以得到跳跃过程，如图所示： 局部放大后如下图所示: 对于强度函数 (t)的/非时齐泊松过程 t 根据定理2.6.2 (2),对于非时齐泊松过程，先由 (t)计算函数m(t) o (s)ds，对于 1的时齐泊松过程 M (u)，令N(t) M (m(t))，则N(t)就是强度函数为 ⑴的泊松过 程。所以采用的方法是先模拟一个参数为 1的时齐泊松过程，然后对于每次跳跃时间， 代入 i m(t)的反函数，可以求得非时齐泊松过程的各次跳跃时间，本例中取 (t) t ，10次跳 跃，900条轨道，对首达时间进行频数统计，得到直方图如下图所示： 2.计算积分 计算积分f (x)dx Ef (X)，我们可以将此积分看成f(X)的数学期望，其中X 计算积分 f (x)dx Ef (X)，我们可以将此积分看成 f(X)的数学期望，其中 X ~U[0,1]。若 Uk,1 k n为独立同分布，则可以取 1 f (Ui)作为 n i i 的估计, 由大数定律当n 时，n ，这表明可以用随机模拟的方法计算积分。 本例中计算标准正态分布概率密度函数在3(0,3)上的积分 本例中计算标准正态分布概率密度函数在 3 (0,3)上的积分o f(x)dx x2 2 dx， 其参考解为 (3) (0) 0.4987。首先生成A~ U (0,1)的随机数，则(0,3)上的随机数为 3 x x 1 3 X 3A，则积分式=3°f(3—)d(—) 3°f(3y)dy f (3U i)。每次蒙特卡洛模拟使 3 3 n i 1 用10000个随机数，并进行100次循环取每次结果的平均值作为最终计算结果。 下图为某次 计算结果： Valu? 田田田田田田1x10000 dcubk 田田田田田田 1.6750e+03 1x100 00 doublo ^1x100 double 3 1DOOO 100 0^987 其中y即为最终计算结果， 与参考值完全符合。为验证算法的稳定性， 运算十次发现与 参考值的误差均在 0.002之内。 3?离散马氏链的模拟 对于离散马氏链的模拟，其关 键是已知第i次取值X，判断下次跳跃后所在位置， 可以 生成一个0到1的随机数a,根据随机变量的值和一步转移矩阵第 X行的值，来确定跳跃后， 即第i + 1次的取值。当a大于第X行前j-1项和，小于或等于前 j项和时，第i + 1跳跃的 取值为j。初值取法类似，只是改由初始分布 n (0来确定。 0.7 0.3 0 0 0 0.6 0.4 0 P 0.5 0 0 0.5 0 0 0.8 0.2 生成的一条具有19000个数据点的轨道的图像如下 (由于数据点较多轨道路径稠密， 所以只 取局部放大便于观看)： 横坐标表示跳跃次数，纵坐标表示跳跃后的值 4. Metropolis-Hast ings算法实现及应用 实现步骤为： 任取马尔可夫链的初始状态为 Xo x ； 由转移核q(x,x)产生一个尝试移动 x ； 生成 U (0,1)随机数 u，若 u (x,x) min 1, (x)q(x,x)，则令 x1 x，否则 (x)q(x,x) 保持当前状态不变，即 X1 X0 x ； 重复上述步骤，依次生成 X2,X3j||,XnJ( |。 应用举例：利用 Metropolis-Hastings算法对复杂积分进行模拟计算。其基本原理是对 于h(x) (x)dx，若(x)为某一概率密度函数，则 h(x) (x)dx E h(X)，即可以通 1 n 过从(x)中产生独立同分布的随机数 X11X21X3^|,Xn^|，根据大数定律一 h(Xk)依 概率收敛到 h(x) (x)dx，所以问题的关键在于如何得到