Ручное моделирование пуассоновского процесса в R

Question

Следующая проблема говорит нам о необходимости поэтапного создания пуассоновского процесса из ρ (время прибытия) и τ (время прибытия).

Один из теоретических результатов, представленных в лекциях, дает следующий прямой метод для моделирования пуассоновского процесса:

• Пусть τ ₀ = 0.
• Генерация iid экспоненциальных случайных величинρ1, ρ2,.,..
• Пусть τ _n = ρ ₁ +.,,+ ρ _n для n = 1, 2,.,,.
• Для каждого k = 0, 1,.,., пусть N _t = k для τ _k ≤ t <τ <sub>k + 1 .

1. Используя этот метод, сгенерируйтереализация пуассоновского процесса (N _t ) _t с λ = 0,5 на интервале [0, 20].
2. Генерация 10000 реализаций пуассоновского процесса (N _t ) _t с λ = 0,5 и используйте ваши результаты для оценки E (N _t ) и Var (N _t ).Сравните оценки с теоретическими значениями.

Мое попытанное решение:

Сначала я сгенерировал значения ρ используя rexp() функцию в R.

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

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

    return (vec)
}

, тогда я создал τ с путем последовательного суммирования ρ с.

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):

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)
}

Теперь, когда у вас есть даос, мы знаем, в какое время происходит каждое событие в интервале времени (0, maxTime). Это приводит нас к генерации ассоциированного Пуассоновского процесса путем нахождения значения Nt для каждого t в интервале времени:

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

Генерирует значение Пуассоновского процесса в каждое целое время интервала. Я подозреваю, что часть вашей проблемы пыталась поместить значения Дао на ось Y вместо количества событий, которые уже произошли. Следующий код помог мне создать случайную лестницу, похожую на ваш пример.

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

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