Simulace nehomogenního Poissonova procesu
V této části si ukážeme praktický způsob, jak generovat trajektorii Poissonova procesu , s danou intenzitou
měnící se v čase. Protože se intenzita v čase mění, nelze přímo použít způsob, kdy se vygenerují délky intervalů z exponenciálního rozdělení a jejich kumulativní součty pak určují okamžiky událostí v Poissonově procesu. Jednou z možností je využití definičních vlastností (i*) a (ii*). Na časové ose postupujeme po velmi malých krocích délky
. V každém kroku, v čase
, je pravděpodobnost příchodu události rovna
a považujeme ji za konstantní během krátké doby
. Vygenerujeme realizaci z alternativního rozdělení s pravděpodobností úspěchu
,
a v případě nastání úspěchu (tj. realizace je rovna 1) zapíšeme nový okamžik příchodu události. Tento postup je ale zdlouhavý a nepřesný, z důvodu nutnosti velmi jemného dělení časové osy.
V praxi se používá jiný, přesnější a svým přístupem zajímavý přístup k simulaci trajektorie nehomogenního Poissonova procesu. Pro proměnnou intenzitu určíme maximální hodnotu funkce na zkoumaném intervalu
a označíme ji
; např. pro periodickou intenzitu (4) je
. Následně, způsobem popsaným v části o homogenním Poissonově procesu, vygenerujeme trajektorii Poissonova procesu pro konstantní intenzitu rovnou
. Označme takto vygenerované okamžiky příchodů událostí
. Tyto časy považujeme za potenciální okamžiky událostí v nehomogenním procesu, přičemž o tom, jestli v čase
skutečně nastala událost, rozhodneme na základě realizace alternativní náhodné veličiny s pravděpodobností úspěchu rovnou podílu
. Celý algoritmus lze shrnout následovně:
- nalezneme
- vygenerujeme vzorek
okamžiků příchodů událostí v homogenním Poissonově procesu s intenzitou
;
- pro každý okamžik
,
, vygenerujeme realizaci
z alternativního rozdělení pravděpodobnosti s pravděpodobností úspěchu
;
- výsledný vzorek okamžiků příchodů událostí nehomogenního Poissonova procesu sestává právě z těch časů
.
Nyní ukážeme, jak trajektorii vygenerovat ve statistickém softwaru R [3] pro periodickou intenzitu (4). Využijeme přitom postup z kapitoly Statistické metody pro Poissonův proces. Po úvodním uložení hodnot do proměnných spočítáme maximální intenzitu Lambda a vygenerujeme okamžiky S homogenního Poissonova procesu s touto intenzitou. Spočítáme hodnoty proměnné intenzity lambda v těchto časech S a následně pomocí podílu pravděpodobnosti p. Ty využijeme jako parametry funkce rbinom pro generování náhodného výběru Z z alternativního rozdělení (které je speciálním případem binomického s parametrem size=1). Jako okamžiky událostí v nehomogenním Poissonově procesu Sn vybereme pomocí indexování S[ ] a funkce which právě ty z časů S, pro něž je odpovídající Z rovno jedné (Z == 1). Dále již analogicky jako u homogenního Poissonova procesu doplníme hodnoty čítacího procesu N (od nuly) a vytvoříme schodovitou funkci X odpovídající trajektorii procesu. Graf takto vygenerované trajektorie nehomogenního Poissonova procesu je na Obr.1.
K <- 150
a <- 1.5
b <- 1.5
omega <- 2 * pi / 24
phi <- -pi / 2
Lambda <- a + b
T <- rexp (K, rate=Lambda)
S <- cumsum (T)
lambda <- a + b * sin (phi + omega * S)
p <- lambda / Lambda
Z <- rbinom (K, size=1, prob=p)
Sn <- S[ which (Z == 1) ]
N <- seq (0, length (Sn), by=1)
X <- stepfun (Sn, N)
plot (X, verticals=FALSE, pch=20, col=2, lwd=1.5,
xlab="t", ylab="N(t)", main="")
