I'm trying to replicate some published research ( DOI: 10.1021/ja908004w ), and part of that involves generating a list of random Poisson-distributed numbers with additional sine weighting.
The supporting information for the paper provides their C source code, and I'm trying to replicate/rewrite it in Python. However, I'm stuck with some for/while/if loops.
Essentially, I need to generate a list of unique numbers between 0 and x which are poisson distributed with this formula
where theta varies linearly from 0 to pi/2 for our x random variables
I've written this much:
import numpy as np
import math
seed = 42
np.random.seed(seed)
per = 0.25 #percent of NUS
z = 1024 # total size of indirect dimension
p = int(per*z) #sampled size of indirect dimension
def sinweightedpoisson():
ld = z/p
adj = 2*(ld-1) #initial guess
samples = []
for i in range(p):
lmbd = adj*math.sin((i+0.5)/(z+1)*math.pi/2)
s = np.random.poisson(lam=lmbd,size=1)
while (s in samples):
s = np.random.poisson(lam=lmbd,size=1)
else:
samples.append(int(s))
However, obviously this gets stuck in an infinite loop because of requirement for the numbers to be unique.
The solution is to modify the adj variable slightly until we get the right number of numbers.
i.e.
if len(samples) > p:
adj = adj *1.02
elif len(samples) < p:
adj = adj /1.02
and to repeat the loop until len(samples) == p.
Any ideas?
Aucun commentaire:
Enregistrer un commentaire