4

I am trying generate random variates by trying to generate two standard normal variates r1, r2, by using polar coordinates along with a mean and sigma value. However when I run my code, I keep getting a "-nan(ind)" as my output.

What am I doing wrong here? The code is as follows:

static double saveNormal;
static int NumNormals = 0;
static double PI = 3.1415927;

double fRand(double fMin, double fMax)
{
    double f = (double)rand() / RAND_MAX;
    return fMin + f * (fMax - fMin);
}

static double normal(double r, double mean, double sigma) {
    double returnNormal;
    if (NumNormals == 0) {
        //to get next double value
        double r1 = fRand(0, 20);
        double r2 = fRand(0, 20);
        returnNormal = sqrt(-2 * log(r1)) * cos(2 * PI*r2);
        saveNormal = sqrt(-2 * log(r1)) * sin(2 * PI*r2);
    }
    else {
        NumNormals = 0;
        returnNormal = saveNormal;
    }
    return returnNormal*sigma + mean;
}
random_x_y_z
  • 119
  • 8

2 Answers2

3

So, you're using the Box–Muller method to pseudo randomly sample a normal random variate. For this transform to work, r1 and r2 must be uniformly distributed independent variates in [0,1].

Instead, your r1/r2 are [0,20] supported, resulting in a negative sqrt argument when >1, this will give you nans. Replace with

double r1 = fRand(0, 1);
double r2 = fRand(0, 1);

Moreover, you should use C++11 <random> for better pseudorandom number generation; as of now, your fRand has poor quality due to rand()-to-double conversion and possible spurious correlations between adjacent calls. Moreover, your function lacks some basic error checking and badly depends on global variables and is inherently thread unsafe.

FYI, this is what a C++11 version might look like

#include <random>
#include <iostream>

int main()
{
  auto engine = std::default_random_engine{ std::random_device{}() };
  auto variate = std::normal_distribution<>{ /*mean*/0., /*stddev*/ 1. };

  while(true) // a lot of normal samples ...
    std::cout << variate(engine) << std::endl;
}
Arne Vogel
  • 6,346
  • 2
  • 18
  • 31
Massimiliano Janes
  • 5,524
  • 1
  • 10
  • 22
  • "your fRand has poor quality due to rand()-to-double conversion" – Indeed, everyone needs to stop using ``rand()`` and forget it ever existed. :-) See Stephan T. Lavavej's excellent talk from GoingNative 2013: [rand() Considered Harmful](https://channel9.msdn.com/Events/GoingNative/2013/rand-Considered-Harmful) – Arne Vogel Nov 15 '17 at 11:27
1

r1 can be zero, making log(r1) undefined.


furthermore, don't use rand() except when you need your numbers to look random to a human in a hurry. Use <random> instead

peterchen
  • 40,917
  • 20
  • 104
  • 186
  • Good point about the singularity at zero! log(r1) is only greater than zero if r1 > 1, though. (Which it might well be in the posted code.) Maybe you mistyped? – Arne Vogel Nov 20 '17 at 11:42
  • @ArneVogel: d'oh, you are right, I removed that part. – peterchen Nov 20 '17 at 15:35