如何解决使用rand进行高斯分布
我打算用C语言建立一个基于蒙特卡洛的仿真,因为它超级快,我想知道为什么我的代码会产生均匀的分布。
因此,首先要接您,好吧,我开始编码,我拍摄了几个落下的球的图像,在每个点上它们可以以随机分布向左或向右移动。图片显示在这里:https://www.youtube.com/watch?v=PM7z_03o_kk。
但是我得到的有点奇怪。当我将分散点设置为10(在代码示例中设置为100)时:
while(j < 100) // Number of abitrary change of direction
我得到了类似高斯分布的图像,但是只有bin的图像在起作用。当它的大小足够大(如代码所示)时,每个垃圾箱都会得到大约相同数量的粒子。这是2D情况,一旦按预期工作,它将扩展为3D情况。
其中仍然有些变量并不是真正必要的,只是为了避免我可以想象的任何错误。我正在与gdb一起查找错误。当我仅使用gdb运行distr()时,我发现如果将上面的示例设置为10,它只会产生偶数。当我转到11时,我发现bin [0]开始贡献很小的数值与其余的相比。我还运行了rando()足够的时间来查看它是否真的是伪随机的,我发现它应该可以工作。
我仍然无法弄清错误,所以我真的希望这里有足够多的人比我聪明。
完整代码:
#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;
}
我迫不及待想念你,在此先感谢, 恶意软件
解决方法
在其他几个问题中,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
返回的值的一位,并且最终结果是非零位数的总数。我将所有可能的优化留给读者。
显然,如果+1或-1的偶数相同,那么最终的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;
}
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。