I stumbled upon the following algorithm for Gaussian noise generation and I just can't figure out how this is supposed to work at all. The algorithm is as follows:

Gaussian noise algorithm

First of all: $j = \arg\min_i (q_1 - p[i])$ does not depend on $q_1$ and therefore yields the same result every time - I assume $\arg\min_i (|q_1 - p[i]|)$ is what's really meant here. But then there's more: I can safely assume that $q_1$ is distributed linearly within $[0,1]$, otherwise I could just scale that, round afterwards and wouldn't have so much hassle. In this case I got another problem: As $q_1$ is linearly distributed, choosing $j$ near to zero (which is what $\mu=0$ would implicate) is not particularly likely, so I suspect the resulting distribution is not a normal distribution.

This algorithm reminds me of inverse transform sampling ( http://en.wikipedia.org/wiki/Inverse_transform_sampling); however it does not appear to be a variant of it. I'm a bit confused here because that's somewhat easy to compute and if p was a cumulative distribution here, I'd understand this one right away.

So basically I'm stuck with understanding this one and looking forward to any hints as to where I'm wrong.

Thank you folks in advance!

  • $\begingroup$ From Sonka & al. ? $\endgroup$ – Yves Daoust Oct 21 '15 at 7:24

I agree with you, this looks like a flawed description of an inverse transform generator.

  • The $p[i]$ shouldn't be like described, but should be the prefix sum of these, corresponding to the cumulative distribution.

  • The $\arg\min$ should indeed be taken on the absolute value (or allow positive values only).

  • As the pdf is computed on the positive $i$ only, the range of the uniform drawing should be $[0,\frac12]$ (or the $p[i]$ scaled by $2$.)

(When the $p[i]$ are not cumulated, they are in $[0,\frac1{\sigma\sqrt2\pi}]$, and the range $[0,1]$ doesn't make sense.)

For better efficiency, you can precompute a lookup-table of the $i$ indexes for sufficiently dense $q_1$.

| improve this answer | |

Yes, the algorithm is distinctly poor. I've made a quick attempt to implement it (as stated) in R, and this is the histogram of the resulting values. Please let me know if there's a problem with it (as the algorithm is stated, not trying to fix the algorithm).

It looks no different from a uniformly distributed random variable.

enter image description here

R Code Below


q_21675 <- function()
  sigma_i <- 3
  G <- 256
  p <- rep(0,G)
  for (i in 1:G)
    p[i] <- 1/(sigma_i*sqrt(2*pi))*exp(-i*i/(2*sigma_i*sigma_i))
  q1 <- runif(1,0,1)
  diff <- rep(0,G)
  for (j in 1:G)
    diff[j] <- q1-p[j]
  j <- which.min(diff)
  q2 <- runif(1,-1,1)
  addition <- q2*j

x <- rep(0,1000)
for (t in 1:1000)
  x[t] <- q_21675()

| improve this answer | |
  • $\begingroup$ Any comment by the downvoter? $\endgroup$ – Peter K. Oct 21 '15 at 13:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.