最新消息:雨落星辰是一个专注网站SEO优化、网站SEO诊断、搜索引擎研究、网络营销推广、网站策划运营及站长类的自媒体原创博客

在R中模拟复合泊松过程

SEO心得admin144浏览0评论
本文介绍了在R中模拟复合泊松过程的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧! 问题描述 我正在尝试模拟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$

重复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)))
发布评论

评论列表(0)

  1. 暂无评论