在 R 中手动模拟泊松过程

Question

下面的问题告诉我们从ρ（到达间隔时间）和τ（到达时间）。

One of the theoretical results presented in the lectures gives the following direct method for simulating Poisson process:

• Let τ₀ = 0.
• Generate i.i.d. exponential random variables ρ1, ρ2, . . ..
• Let τ_n = ρ₁ + . . . + ρ_n for n = 1, 2, . . . .
• For each k = 0, 1, . . ., let N_t = k for τ_k ≤ t < τ_k+1.

Using this method, generate a realization of a Poisson process (N_t)_t with λ = 0.5 on the interval [0, 20].

Generate 10000 realizations of a Poisson process (N_t)_t with λ = 0.5 and use your results to estimate E(N_t) and Var(N_t). Compare the estimates with the theoretical values.

我尝试的解决方案：

首先，我使用 R 中的 rexp() 函数生成了 ρ 的值。

rhos <-function(lambda, max1)
{
    vec <- vector()

    for (i in 1:max1) 
    {
        vec[i] <- rexp(0.5)
    }

    return (vec)
}

然后，我通过 ρs 的渐进求和创建了 τs。

taos <- function(lambda, max)
{
    rho_vec <- rhos(lambda, max)
    #print(rho_vec)

    vec <- vector()
    vec[1] <- 0
    sum <- 0
    for(i in 2:max)
    {
        sum <- sum + rho_vec[i]
        vec[i] <- sum
    }

    return (vec)
}

下面的函数是求N_t=k的值，当k 给出。说，就是7，等等

Ntk <- function(lambda, max, k)
{
    tao_vec <- taos(lambda, max)
    val <- max(tao_vec[tao_vec < k])
}

y <- taos(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")

输出：

如您所见，泊松过程的图是空白的，而不是阶梯状的。

如果我将 rexp 更改为 exp，我会得到以下输出：

..这是一个阶梯函数，但所有步骤都是相等的。

为什么我的源代码没有产生预期的输出？

Answer 1

您似乎在使用 max1 来指示在 rhos 函数中对指数分布进行采样的次数。我会推荐这样的东西：

rhosGen <- function(lambda, maxTime){
  rhos <- NULL
  i <- 1
  while(sum(rhos) < maxTime){
    samp <- rexp(n = 1, rate = lambda)
    rhos[i] <- samp
    i <- i+1
  }
  return(head(rhos, -1))
}

这将继续从指数采样，直到这些保持时间的总和大于给定间隔的长度。 head 删除最后一个样本，以便我们跟踪的所有事件都肯定发生在我们感兴趣的时间间隔内。从这里你必须通过总结以前的持有时间 (rhos) 来生成 taos:

taosGen <- function(lambda, maxTime){
  rhos <- rhosGen(lambda, maxTime)
  taos <- NULL
  cumSum <- 0
  for(i in 1:length(rhos)){
    taos[i] <- sum(rhos[1:i])
  }
  return(taos)
}

现在您有了 taos，我们知道时间间隔 (0,maxTime) 中的每个事件在什么时间发生。这导致我们通过找到时间间隔内每个 t 的 Nt 值来生成相关的泊松过程：

ppGen <- function(lambda, maxTime){
  taos <- taosGen(lambda, maxTime)
  pp <- NULL
  for(i in 1:maxTime){
    pp[i] <- sum(taos <= i)
  }
  return(pp)
}

这会在间隔中的每个整数时间生成泊松过程的值。我怀疑您的部分问题是试图将 tao 值放在 y 轴上，而不是已经发生的事件数。以下代码对我有用，可以生成一个外观随机的楼梯，类似于您的示例。

y <- ppGen(0.5, 20)
x <- seq(0, 20-1, by=1)

plot(x,y, type="s")

Answer 2

这是另一种可能的实现方式。我们的想法是生成一个等待时间向量 (tau)，并根据我们正在等待的事件列表 (max1)

绘制它

poi.process <- function(lambda,n){
  
  # initialize vector of total wait time for the arrival of each event:
  s<-numeric(n+1)
  # set S_0 = 0
  s[1] <-0
  # generate vector of iid Exp random variables:
  x <-replicate(n,rexp(1,lambda))
  # assign wait time to vector s in for loop:
  for (k in 1:n){
    
    s[k+1] <-sum(x[1:k])
    
  }
  # return vector of wait time
  return(s)
  
}

使用 stepfun 绘制它会得到这样的结果：


n<-20

lambda <-3

# simulate list of wait time:
s_list <-poi.process(lambda,n)

# plot function:
plot(stepfun(0:(n-1), s_list), 
     do.points = TRUE,
     pch = 16,
     col.points = "red",
     verticals = FALSE,
     main = 'Realization of a Poisson process with lambda = 3',
     xlab = 'Time of arrival',
     ylab = 'Number of arrivals')

示例泊松过程：

在 R 中手动模拟泊松过程

Manually simulating Poisson Process in R

r

poisson

montecarlo

stochastic-process