2

I am attempting to generate a random distribution that follows an upside-down gaussian distribution, shifted uo so that it is still in range(0,1). I need to do this with as few special functions as possible and can only use a flat random number generator.

I am able to generate according to a Gaussian by putting the flat random numbers through the inverse Gaussian CDF. This works and gives me the gaussian dist that I would expect. In python, this looks like this:

def InverseCDF(x, mu, sigma):
  return mu + sigma * special.erfinv(2*x - 1)

Now when I am trying to generate a distribution that follows 1-e^(-x^2), I believe the inverse CDF of this function is the same as for the regular gaussian with the argument of the inverse error function now 2*p + 1. So it would look like below:

   def InverseCDF(x, mu, sigma):
      return mu + sigma * special.erfinv(2*x + 1)

The problem here is that erfinv is only defined from (-1,1) and the argument is now greater than 1. I have tried scaling this and flipping in all sorts of ways, putting negatives almost everywhere I can, and I can never seem to generate a histogram that follows an upside-down gaussian. In most cases, I actually get back a regular gaussian distribution.

Any idea what I'm doing wrong, or any tips on how to generate this upside-down gaussian? Thanks in advance for any help.

SClark
  • 21
  • 2
  • 1
    Interesting question, but I'm a little confused about the definition of the pdf. Let g(x) be the ordinary normal pdf with mean 0 and variance 1. I guess you are looking at max(g) - g(x) where max(g) is the maximum value of the bump (which I guess is 1/sqrt(2 pi)). But max(g) - g(x) has infinite area under the curve -- are you limiting the range of x in some way? – Robert Dodier May 08 '20 at 15:48
  • Yes I should have stated that there is a limit on x. The limit is determined by my specific problem but we can take it to be (0,1) for this toy example. – SClark May 08 '20 at 16:46
  • Is your truncated normal centered on the range, or is the mode at zero? Are you willing to use acceptance/rejection techniques? – pjs May 08 '20 at 18:07
  • Does this answer your question? [Inverse normal random number generation in python?](https://stackoverflow.com/questions/62899379/inverse-normal-random-number-generation-in-python) – Peter O. Oct 14 '20 at 10:41

1 Answers1

1

OK, with x between 0 and 1, I get this for the cdf: 

-(sqrt(%pi)*(sqrt(2)*sigma*erf((sqrt(2)*x-sqrt(2)*mu)/(2*sigma))
            +sqrt(2)*erf(mu/(sqrt(2)*sigma))*sigma)
 -2*x)
 /(sqrt(%pi)*(sqrt(2)*erf((sqrt(2)*mu-sqrt(2))/(2*sigma))
             -sqrt(2)*erf(mu/(sqrt(2)*sigma)))*sigma
  +2)

Maybe some algebra will make it possible to figure out a formula for the inverse, if not, I guess a numerical root search will work. I guess it will be simpler for specific values of mu and sigma.

I did that with Maxima (http://maxima.sourceforge.net), by constructing the pdf and integrating it. Plotting the expression above yields a plausible picture.

Robert Dodier
  • 16,905
  • 2
  • 31
  • 48
  • This looks promising, thanks for the help. I don't think the inverse will be too hard to get to, but it's a bit too complex/hard to read write now. We can say mu=0 which will simplify things quite a bit. In that case the second terms go away in both numerator and denominator – SClark May 08 '20 at 17:36
  • With mu=0 this gives: -(sqrt(π) (sqrt(2) σ erf((sqrt(2) x)/(2 σ))) - 2 x) / (sqrt(π)×erf(-sqrt(2))/(2 σ) + 2) – SClark May 08 '20 at 17:44
  • 1
    I looked at it for a few minutes but I don't see an easy solution. I think a numerical root finder (even just bisection) is a workable way to solve it. That would work for any values of mu and sigma (and also any range, I'm pretty sure). – Robert Dodier May 08 '20 at 19:47