mercredi 10 novembre 2021

Reproducibility of a poisson noise realisation over a sequence

Here is my problem : I have a sequence of float int values (it's a simulated astronomical image but that is not important). I need to get a poisson realisation of this sequence. Then I want to change a few values/pixels and get again a poisson realisation of this new sequence. I need that the unchanged pixel get the same noise.

import numpy
from numpy.random import default_rng,randint
import copy
import pylab

# Create the two lists (with only few pixels changing)
a0 = numpy.random.randint(low=0, high=1e6, size=10000).astype('float')
a = numpy.random.randint(low=0, high=30, size=10000).astype('float')
b = copy.copy(a)
b[a0>9.9e5] += 1e6

# Run poisson on them
gen1 = default_rng(seed=1234)
gen2 = default_rng(seed=1234)
resa = gen1.poisson(a)
resb = gen2.poisson(b)

# Plot the difference (reshape to 2D image to better see the differences)
pylab.imshow(resa.reshape((100,100))-resb.reshape((100,100)), vmin=-10,vmax=10)

With the numpy implementation of poisson (numpy.random.poisson), even if the Seed is the same, the sequence will differ even for the pixels that are identical. The reason is that the Poisson noise is calculated in a conditional loop, which gets a new random value and returns when this value is above some threshold based on lambda. Therefore, after two different pixels, the default_rng (or RandomState) will be different in the two sequences and will give different results even for identical pixels.

Is there an implementation of poisson that would garanty the same number of draw on the RandomState generator whatever the signal level ? Or do I have to tweak the code to do it myself (by looping of draws to always reach a predefined loop iterations) ?




Aucun commentaire:

Enregistrer un commentaire