在R中模拟复合泊松过程

Question

我正在尝试在R中模拟复合泊松过程。该过程由$ \sum_{j=1}^{N_t} Y_j $定义，其中$Y_n$是独立同分布的序列，其值为$N(0,1)$，而$N_t$是参数为$1$的泊松过程。我一直在努力模拟这个过程，但没有成功。我有一个算法来计算这个过程，如下所示：

从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) = sup(k : \sum_{i=1}^k T_i \leq t )$

请问有人能帮我在R中模拟这个过程，以便我可以绘制出它的图像吗？我已经尝试过了，但无法完成。

- Slangers

---

你能否使用 rnorm，rexp 和 while？虽然可能会慢一些，但与其他编程语言没有区别。你尝试过什么？ - AdamO

嗨，我这里也有同样的问题，您如何从任何分布中模拟Y？ - G-09

Answer 1

使用 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年前的系统每秒可以生成超过三百万个（时间，值）对。

---

这通常足够进行模拟，但它并不能完全满足问题，因为问题要求在时间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)))

- whuber

---

非常感谢，这对我帮助很大。是否可能得到一个更好的图形呢？比如不显示每个点周围的圆圈？ - Slangers

是的：请参阅?plot.stepfun帮助页面。它告诉您在调用最后一行上的plot时包括参数do.points=FALSE，并描述了其他控制绘图外观的选项。 - whuber