diff options
author | Thomas White <taw@physics.org> | 2009-12-04 17:12:10 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2009-12-04 17:12:10 +0100 |
commit | 8ca26b134904a1d998a713e4215f1a50a5613678 (patch) | |
tree | c76527c8874d51d88e79c8b9dbc5abb01428723e /src | |
parent | b7741b3c35045c1d1d4b0cab3aa2c2b9157247e4 (diff) |
Fix serious numerical problem with Poisson noise generation
Diffstat (limited to 'src')
-rw-r--r-- | src/utils.c | 28 |
1 files changed, 27 insertions, 1 deletions
diff --git a/src/utils.c b/src/utils.c index 5942299d..d3499463 100644 --- a/src/utils.c +++ b/src/utils.c @@ -74,6 +74,28 @@ void progress_bar(int val, int total, const char *text) } +static int fake_poisson_noise(double expected) +{ + double x1, x2, w; + double noise, rf; + + do { + + x1 = 2.0 * ((double)random()/RAND_MAX) - 1.0; + x2 = 2.0 * ((double)random()/RAND_MAX) - 1.0; + w = pow(x1, 2.0) + pow(x2, 2.0); + + } while ( w >= 1.0 ); + + w = sqrt((-2.0*log(w))/w); + noise = w * x1; + + rf = expected + noise*sqrt(expected); + + return (int)rf; +} + + int poisson_noise(double expected) { double L; @@ -82,12 +104,16 @@ int poisson_noise(double expected) L = exp(-expected); + /* For large values of the mean, we get big problems with arithmetic. + * In such cases, fall back on a Gaussian with the right variance. */ + if ( !isnormal(L) ) return fake_poisson_noise(expected); + do { double r; k++; - r = (double)random()/RAND_MAX; + r = (double)random()/(double)RAND_MAX; p *= r; } while ( p > L ); |