使用 R 和 ggplot2 模拟泊松过程答案

【问题标题】：Simulation Poisson Process using R and ggplot2使用 R 和 ggplot2 模拟泊松过程
【发布时间】：2020-08-15 22:40:29
【问题描述】：

使用速率 lambda = 0.7 的泊松过程模拟。显示泊松过程的样本运行，纵轴为 N(t)，横轴为时间 t。模拟来自范围 t[0: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))

我该怎么做？

非常感谢！

【问题讨论】：

一些简单的东西，在您的第一个matplot() 中，您使用x 和y，它们是您尚未定义的对象。在第二个matplot() 中，这可能是一个错字，当您实际定义的函数称为Process_poisson() 时，您使用Proceso_poisson()。第三，N_simulations() 不返回任何内容，您需要在函数末尾添加return(s_poisson)。
@OTStats 在大多数情况下都带来了很多好处，虽然我同意您应该将 return(s_poisson) 添加到 N_simulations() 函数的末尾，但事实并非如此必要的;您已经隐式返回s_poisson，但不可见。 OTStats，您可以通过查看str(N_simulations(100, t, lambda)) 来确认这一点
@duckmayr 说得好，感谢您提及

标签： r simulation montecarlo poisson

【解决方案1】：

我认为你可以让这更简单。这是一个ggplot解决方案。

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

当达到目标时，该函数返回向量的累积和，它表示适当 lambda 的泊松过程的到达时间。请注意，为了使绘图更容易，它实际上返回一个包含累积时间、累积计数和分组变量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)

【讨论】：