I'm trying to simulate a system of reactions over time. In order to do this I have to multiply the value of probability of a reaction occurring with a pre-calculated time step that it can occur in, save this result in new variable and use the new variable to sample from the poisson distribution.
This is a snippet of my code:
lam = (evaluate_propensity*delta_t)
rxn_vector = np.random.poisson(lam) # probability of a reaction firing in the given time period
I've written a function to calculate the value of delta_t based on system specific parameters, the value calculated is very small 0.00014970194372884217 and I think this is having an impact on the np.random.poisson function.
The evaluate_propensity variable is an array that details the probability of a reaction occurring based on the number of molecules in the system and the ratios between molecules in a reaction. This is calculated dynamically and changes after each iteration as the molecule numbers change, but the values for the first iteration are:
evaluate_propensity = np.array([1.0, 0.002, 0.0, 0.0])
The documentation states that lam must be >= 0 and mine is (just) but rxn_vector just always returns an array of zeros.
rxn_vector = [0 0 0 0]
I know that the last two elements of the array will evaluate to zero. But didn't think that the first two would as well. Is there a way to make it more sensitive or amplify my results somehow or am I doing something wrong?
Cheers
The probability to draw a non-zero number for lambda = 1.5e-4 is tiny, it is P(k>0) = 1 - P(k=0) = 1.5e-4. On average you'd have to draw a lot more than four samples to get a non-zero value, 1 / 1.5e-4 = 6667 samples for propensity = 1. For smaller values the number of necessary samples is obviously even larger.
You can confirm this with scipy.stats
from scipy.stats import poisson
pdist = poisson(1.5e-4)
prob = 1 - pdist.pmf(0)
print(prob) # 0.00014998875056249084