unbiased random number generator

Question

I am generating the random numbers in a domain by using the following code. When I plot them they look grouped at the right. I could show you my plot but I do not know how to upload it. Basically I associate a some data value to the respective point. May you tell me how can I correct it please? My complete code is

#include <iostream>
#include <cmath>
#include <fstream>
#include <sstream>
#include <string>
#include <cstdlib>
#include <cstdio>
#include <time.h>

using namespace std;
string int2string1( int l );
string int2string2( int m );
int main ()
{
    ofstream outFile;
    ofstream myimp;
    string filename;
    srand((unsigned)time(0));
    int nx = 400;
    int ny = 200;
    int i,j,ix,iy,xm,ym,nimp,nfam[nx][ny];
    float vo,rnd,rr,rad,sig,vimp[nx][ny];
    for (i=0; i<nx; i++)
    {
        for (j=0; j<ny; j++)
        {
            vimp[i][j] = 0.0;
        }
    }
    rad = 5.0;
    xm = 0;
    ym = 0;
    vo = 0.08;
    sig = 4.0;
    myimp.open("imp.dat");
    for(i=1; i<nx-1; i++)
    {
        for(j=1; j<ny-1; j++)
        {
            rnd = (random() %1000 + 1)*1.0/1000.0;
            if(rnd>0.99)
            {
                xm = random() % 398 + 1;              /***1 through 399 ***/
                ym = random() % 198 + 1;              /***1 through 199 ***/
                for(ix=xm-5; ix<=xm+5; ix++)
                {
                    for(iy=ym-5; iy<=ym+5; iy++)
                    {
                        rr = sqrt(pow(ix-xm,2.)+pow(iy-ym,2.));
                        if(rr<=rad)
                        {
                            vimp[ix][iy] = vo*1.6e-19;
                        }
                    }
                }
            }
            myimp<<i<<"\t\t"<<j<<"\t\t"<<xm<<"\t\t"<<ym<<"\t\t"<<nfam[i][j]<<"\t\t"<<vimp[i][j]*6.23e18<<"\n";
        }
    }
    myimp.close();
    return 0;
}

Answer 1

int r = rand() % N;

Does not lead to a uniform distribution¹

Instead, I recommend just using C++ TR1 (or boost) random:

#include <random>

std::mt19937 rng(seed);
std::uniform_int_distribution<int> gen(0, N); // uniform, unbiased

int r = gen(rng);

Or to generate floating point numbers of any kind:

std::uniform_real_distribution<double> gen(-2*PI, +2*PI); // uniform, unbiased
double r = gen(rng);

---

_{¹ Backgound, e.g.: Using rand()}

If you're really stuck with using rand() and N which doesn't evenly divide MAX_RAND, the page has some hints on how to achieve somewhat better integer distributions using other formulae. Note that I'd steer you to André Caron's answer instead.

Answer 2

Basically, the rand() % N expression introduces a bias if RAND_MAX is not a multiple of N. It projects numbers in [0,RAND_MAX] onto the range [0,N] with a non-uniform fashion.

Suppose that RAND_MAX=4 and N=2. Then, there are 3 numbers that produce 0 (0, 2 and 4) and 2 numbers that produce 1 (1 and 3). Thus, you have 60% change of getting 0 and 40% chance of getting 1.

The correct way to implement an unbiased projection from [0,RAND_MAX] onto [0,N] is by invoking rand() repeatedly until the random value is in the desired interval. See the documentation for Random.nextInt() in Java (Credits to Oli Charlesworth for the link).

Assuming that, for sheer execution speed, you want to avoid calling rand() multiple times, the way to generate the least bias possible is to use an intermediate double number, such as:

double myrand ()
{
    return double(rand()) / double(RAND_MAX);
}

int myrand ( int max )
{
    return int(myrand() * double(max));
}

Edit: Here's a simple class that will project outputs of the rand() function in to a range [0,N] with no less bias than rand().

class BoundedRandom
{
    const int m_maximum;
    const int m_divisor;
public:
    BoundedRandom ( int maximum )
        : m_maximum(maximum),
          m_divisor(RAND_MAX/(maximum+1))
    {}

    int operator() ()
    {
        int result = rand() / m_divisor;
        while ( result > m_maximum ) {
            result = rand() / m_divisor;
        }
        return (result);
    }
};

Caution: not tested or debugged.

You can use this generator like this:

BoundedRandom randomx(398);
BoundedRandom randomy(198);
// ...
 xm = randomx() + 1; // 1 through 399
 ym = randomy() + 1; // 1 through 199

Answer 3

C++11 introduced random number generators that will almost certainly do a better job than rand. Here's an example using the mersenne twister algorithm, but there are others to choose from depending on the characteristics you need.

// [1, 399]
auto random_int = std::bind(std::uniform_int_distribution<int>(1,399),std::mt19937()); 

// [0, 1.0)
auto random_double = std::bind(std::uniform_real_distribution<double>(0.0,1.0),std::mt19937());

Answer 4

rand() % 398 + 1; /*** 1 through 399 ***/ will generate numbers 1 through 398, because `rand() % 398 will be 0-397 (398 % 398 is 0). Same for the next line.

As an aside, note that using pow with a constant power of two can cost an order of magnitude more CPU than simply writing out the multiplication and should usually be avoided.

Also since you only use rnd in a single comparison against the constant 0.99 you should just work that instead as integer math, as the conversion and comparison in floating point will cost rather more than just doing an integer comparison. For example instead the two lines (rnd = and if) use if((rand() % 100) == 0) which while slightly biased should accurately indicate your intention.

Answer 5

I've read through your code and cannot find anything to bias it to the right. If anything, there's a statistical bias to the left (the left two thirds should be favored by 0.6% ish).