使用 rand() 的高斯分布

Question

我准备在 C 中构建一个基于 Monte Carlo 的模拟，因为它非常快，我想知道为什么我的代码会产生均匀分布。

所以首先要接你，在我开始编码之前，我想象了几个球掉下来的图片，在每个点上它们可以随机分布向左或向右移动。图片显示在这里：https://www.youtube.com/watch?v=PM7z_03o_kk。

但我得到的有点奇怪。当我将散点设置为 10（在代码示例中设置为 100）时：

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

我得到了一张像高斯分布的图片，但只有 bin 对它有贡献。当它像代码中显示的那样足够大时，每个 bin 都会得到大约相同数量的粒子。这是 2D 案例，一旦按预期工作，它将扩展到 3D 案例。

其中仍然有一些变量不是真正必要的，只是为了避免我能想象到的任何可能的错误。我正在使用 gdb 来查找错误。当我刚刚用 gdb 运行 distr() 时，我发现如果上面的示例设置为 10，它只会产生偶数。当我转到 11 时，我发现 bin[0] 开始贡献很小的数量与其他人相比。我还跑了足够多的时间来看看它是否真的是伪随机的，我发现它应该可以工作。

我仍然无法弄清楚错误，所以我真的希望这里有足够多的人比我更聪明oO。

完整代码：

#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;
}

我迫不及待地想读你的书，提前谢谢你， 恶意软件

Answer 1

在其他几个问题中，distr 函数及其调用者的实现中存在一些逻辑错误。

在链接的视频中，“障碍物”呈直线形，形成一个三角形：

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

请注意，障碍物以交错线分布，每个箱子都可以被两个障碍物“喂食”，并且一些球最终会在外面 允许的垃圾箱。

贴出的代码似乎实现了另一个模型：

                     .
               -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  |     |     |     |     |   
|_____|_____|_____|_____|_____|_____|_____|        

没有拒绝，只有无法到达的垃圾箱，奇数或偶数取决于障碍物的行数。

如果有足够多的障碍线，超过箱子的数量，“球”就无法逃到外面，它们会在连续的线中落入相应的障碍物（实际上是0运动）。现在它们开始向中心扩散，填满所有垃圾箱，这绝对是非正态分布。

我建议考虑一个不同的随机游走，一个提前或不提前一个：

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

可以像这个未优化的函数一样轻松实现

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

但请注意，它仅使用rand 返回的值的一位，并且最终结果是非零位的总数。我将所有可能的优化留给读者。

Answer 2

显然，对于相同的偶数 +1 或 -1，您将始终得到最后一个 bin 的偶数。所以如果你从一个奇数的 bin 起点开始，最终的 bin 将是奇数。从偶数 bin # 开始，就像我做的那样，您将得到一个偶数的最终 bin。所以我所做的是随机从 j = 1 或 0 开始，所以你得到大约一半奇数的运动和一半偶数。我将迭代次数减少到 50 次，并增加了 bin 的数量，以便捕获大部分结果 (99%)。你现在得到了一个很好的正态分布。

#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; 

}