#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu May  4 14:37:10 2023

@author: robertwebber
"""

import numpy as np
import matplotlib.pyplot as plt

# pretty plots
font = {'family' : 'sans-serif',
        'weight' : 'regular',
        'size'   : 18}
plt.rc('font', **font)
plt.rcParams['axes.linewidth'] = 1.5

#%%%

T = int(1e6)
lam = 3
states = np.zeros(T, dtype = int)
now = 0
for i in range(T):
    if i % 1000 == 0:
        print(i)
    states[i] = now
    if np.random.uniform() < .5:
        if np.random.uniform() < now / lam:
            now = now - 1
    else:
        if np.random.uniform() <= lam / (now + 1):
            now = now + 1

plt.hist(states, bins = int(20*lam), color = 'red')
x = np.arange(int(10*lam))
y = np.zeros(x.size)
now = np.exp(-lam)
for i in range(x.size):
    y[i] = now
    now = now * lam / (i + 1)
plt.plot(x, y*T, color = 'black', linestyle = '--')
plt.xlim([-.5, 5*lam])