Gaussian noise generator algorithm trouble

Question

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:

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!

Answer 1

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$.

Answer 2

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.

---

R Code Below

#21675


q_21675 <- function()
{
  #1 
  sigma_i <- 3
  #2
  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))
  }
  #3 
  q1 <- runif(1,0,1)
  diff <- rep(0,G)
  for (j in 1:G)
  {
    diff[j] <- q1-p[j]
  }
  j <- which.min(diff)
  #4
  q2 <- runif(1,-1,1)
  addition <- q2*j
  return(addition)
}


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

hist(x)