Gaussian noise generator algorithm trouble

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.

• From Sonka & al. ?
– user7657
Oct 21, 2015 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$.

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)
}

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

hist(x)

• Any comment by the downvoter?
– Peter K.
Oct 21, 2015 at 13:46