启动:$k=0$
重复While$sum_{i=1}^k T_i<;T$
设置$k=k+1$
模拟$T_k sim exp(Lambda)$(在我的例子中,$lambda=1$)
模拟$Y_k sim N(0,1)$(这只是一个特例,我希望能够将其更改为任何分布)
轨迹由$X_t=sum_{j=1}^{N_t}Y_j$给出,其中$N(T)=sup(k:sum_{i=1}^k T_i leq t)$
有人可以帮助我在r中模拟这一过程,以便我可以绘制过程图吗?我试过了,但做不到。
推荐答案将cumsum用于确定时间N_t和X_t的累计和。此说明性代码指定要模拟的次数,n,模拟n.t中的时间和x中的值,并(显示其所做的)绘制轨迹。
n <- 1e2 n.t <- cumsum(rexp(n)) x <- c(0,cumsum(rnorm(n))) plot(stepfun(n.t, x), xlab="t", ylab="X")由于该算法依赖于低级优化函数,因此速度很快:我测试了6年的系统每秒将生成300多万对(时间、值)。
这对于模拟来说通常是足够好的,但它不能完全满足问题,该问题要求在时间T之前生成模拟。我们可以利用前面的代码,但是解决方案有点棘手。它计算时间T之前在泊松过程中发生的次数的合理上限。它生成到达间隔时间。它被包装在一个循环中,该循环将在实际未达到时间T时重复(罕见)事件中的过程。
增加的复杂性不会改变渐近计算时间。
T <- 1e2 # Specify the end time T.max <- 0 # Last time encountered n.t <- numeric(0) # Inter-arrival times while (T.max < T) { # # Estimate how many random values to generate before exceeding T. # T.remaining <- T - T.max n <- ceiling(T.remaining + 3*sqrt(T.remaining)) # # Continue the Poisson process. # n.new <- rexp(n) n.t <- c(n.t, n.new) T.max <- T.max + sum(n.new) } # # Sum the inter-arrival times and cut them off after time T. # n.t <- cumsum(n.t) n.t <- n.t[n.t <= T] # # Generate the iid random values and accumulate their sums. # x <- c(0,cumsum(rnorm(length(n.t)))) # # Display the result. # plot(stepfun(n.t, x), xlab="t", ylab="X", sub=paste("n =", length(n.t)))