Gauss distribution using rand()

Question

i am up to build a Monte Carlo based simulation in C since it is super fast and I was wondering why my code is producing a even distribution.

So first to pick you up, bevor I started coding I imagened the picture of a few balls falling down and on every dot they can move left or right with random distribution. The picture is showen here: https://www.youtube.com/watch?v=PM7z_03o_kk.

But what I got is kind of strange. When I set the scatterpoints to 10 (which is in the code example setted to 100):

 while(j < 100) // Number of abitrary change of direction

I got a picture like a Gauss-Distribution but only even bin´s are contributing to it. When its large enough like shown in the code, every bin got about the same amount of particles. This is the 2D case which will be expanded to the 3D case once it is working as expected.

There are still some Variables in it which are not really neccessary just to avoid any possible mistake I can imagine. I was working with gdb to find the error. When I just run distr() with gdb I figured out that it is producing only even numbers if the example above is setted to 10. When I go to 11 I found out, that bin[0] starts to contribute with a very small amout compared to the rest. I also ran rando() enough times to see if it's really pseudo-random and I found out that it should work.

I was still not able to figure out the mistake, so I really hope that here are enough people smarter than me oO.

Full code:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#define SEED time(NULL)


int rando(void) // Monte Carlo method for choosing weather left or right
{
    double temp=0;

    temp = (double) rand()/RAND_MAX; // pseudo random number between 0-1

    if (temp >= 0.5)
        return 1; // Particle to the right
    if (temp < 0.5)
        return 0; // Particle ot the left
    else
        return EXIT_FAILURE;
}

int distr(void) // Binning particle
{
    int i=10; // Center of bin
    int j=0;
    int k=0;

    while(j < 100) // Number of abitrary change of direction
    {
        k=rando();

        if(k == 1)
        {
            if ( i < 13) // Choose upper bound of bins
                i++;
            j++;
        }
        if(k == 0)
        {
            if (i > 7) // Choose lower bound of bins
                i--;
            j++;
        }
    }
    return i;
}

int main(void)
{
    srand ( SEED );
    int* bin;
    int binning;
    int k=0;
    int l=0;
        int iter;
    fprintf(stdout, "\nIterations: ");
    fscanf(stdin, "%d", &iter);

    bin = malloc(21*sizeof(int));
    
    while (k < 20)
    {   
        bin[k] = 0;
        k++;
    }

    k = 0;

    while(k < iter) // Count of particle ot distribute
    {
        binning = distr(); // Choosing the bin
        bin[binning]+=1; // Counting binned particle per bin
        k++;
    }

    while(l < 20)
    {
        fprintf(stdout, "\n %d", bin[l]);
        l++;
    }
    return EXIT_SUCCESS;
}

I can't wait to read you and thanks in advance, Maleware

Answer 1

Among several other issues, there are a couple of logical errors in the implementation of the distr function and its caller.

In the linked video, the "obstacles" are shaped in lines forming a triangle:

          .              Falling balls.
           .        
         . *
       . * . *
     . *   * . *          Obstacles.
   . *   * . *   *
 . *   * . *   *   * 
 .       .
 . |   | o |   |   |      Allowed result.  
 . |___|___|___|___|        
 .
 o                        Rejected.

Note that the obstacles are spread in staggered lines, every bin can be "fed" by two obstacle, and that some balls will end up outside the allowed bins.

The posted code seems to implement another model:

                     .
               -1    .    +1
                  <--.-->        
                   .[*].                If the choice is between +1 and -1...
                .         .       
               .           .
               .           .
             .[*].  [X]   [*].          One out of two obstacles can't be reached.
          .         .           .                           
         .           .           .               
         .           .           .
        [*].  [x]  .[*]   [X]   [*].           
              . .                     .             
               .                       . 
               .                       .
               .                       .
|     |  X  |  o  |  X  |     |  X  |  .  |     The same happens to the bins.        
|     |     |  .  |     |     |     |  o  |   
|     |     |  o  |     |     |     |     |   
|_____|_____|_____|_____|_____|_____|_____|        

There are no rejections, only unreachable bins, the odd ones or the even ones depending on the number of lines of obstacles.

If there are enough lines of obstacles, more than the number of bins, the "balls" can't escape outside and they fall into the corresponding obstacle in the successive line (effectively a 0 movement). Now they start spreading towards the center, filling all the bins, in what is definetely not a normal distribution.

I'd suggest to consider a different random walk, one that advance or not by one:

    |
[1] .\.               
[0] . | . 
[0] . |   .
[1] .  \    .
[0] .   |     .
[1] .    \      .
[1] .      \      .
[0] .       |       .
            v
   | | | | | | | | | |
    0 1 2 3 4 5 6 7 8 
            ^             Result

Which could be implemented as easily as this unoptimized function

int random_walk(int steps)   
{
    int count = 0;
    while ( steps-- )
    {
        count += rand() > RAND_MAX / 2;
    }
    return count;
}

Note, though, that it uses only one bit of the value returned by rand and that the final result is the total number of non-zero bits. I leave all the possible optimization to the reader.

Answer 2

Apparently with the same even number of +1 or -1's you will always get an even number for the final bin. so if you start as you did with an odd bin starting point the final bin will be odd. Start with an even bin # as I did you will get an even final bin. So what I did was randomly start with j = 1 or 0 so you get about half odd number of movements and half even. I reduced the number of iterations to 50 and increased the number of bins so that most of the results (99%) are captured. You now get a nice normal distribution.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#define SEED time(NULL)


int rando(void) // Monte Carlo method for choosing weather left or right
{
    return rand() & 1;
}  


int distr(void) // Binning particle
{
    int i=20; // Center of bin  
    int j=rand() & 1; // makes the number if movements either odd or even
    int k=0;

    while(j < 50) // Number of 50 or 49 changes of direction 
    {
        k=rando();
        if(k == 1)
        {
            i++;
            j++;
            printf(" i = %d ", i);
        }
    
        if(k == 0)
        {
            i--;
            j++;
            printf(" i = %d ", i);
        }
    }
    return i;
}

int main(void)
{
    srand ( SEED );
    int* bin;
    int binning;
    int k=0;
    int iter;
    printf("\nIterations: ");
    scanf("%d", &iter);

    bin = malloc(21*sizeof(int));

    while (k < 40)  // zero's out bin[0] to bin[39]?
    {   
        bin[k] = 0;
        k++;
    }

    k = 0;

    while(k < iter) // Count of particle of distribute
    {
        binning = distr(); // Choosing the bin
        printf("binning = %d ", binning);
        bin[binning]+=1; // Counting binned particle per bin
        k++;
    }
    int l = 0;
    while(l < 40)
    {
        printf("\n %d", bin[l]);
        l++;
    } 
    /* total the mumber of iterations of distr() */
    int total = 0;
    l = 0;
    while(l < 40)
    {
        total += bin[l];
        l++;
    } 
    printf("\n total number is %d\n\n", total);

    return 0; 

}