Лучшая реализация симуляции неоднородного пуассоновского процесса в Python - PullRequest
Новая
       105

Лучшая реализация симуляции неоднородного пуассоновского процесса в Python

0 голосов
/ 02 апреля 2020

Я реализовал метод прореживания , чтобы моделировать неоднородный пуассоновский процесс. См. Wiki Link , и алгоритм можно найти на странице 92 этого PDF примечания . Мой вопрос заключается в том, эффективна ли моя реализация с точки зрения реализации Python.

В качестве примечания, этот вопрос похож на тот, что был в в этом старом посте где ответ также указал, что истончение лучше с точки зрения алгоритма 1016 *, по крайней мере, в моем понимании.

Вот MWE. Пожалуйста, прокомментируйте, если есть какой-либо способ улучшить реализацию. Спасибо! 

import numpy as np
import matplotlib.pyplot as plt


mytime = np.arange(0,24,0.001)
lam1 = 1+0.6*np.sin(mytime)
lam2 = 1-0.6*np.sin(mytime)

max_lam = 1.8
P = []
P.append(lam1/max_lam)
P.append(lam2/max_lam)

np.random.seed(0)

def nhpp_next_arrival(t,Lambda, h_lam):
    # current starting time t
    # rate handle h_lam
    # Lambda a majorizing function/constant

    U = np.random.uniform()
    V = np.random.uniform()

    t = t - np.log(U)/Lambda


    while V > h_lam(t)/Lambda:

        t = t - np.log(U)/Lambda
        U = np.random.uniform()
        V = np.random.uniform()

    return t


def nhpp(h_lam,max_lam,T_cap=24):
# lam is a function 

    t = 0
    # n = 0
    T = [0]

    while t <= T_cap:
        t = T[-1]

        T.append(nhpp_next_arrival(t=t, Lambda=max_lam, h_lam=h_lam))

    return T

if __name__ == '__main__':
    n = 100
    h_lam = lambda t: (1+0.6*np.sin(t))*n

    T = nhpp(h_lam, max_lam=n*max_lam)
    plt.plot(T,np.arange(0,len(T)), label='sim')
    plt.plot(mytime, np.cumsum(h_lam(mytime))*0.001,label='ref')
    plt.title('scale = %s' % str(n))
    plt.legend()
...