我正在尝试在 r 中模拟复合泊松过程。该过程由 $ \sum_{j=1}^{N_t} Y_j $ 定义,其中 $Y_n$ 是独立同分布序列独立的 $N(0,1)$ 值,$N_t$ 是参数为 $1$ 的泊松过程。我试图在 r 中模拟这个,但没有运气。我有一个算法来计算它,如下所示:模拟从 0 到 T 的 cPp:
启动:$ k = 0 $
当 $\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) = su(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")
Run Code Online (Sandbox Code Playgroud)
该算法由于依赖于低级优化函数,因此速度很快:我测试它的六年前的系统每秒将生成超过三百万个(时间,值)对。
这对于模拟来说通常已经足够好了,但它并不能完全满足问题,该问题要求在时间 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)))
Run Code Online (Sandbox Code Playgroud)