mercredi 1 avril 2020

Better implementation of simulating a non-homogeneous Poisson process in Python

I have implemented the thinning approach to simulate a non-homogeneous Poisson process. See Wiki Link and the algorithm can be found on page 92 of this PDF notes. My question is if my implementation is efficient in terms of Python realization.

As a side note, this question is similar to the one in this old post, where the answer also pointed out that thinning is better from an algorithm perspective, at least in my understanding.

Here is an MWE. Please comment if there is any way to improve the implementation. Thank you!

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




Aucun commentaire:

Enregistrer un commentaire