Как мне генерировать пуассоновский процесс?

Question

Оригинальный вопрос:

Я хочу сгенерировать процесс Пуассона. Если число поступлений по времени t равно N (t) и у меня есть распределение Пуассона с параметром λ , как я могу сгенерировать N (t ) ? Как бы я сделал это в C ++?

Пояснение:

Изначально я хотел сгенерировать процесс, используя распределение Пуассона. Но я был смущен тем, какой параметр из процесса мне нужен; Я думал, что смогу использовать N (t) , но это говорит мне, сколько прибытий произошло на интервале (0, t] , что было не то, что я хотел. Итак, я думал, что я мог бы использовать N (t2) -N (t1) , чтобы получить количество поступлений на интервале [t1, t2] . Так как N (t) ~ Пуассон (tx λ) Я мог бы использовать Пуассона (t2 x λ) -Пуассона (t1 x λ) , но я не хочу количество прибытий в интервале.

Скорее, я хочу сгенерировать явное время прибытия.

Я мог бы сделать это, сделав интервал [t2, t1] достаточно маленьким, чтобы у каждого интервала было только одно прибытие (что происходит как | t2-t1 | -> 0 ).

Answer 1

Если у вас есть процесс Пуассона с параметром скорости L (что означает, что в долгосрочной перспективе существует L прибытий в секунду), то время между прибылями экспоненциально распределяется со средним значением 1 / L. Таким образом, PDF это f (t) = -L * exp (-Lt), а CDF это F (t) = Prob (T

Предполагая, что в используемом вами языке есть функция (назовем ее rand()) для генерации случайных чисел, равномерно распределенных между 0 и 1, метод обратного CDF сводится к вычислению:

-log(rand()) / L

Поскольку python предоставляет функцию для генерации экспоненциально распределенных случайных чисел, вы можете смоделировать первые 10 событий в пуассоновском процессе со средней скоростью 15 поступлений в секунду, например:

import random
for i in range(1,10):
   print random.expovariate(15)

Обратите внимание, что это будет генерировать * меж * время прибытия. Если вам нужно время прибытия, вам нужно будет продолжать перемещать переменную времени вперед, как это:

import random
t= 0
for i in range(1,10):
   t+= random.expovariate(15)
   print t

Answer 2

Вот пример кода для генерации сэмплов Пуассона с использованием C ++ TR1 .

Если вы хотите процесс Пуассона , времена между поступлениями экспоненциально распределены, и экспоненциальные значения могут быть сгенерированы тривиально с помощью обратного метода CDF: -k * log (u) где u - равномерная случайная величина и k - среднее экспоненты.

Answer 3

Я был бы очень осторожен с использованием обратного CDF и прокачкой через него равномерного случайного числа. Проблема здесь состоит в том, что часто обратный CDF является численно нестабильным или функции, которые его производят, могут давать нежелательные колебания вблизи концов интервала. По этой причине я бы порекомендовал что-то вроде метода отбраковки, используемого в «Числовых рецептах в Си». См. Функцию poidev, приведенную в разделе 7.3 NRC: http://www.nrbook.com/a/bookcpdf/c7-3.pdf

Answer 4

Если вы используете python, вы можете использовать random.expovariate (rate) для генерации времени прибытия с частотой событий за интервал времени

Answer 5

Чтобы выбрать выборку из распределения, вам необходимо вычислить обратную интегральную функцию распределения (CDF). Сначала вы выбираете случайное число равномерно на реальном интервале [0, 1], а затем берете обратный CDF этого значения.

Answer 6

Генерирование времени прибытия через процесс Пуассона не означает использование распределения Пуассона. Это делается путем создания экспоненциального распределения на основе коэффициента пуассоновского прихода лямда.

Короче говоря, вам нужно сгенерировать экспоненциальное распределение со средним значением = 1 / лямда, см. Следующий пример:

#include <iostream>
#include <iterator>
#include <random>

int
main ()
{
 // seed the RNG
 std::random_device rd; // uniformly-distributed integer random number generator
 std::mt19937 rng (rd ()); // mt19937: Pseudo-random number generation

 double averageArrival = 15;
 double lamda = 1 / averageArrival;
 std::exponential_distribution<double> exp (lamda);

double sumArrivalTimes=0;
double newArrivalTime;


 for (int i = 0; i < 10; ++i)
  {
   newArrivalTime=  exp.operator() (rng); // generates the next random number in the distribution 
   sumArrivalTimes  = sumArrivalTimes + newArrivalTime;  
   std::cout << "newArrivalTime:  " << newArrivalTime  << "    ,sumArrivalTimes:  " << sumArrivalTimes << std::endl;  
  }

}

Результат выполнения этого кода:

newArrivalTime:  21.6419    ,sumArrivalTimes:  21.6419
newArrivalTime:  1.64205    ,sumArrivalTimes:  23.2839
newArrivalTime:  8.35292    ,sumArrivalTimes:  31.6368
newArrivalTime:  1.82962    ,sumArrivalTimes:  33.4665
newArrivalTime:  34.7628    ,sumArrivalTimes:  68.2292
newArrivalTime:  26.0752    ,sumArrivalTimes:  94.3045
newArrivalTime:  63.4728    ,sumArrivalTimes:  157.777
newArrivalTime:  3.22149    ,sumArrivalTimes:  160.999
newArrivalTime:  1.64637    ,sumArrivalTimes:  162.645
newArrivalTime:  13.8235    ,sumArrivalTimes:  176.469

Итак, на основании вашего эксперимента вы можете использовать: newArrivalTime или sumArrivalTimes.

ref: http://www.math.wsu.edu/faculty/genz/416/lect/l05-45.pdf

Answer 7

В Python вы можете попробовать код ниже.

Если вы хотите создать 20 случайных показаний за 60 секунд. то есть (20 лямбда)

 def poisson_job_generator():
    rateParameter = 1.0/float(60/20) 
    while True:
        sl = random.expovariate(rateParameter)

Answer 8

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

https://stackoverflow.com/a/15307412/1650437