用R和ggplot2模拟泊松过程

Question

使用速率λ= 0.7的泊松过程的模拟。展示了一个泊松过程的运行示例，其中N(t)在垂直轴上，时间t在水平轴上。模拟范围为t0:100。生成包含10条轨迹的第一张图和包含100条轨迹的第二张图。

我已经尝试了下面的代码，但是我不能同时生成两个图。

library(plyr)
library(ggplot2)

Process_poisson<- function(t, lambda){
distr_poisson<- rpois(1, t*lambda)
s_poisson<- sort(runif(distr_poisson, 0, t))
data.frame(x = c(0, 0, s_poisson),y = c(0, 0:distr_poisson))
}

N_simulations<- function(n,t,lambda){
s_poisson<- lapply (1:n, function(n) data.frame(Process_poisson(t, lambda), simulation = n))
s_poisson<- ldply (s_poisson, data.frame)
s_poisson$simulation<- factor(s_poisson$simulation)
}

t<- 0:100
lambda<- 0.7
N_simulations(10, t, lambda)
N_simulations(100, t, lambda)

par(mfrow = c(1,2))

matplot(x, y, type = "l", lty = 0:5, lwd = 1, lend = par("lend"),
     pch = NULL, col = simulation, cex = 0.5, bg = NA, main =sprintf("Nº simulations of trajectories of Poisson Process",10,lambda), xlab = "Time", ylab = "N(t)",
   xlim = c(0,100), ylim = c(-10,0))

matplot(Proceso_poisson(t, lambda), n, y, type = "l", lty = 0:5, lwd = 1, lend = par("lend"),
     pch = NULL, col = simulacion, cex = 0.5, bg = NA, main =sprintf("Nº simulations of trajectories of Poisson Process",10,lambda), xlab = "Time", ylab = "N(t)",
     xlim = c(0,100), ylim = c(-10,0))

我怎么能做到呢？

非常感谢！

Answer 1

我想你可以让这件事变得更简单。这里有一个ggplot解决方案。

首先，创建一个函数，通过使用适当的λ从指数分布中提取样本来模拟泊松过程。在本例中，我使用了一个以第一个元素为0的向量x开头的while循环。该函数通过添加随机样本来增长此向量，直到其总和达到目标持续时间tmax。这不是最有效的方法，但应该使示例更清晰。

当达到目标时，该函数返回向量的累积和，该向量表示适当λ的泊松过程的到达时间。请注意，为了简化绘图，它实际上返回了一个数据框，其中包含累计时间、累计计数和一个分组变量run，该变量允许我们在单个绘图上轻松地绘制多个运行。

make_sample_df <- function(run, tmax, lambda)
{
  x <- 0
  while(sum(x) < tmax) x <- c(x, rexp(1, lambda))
  data.frame(t = cumsum(x), N = seq_along(x), run = rep(run, length(x)))
}

现在，我们可以在实际的绘图函数中使用此函数：

plot_poisson <- function(runs, tmax, lambda)
{
  # Creates one data frame for each run, this sticks them all together:
  df <- do.call("rbind", lapply(seq(runs), make_sample_df, tmax, lambda))

  ggplot2::ggplot(df, aes(t, N, group = run)) + 
    geom_step(alpha = 0.25) + 
    labs( title = paste(runs, "runs of Poisson process with lambda", lambda)) +
    theme(legend.position = "none") +
    coord_cartesian(xlim = c(0, tmax))
}

所以你可以这样做：

plot_poisson(runs = 10, tmax = 100, lambda = 0.7)

​

​

plot_poisson(runs = 100, tmax = 100, lambda = 0.7)

​

​