7.3.2 非齐次泊松过程的生成
本小节内容参考文献[15]第69页.
假设我们希望模拟强度函数为λ(t)的非齐次泊松过程在时刻T前的状态.第一种方法就是文献[15]给出的稀松或随机抽样法.在选取满足
的λ之后,如在文献[15]第2章证明的那样,我们可以通过随机选取参数为λ的泊松过程的事件发生时间来生成非齐次泊松过程.也就是说,如果一个参数为λ的泊松过程在时刻t发生的事件以概率λ(t)/λ被记录,则被记录的事件列就构成一个强度为λ(t),0≤t≤T的非齐次泊松过程.于是,我们可以通过对泊松过程的事件进行随机计数来产生所需的非齐次泊松过程,其算法如下：
步骤1：t=0,I=0;
步骤2：生成一个随机数U;
步骤4：生成一个随机数U;
步骤5：如果U≤λ(t)λ,则令I=I+t,S(I)=t;
步骤6：转至步骤2.
在上述算法中,λ(t)是强度函数,λ满足λ(t)≤λ;I的最终值表示时刻T前发生的事件数;S(1),···,S(I)是事件发生的时刻.
上述算法之所以被称为稀松法,是由于它仅取齐次泊松过程的部分事件.如果以被记录的事件数为有效与否的标准,则当λ(t)与λ在整个区间上相差不多时,上述算法是非常有效的.于是,若把整个区间划分为若干子区间,然后在子区间上应用上述算法,则可以提高上述算法的效率.现选取适当的k,ti和λi,使得0=t0＜t1＜t2＜···＜tk＜tk+1=T,且
如果ti-1≤s＜ti,i=1,···,k+1为生成区间(ti-1,ti)上的非齐次泊松过程,我们先生成参数为λi的指数随机变量,然后以概率λ(s)/λi接受在s∈(ti-1,ti)时刻发生的事件.由指数分布的无记忆性和其参数的数乘性可知,从一个子区间到另一个子区间,其效率并没有损失.也就是说,如果现在我们考虑的区间为(ti-1,ti),且X表示参数为λi的指数分布随机数,即t+X＞ti,则λi[X-(ti-t)]/λi+1可作为下一区间的参数为λi+1的指数分布随机数.
对于强度为λ(s)且满足关系式(7.3.8)的非齐次泊松过程,我们可用下述算法来生成它的随机数.在下面算法中,t表示当前时刻;J表示区间(若tj-1≤t≤tj,则J=j);I表示前面发生的事件数;S(1),···,S(I)表示事件发生的时间.
步骤1：t=0,J=1,I=0;
步骤3：如果t+X＞tJ,则转至步骤8;
步骤4：t=t+X;
步骤5：生成一个随机数U;
步骤6：如果U≤λ(t)/λJ,则I=I+1,S(I)=t;
步骤7：转至步骤2;
步骤8：如果J=k+1,则算法停止;
步骤9：X=(X-tJ+t)λJ/λJ+1,t=tJ,J=J+1;(https://www.daowen.com)
步骤10：转至步骤3.
对于区间(ti-1,ti),记
且假设λi＞0.此时,我们不直接应用稀松算法,而是先在第i个区间产生参数为λi的泊松过程,之后再模拟强度为λ(s)=λ(s)-λi,s∈(ti-1,ti)的非齐次泊松过程.注意,虽然为模拟此泊松过程而生成的最后一个指数随机变量超出了要求区间,但也不必浪费掉,因为它经适当变换后还可被重新使用.把上面产生的两个过程叠加起来,则可得到此区间所需的过程.我们这样做的原因是它将节省产生均值为λi(ti-ti-1)的泊松分布随机数所需的均匀分布随机变量.例如,我们考虑
利用λ=11的稀松算法,可能生成其期望数为11的事件列,其中每一个事件都需要用一个随机数来确定此事件是否被接受.换言之,先产生参数为10的泊松过程与参数为λ(s)=s,0＜s＜1的非齐次泊松过程(利用参数为λ=1的稀松算法),再把二者合并后得到的事件数与所要求的同分布,但需要验证被接受事件数的期望等于11.
产生强度为λ(t),t＞0的非齐次泊松过程的另一个方法就是直接产生相继事件的发生时间.记S1,S2,···为此过程相继事件的发生时间.由于这些随机变量显然是相关的,故我们在生成它们时是先生成S1,然后再利用S1的生成值去模拟S2,以此类推.
注意,如果一个事件在时刻s发生,那么它独立于时刻s前出现的事件,且到下一个事件发生的时间间隔的分布Fs为
于是,我们可以由分布F0产生S1,如果记其值为s1,则S2的模拟值为s1与分布Fs1产生的随机数之和,如果记此和为s2,则S3的模拟值为s2与分布Fs2产生的随机数之和,以此类推,就得到了事件发生时间S1,S2,···的模拟值.上述由分布进行随机数模拟的方法依赖其分布形式.在下面的例子中,因为分布Fs的逆函数容易求得,故我们可以利用逆变换法.
例7.3.2设λ(t)=1/(t+a),t≥0,其中a为正常数,则
于是,由式(7.3.9)有
或等价地有
即
因此,事件相继的发生时间S1,S2,···可用如下方式产生：生成随机数U1,U2,···后,依次取
一般地,
