使用 R 和 ggplot2 模拟泊松过程

Question

使用速率 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))

我该怎么做？

非常感谢！

Answer 1

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

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

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

使用 R 和 ggplot2 模拟泊松过程

Simulation Poisson Process using R and ggplot2

r

simulation

montecarlo

poisson