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

1、随机模拟报告2011211173航院陈金秀在本次所有的模拟过程中，下标为 1作为开始时刻，如在泊松过程中，是用 S(1)表示 第0次跳跃的时间，取为 0。 S(i)表示第i 1次跳跃的时间，S(3)表示N(t) 2的时刻, X(i)也是采取类似的取法。1泊松过程的模拟1)对于参数为时齐泊松过程的模拟。时齐泊松过程的模拟是基于，N(t)每次增加的时间间隔， 即每次跳跃的时间间隔是独立且服从参数为的指数分布。所以时齐泊松过程的模拟本质上是模拟每次跳跃的时间间隔。方法为生成一个在 0到1上的随机数X，然后取丫In(X)，则丫就是参数为的指数分布。对于3000次的跳跃的模拟过程，就是用此方法生成300

2、0个Y(i)，则第N次跳跃的N时刻SN Y(i)。Y(0)表示表示初始时间取为0，在编程的时候，由于不指定如何用向i 0量表示Y(0),都是采用Y(i)表示第i 1次的跳跃时间间隔，其余相关符号也采用类似方法 详情可以见代码中的注释。以S为横坐标，纵坐标表示跳跃后的值，就可以得到跳跃过程，如图所示：局部放大后如下图所示:2)对于强度函数(t)的/非时齐泊松过程t根据定理2.6.2 (2),对于非时齐泊松过程，先由(t)计算函数m(t) o (s)ds，对于1的时齐泊松过程 M (u)，令N(t) M (m(t)，则N(t)就是强度函数为的泊松过程。所以采用的方法是先模拟一个参数为1的时齐泊松过

3、程，然后对于每次跳跃时间，代入im(t)的反函数，可以求得非时齐泊松过程的各次跳跃时间，本例中取(t) t ，10次跳跃，900条轨道，对首达时间进行频数统计，得到直方图如下图所示：2.计算积分10计算积分f (x)dx Ef (X)，我们可以将此积分看成f(X)的数学期望，其中X U0,1。若 Uk,1k n为独立同分布，则可以取1 f (Ui)作为n i i的估计,由大数定律当n 时，n，这表明可以用随机模拟的方法计算积分。本例中计算标准正态分布概率密度函数在3(0,3)上的积分o f(x)dxx22 dx，对于状态空间S 0,1,2,3，初始分布(0)0.2,0.3,0.2,0.3，一步

4、概率转移矩阵其参考解为(3)(0)0.4987。首先生成A U (0,1)的随机数，则(0,3)上的随机数为3 x x13X 3A，则积分式=3f(3)d() 3f(3y)dyf (3U i)。每次蒙特卡洛模拟使33n i 1用10000个随机数，并进行100次循环取每次结果的平均值作为最终计算结果。下图为某次计算结果：Valu田田田田田田1.6750e+0331DOOO1000987其中y即为最终计算结果， 与参考值完全符合。为验证算法的稳定性，运算十次发现与参考值的误差均在 0.002之内。3离散马氏链的模拟对于离散马氏链的模拟，其关 键是已知第i次取值X，判断下次跳跃后所在位置，可以生成

5、一个0到1的随机数a,根据随机变量的值和一步转移矩阵第X行的值，来确定跳跃后，即第i + 1次的取值。当a大于第X行前j-1项和，小于或等于前 j项和时，第i + 1跳跃的 取值为j。初值取法类似，只是改由初始分布n (0来确定。0.70.30000.60.40P0.5000.5000.80.2生成的一条具有19000个数据点的轨道的图像如下(由于数据点较多轨道路径稠密，所以只取局部放大便于观看)：横坐标表示跳跃次数，纵坐标表示跳跃后的值4. Metropolis-Hast ings算法实现及应用实现步骤为：1) 任取马尔可夫链的初始状态为Xo x ；2) 由转移核q(x,x)产生一个尝试移动

6、 x ；3) 生成 U (0,1)随机数 u，若 u (x,x) min 1,(x)q(x,x)，则令 x1 x，否则(x)q(x,x)保持当前状态不变，即 X1 X0 x ；4) 重复上述步骤，依次生成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，所以问题的关键在于如何得到(x)的独立同分

7、布的样本点。 而Metropolis-Hastings方法可以得到收敛到(x)的马尔可夫链，例如模拟产生10000条长度为10000的马氏链，那么取每条链的最后一个值就可以视为得到了一个来自于极限分布(X)的有10000个独立样本点的样本，然后就可以进行随机模拟。本例为计算标准正态分布的高阶原点矩E(Xk)，当k为奇数时易知其为0，当k为偶数时，本例中仅以4阶原点矩为例。采用独立抽样,q(x,x) q(x)。为使独立抽样能有较好的结果，选取的q(x)应与(x)接近，本例中选用类指数分布作为转移核。F图为抽样分布与目标分布的比较，可以看到符合程度很好。700600600400003001005-4-3-2-101234具体程序参见源代码，下图为某一次的计算结果:ManneValueM4SXY hIIIJrratio2.97&了2r9767e+04*150C0zl doublelxl0(00 double0150310000LOOOOLOOOO1.0300151 doubitM4即为四阶原点矩计算值。为检验算法的准确性，查阅相关文献可知标准正态分布的n阶原点矩为E(Xn)0(n 1) (n 3)|ll 3 1 (n 1)!n为奇数n为偶数故四阶原点矩的精确值为3，计算值与其十分