【问题标题】:Simulate Compound poisson process in r在r中模拟复合泊松过程
【发布时间】:2017-06-27 18:36:27
【问题描述】:

我正在尝试在 r 中模拟复合泊松过程。该过程由 $ \sum_{j=1}^{N_t} Y_j $ 定义,其中 $Y_n$ 是独立于 i.i.d 序列的 $N(0,1)$ 值,$N_t$ 是带有参数 $1$ 的泊松过程。我试图在没有运气的情况下在 r 中模拟这个。我有一个算法来计算它,如下所示: 模拟cPp从0到T:

启动:$ k = 0 $

在 $\sum_{i=1}^k T_i 时重复

设置$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 中模拟这个,以便我可以绘制过程吗?我已经尝试过,但无法完成。

【问题讨论】:

  • 你能使用rnormrexp,while吗?它可能很慢,但与任何其他编程语言都没有区别。你试过什么?
  • 您好,我也有同样的问题,您如何从任何分布中模拟 Y?

标签: r simulation


【解决方案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")

这个算法,因为它依赖于低级优化函数,所以速度很快:我测试它的六年历史的系统每秒会生成超过 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)))

【讨论】:

  • 非常感谢,这对我帮助很大。有没有可能得到一个更好的情节?即不是在每个点周围都显示一个圆圈?
  • 是的:请参阅位于?plot.stepfun 的帮助页面。它告诉您在最后一行调用plot 时包含参数do.points=FALSE,并描述了用于控制绘图外观的其他选项。
猜你喜欢
  • 1970-01-01
  • 2019-09-15
  • 1970-01-01
  • 2020-08-15
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2023-03-23
  • 2011-11-01
相关资源
最近更新 更多