Моделирование пуассоновского процесса с использованием R и ggplot2

Question

Использование имитации пуассоновского процесса со скоростью лямбда = 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.

Сначала создайте функцию, которая будет имитировать процесс Пуассона, беря выборки, взятые из экспоненциального распределения, с соответствующей лямбдой. В этом примере я использовал while l oop, который начинается с вектора x, чей первый элемент равен 0. Функция увеличивает этот вектор, добавляя случайные выборки, пока его сумма не достигнет целевой продолжительности 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)