renewal_process_dead_time.py

"""
directly simulate N(t) for a renewal process
Include dead time effect
"""

from random import expovariate,uniform
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

from joblib import Parallel, delayed
import multiprocessing

def Nalt_runner(i):
    print i
    pc = Nalt()
    return pc

def Nalt(t_obs = 0.001, lam=1.0e9, taud = 100e-7, sigd=10e-7):
    last_photon_time = 0
    last_count_time = -np.inf
    photon_count = 0


    while last_photon_time <= t_obs:
        last_photon_time = last_photon_time+expovariate(lam)

        fuzzy_dead_time = uniform(taud-np.sqrt(3)*sigd, taud+np.sqrt(3)*sigd)

        if last_photon_time >= last_count_time+fuzzy_dead_time:
            last_count_time = last_photon_time
            photon_count += 1

    return photon_count


def N(t_obs=0.001, lam=1.0e9, taud=100e-7, sigd=10e-7):
    Nt = 0
    t_cur = 0.0
    while(t_cur<=t_obs):
        t_cur += expovariate(lam) + uniform(taud-np.sqrt(3)*sigd,taud+np.sqrt(3)*sigd)
        Nt += 1

    return Nt

if __name__=='__main__':
    Nsim = 1000

    # Method 1
    Nt_nounc = np.zeros(Nsim)
    for i in tqdm(range(Nsim)):
        Nt_nounc[i] = N()
    
    # Method 2
    Nalt_nounc = np.array( Parallel(n_jobs=35)(delayed(Nalt_runner)(i) for i in range(Nsim)) )
    
    print np.mean(Nt_nounc), np.mean(Nalt_nounc)
    print np.std(Nt_nounc), np.std(Nalt_nounc)