NHPP jest definiowany w terminach monotonicznej, niemalejącej, prawostronnie ciągłej funkcji \(\displaystyle{ \Lambda(x) }\), ograniczonej na skończonym przedziale.
Liczba punktów na dowolnym przedziale \(\displaystyle{ (0, \ \ x_{0}] }\) ma rozkład Poissona z parametrem \(\displaystyle{ \mu_{0}= \Lambda(x) - \Lambda(x_{0}).}\)
Pochodna prawostronna funkcji \(\displaystyle{ \Lambda(x) }\) jest zwykle oznaczana przez \(\displaystyle{ \lambda(x) }\) i nazywana funkcją intensywności procesu.
Za pomocą funkcji \(\displaystyle{ \lambda(x) }\) realizacje (zdarzenia) NHPP - \(\displaystyle{ N(t) }\) i ich własności opisują równania:
\(\displaystyle{ Pr( \{ N(0) =0\}) = 1, }\)
\(\displaystyle{ Pr( \{N(t_{1}) \cap N(t_{2}) \}) = Pr(\{N(t_{1}\}) \cdot Pr(\{N(t_{2} \}) \ \ \mbox{dla} \ \ t_{1}\neq t_{2} }\) (zdarzenia są niezależne),
\(\displaystyle{ Pr(\{ N(t+s) –N(t)= k \} ) = \frac{ \left( \int_{t}^{t+s} \lambda(x)dx \right )^{k}}{k!} e^{-\int_{t}^{t+s}\lambda(x)dx}. }\)
W praktyce stosuje się pięć różnych metod symulacji NHHP.
Do nich należą:
- metody skalowania przedziałów czasowych jednorodnego procesu Poissona,
- metody indywidualnego generowania przedziałów czasowych,
- metody estymacji statystyki pozycyjnej (porządkowej),
- metody symulacji metodą rozrzedzania (Tainning Simulation),
- metody indywidualnego generowania przedziałów przez rozrzedzanie.
Zastosowano pierwszą z tych metod.
Jej idea wynika ze związku między jednorodnym procesem Poissona (Homogeneus Poisson Process (HPP)), w którym intensywność \(\displaystyle{ \lambda }\) jest stała-niezależna od czasu i niejednorodnym procesem Poissona NHPP.
Jeżeli \(\displaystyle{ N_{1}(t) }\) jest jedną z realizacji HPP, to odpowiadający jej wewnątrz przedziałów rozkład \(\displaystyle{ I_{0} }\) jest rozkładem wykładniczym
\(\displaystyle{ P(I_{0}\geq(t))= \exp(-t) \rightarrow }\)
\(\displaystyle{ \rightarrow \ \ P(I_{0}\geq \Lambda(t))= \exp(-\Lambda(t)) \rightarrow }\)
\(\displaystyle{ \rightarrow P(I_{0}\geq \Lambda^{-1}(t))= \exp(-\Lambda^{-1}(t)) }\)
Niech \(\displaystyle{ I'_{0} }\) oznacza rozkład wewnątrz przedziału NHHP, wtedy \(\displaystyle{ P(I'_{0}\geq t)= exp(-\Lambda^{-1}(t)) }\)
Z równości tych wnioskujemy, że jeżeli
\(\displaystyle{ X'_{1}, X'_{2},...}\) są punktami NHHP, z ciągłą i całkowalną funkcją \(\displaystyle{ \Lambda(x) }\), to wtedy punkty \(\displaystyle{ X_{1}= \Lambda(X'_{1}), X_{2}= \Lambda(X'_{2})...}\) są punktami HPP .
Program get_realization_nhpp w kodzie języka programu R 3.6.0, przedstawia symulację dwóch funkcji - funkcji liniowej postaci \(\displaystyle{ f(t) = 1 + bt }\) i funkcji okresowej intensywności \(\displaystyle{ g(t) = b\cdot\sin(2\pi t).}\)
Dlaczego akurat przeprowadzamy symulację tych funkcji intensywności - funkcji liniowej i funkcji okresowej sinus?
Liniowa funkcja intensywności NHHP służy do modelowania ilości klientów, przybywających w określonym przedziale czasowym do obiektów masowej obsługi takich jak: galerie, hipermarkety, MacDonaldy, dworce ... itp. .
Okresowa funkcja sinus służy do modelowania zjawisk okresowych, takich jak trzęsienia Ziemi, powodzie, wylewy mórz I oceanów, tsunami ...itp.
Kod i wyniki programu get_nhpp_realization dla \(\displaystyle{ t \in [ 0, t_{max}] = [0, 10] }\) i wartości parametru \(\displaystyle{ b}\) należących do zbiorów \(\displaystyle{ \{0.01, 0.1, 1 \} , \ \ \{0.1, 1, 10 \} }\) odpowiednio dla funkcji liniowej i okresowej funkcji sinus.
Kod: Zaznacz cały
> bs <- c(0.01, 0.1, 1)
> get_nhpp_realization <- function(lambda){
+ set.seed(1)
+ t_max <- 10
+ t <- 0
+ s <- 0
+ Lambda <- function(tupper) integrate(f = lambda, lower = 0, upper = tupper)$value
+ Lambda_inv <- function(s){
+ v <- seq(0,t_max+3, length.out = 1000)
+ min(v[Vectorize(Lambda)(v)>=s])
+ }
+ X <- numeric(0)
+ while(t <= t_max){
+ u <- runif(1)
+ s <- s -log(u)
+ t <- Lambda_inv(s)
+ X <- c( X, t)
+ }
+ return(X)
+ }
> b <- bs[1]
> lambda <- function(t) 1 + b*t
> res_1 <- get_nhpp_realization(lambda)
> n_1 <- length(res_1)
> b <- bs[2]
> lambda <- function(t) 1 + b*t
> res_2 <-get_nhpp_realization(lambda)
> n_2 <- length(res_2)
> b <- bs[3]
> lambda <- function(t) 1 + b*t
> res_3 <- get_nhhp_realization(lambda)
> n_3 <- length(res_3)
> get_nhpp_realization(lambda)
[1] 0.9239239 1.3793794 1.6006006 1.6396396 2.1861862 2.2252252
[7] 2.2382382 2.3683684 2.4984985 3.2272272 3.5785786 3.9429429
[13] 4.0210210 4.2162162 4.2552553 4.3853854 4.4504505 4.4504505
[19] 4.6326326 4.6716717 4.6846847 4.9579580 5.0230230 5.3613614
[25] 5.5695696 5.7127127 6.3243243 6.4544545 6.4674675 6.6106106
[31] 6.7147147 6.7797798 6.8708709 7.0790791 7.1051051 7.1571572
[37] 7.1831832 7.4434434 7.4824825 7.5865866 7.6126126 7.6646647
[43] 7.6906907 7.7557558 7.8338338 7.8598599 8.2762763 8.3543544
[49] 8.3803804 8.4194194 8.4974975 8.5235235 8.6016016 8.7447447
[55] 9.0180180 9.2522523 9.3563564 9.4214214 9.4604605 9.5515516
[61] 9.5515516 9.6686687 9.7467467 9.8508509 9.8898899 10.0070070
> bs <- c(1, 10, 20)
> b <- bs[1]
> lambda <- function(t) 1 + b*(1+sin(2*pi*t))
> res_1 <- get_hhp_realization(lambda)
> n_1 <- length(res_1)
> b <- bs[2]
> lambda <- function(t) 1 + b*(1+sin(2*pi*t))
> res_2 <- get_nhpp_realization(lambda)
> n_2 <- length(res_2)
> b <- bs[3]
> lambda <- function(t) 1 + b*(1+sin(2*pi*t))
> res_3 <- get_nhpp_realization(lambda)
> n_3 <- length(res_3)
> get_nhhp_realization(lambda)
[1] 0.06506507 0.09109109 0.11711712 0.11711712 0.15615616 0.15615616
[7] 0.16916917 0.16916917 0.18218218 0.26026026 0.29929930 0.33833834
[13] 0.35135135 0.37737738 0.37737738 0.40340340 0.41641642 0.41641642
[19] 0.44244244 0.45545546 0.45545546 0.52052052 0.54654655 0.92392392
[25] 1.00200200 1.04104104 1.18418418 1.19719720 1.21021021 1.23623624
[31] 1.24924925 1.26226226 1.27527528 1.32732733 1.32732733 1.34034034
[37] 1.34034034 1.40540541 1.41841842 1.44444444 1.45745746 1.47047047
[43] 1.48348348 1.50950951 1.53553554 1.54854855 2.03003003 2.05605606
[49] 2.06906907 2.08208208 2.10810811 2.10810811 2.13413413 2.17317317
[55] 2.23823824 2.29029029 2.31631632 2.34234234 2.35535536 2.36836837
[61] 2.38138138 2.40740741 2.43343343 2.48548549 2.49849850 2.57657658
[67] 2.66766767 2.82382382 3.03203203 3.03203203 3.07107107 3.08408408
[73] 3.11011011 3.14914915 3.16216216 3.16216216 3.17517518 3.20120120
[79] 3.20120120 3.20120120 3.22722723 3.22722723 3.25325325 3.27927928
[85] 3.29229229 3.33133133 3.34434434 3.39639640 3.43543544 3.52652653
[91] 3.64364364 4.02102102 4.04704705 4.04704705 4.06006006 4.06006006
[97] 4.09909910 4.12512513 4.12512513 4.13813814 4.15115115 4.17717718
[103] 4.21621622 4.21621622 4.22922923 4.26826827 4.30730731 4.33333333
[109] 4.33333333 4.34634635 4.34634635 4.35935936 4.38538539 4.41141141
[115] 4.47647648 4.98398398 4.99699700 5.08808809 5.11411411 5.12712713
[121] 5.12712713 5.14014014 5.16616617 5.20520521 5.21821822 5.23123123
[127] 5.25725726 5.29629630 5.32232232 5.33533534 5.36136136 5.42642643
[133] 5.59559560 5.66066066 5.68668669 5.88188188 5.94694695 5.98598599
[139] 5.98598599 6.01201201 6.03803804 6.05105105 6.10310310 6.14214214
[145] 6.15515516 6.16816817 6.22022022 6.22022022 6.27227227 6.28528529
[151] 6.29829830 6.31131131 6.33733734 6.36336336 6.37637638 6.42842843
[157] 6.45445445 6.57157157 6.88388388 7.00100100 7.05305305 7.05305305
[163] 7.09209209 7.09209209 7.09209209 7.11811812 7.19619620 7.22222222
[169] 7.23523524 7.26126126 7.27427427 7.27427427 7.27427427 7.30030030
[175] 7.32632633 7.32632633 7.33933934 7.35235235 7.36536537 7.36536537
[181] 7.40440440 7.45645646 7.45645646 7.48248248 7.49549550 7.58658659
[187] 7.62562563 7.67767768 7.70370370 7.89889890 7.93793794 7.98998999
[193] 8.08108108 8.09409409 8.12012012 8.14614615 8.19819820 8.19819820
[199] 8.23723724 8.23723724 8.27627628 8.30230230 8.32832833 8.35435435
[205] 8.40640641 8.43243243 8.44544545 8.53653654 8.69269269 8.86186186
[211] 8.86186186 9.03103103 9.04404404 9.04404404 9.04404404 9.09609610
[217] 9.10910911 9.10910911 9.10910911 9.13513514 9.17417417 9.22622623
[223] 9.25225225 9.26526527 9.26526527 9.27827828 9.31731732 9.39539540
[229] 9.42142142 9.43443443 9.47347347 9.57757758 9.81181181 9.90290290
[235] 9.98098098 9.99399399 10.02002002
Załącznik 1 - wykres symulacji funkcji liniowej
Załącznik 2 - wykres symulacji funkcji sinus
Wnioski
Własności funkcji intesywności NHHP potwierdza przeprowadzona symulacja. Funkcje te są jest monotonicznie niemalejące. Wartości teoretyczne tych funkcji \(\displaystyle{ \lambda(t) }\) odpowiadają w przybliżeniu wysymulowanym wartością \(\displaystyle{ \lambda_{w}}\). Na przykład wartość liniowej funkcji intensywności dla przedziału czasowego \(\displaystyle{ [1, 2] }\)
\(\displaystyle{ \lambda(t) = \int_{1}^{2} ( 1 + t) dt = \left[ t + \frac{1}{2}t^2\right]_{1}^{2} = 2,5. }\)
\(\displaystyle{ \lambda_{w} = 2.4984985.}\)
Jeśli wygenerujemy NHHP o liniowej funkcji intensywności \(\displaystyle{ \lambda }\), to z wykresu widać, że porządek liniowy został zaburzony i nie można dokładnie uzyskać losowej prostej intensywności.
Fakt ten można formalnie udowodnić, Dowód, czy precyzyjne sformułowanie tego stwierdzenia wykracza poza ramy tego zadania.