我尝试通过Box-Muller生成随机数有什么问题?

Question

我正在尝试使用Box-Muller方法从高斯分布中生成一个随机数。然而，我的结果与事实相去甚远。

你能告诉我我做错了什么吗？

.stat.GetOneGaussianByBoxMuller:{ 
sq:{
a:2.0*rand[abs[system"S"]]%abs[system"S"]-1;
b:2.0*rand[abs[system"S"]]%abs[system"S"]-1;
sq:(a*a)+(b*b);
x:sq
}/[{x>=1};1]; 
:(2.0*rand[abs[system"S"]]%abs[system"S"]-1)*sqrt[(neg[2]*log[sq])%sq]
};

Answer 1

你的答案有三个问题：q运算的优先级，如何在[-1]中生成随机变量;1]，以及变量sq上的条件。

操作优先级：

请注意，下面，a1和a2是不同的：

u:rand[1f];
a1:-1+2.0*u;
a2:2.0*u-1;
u:0.5
a1:0f
a2:-1f
u:0
a1:-1f
a2:-2f

在 [-1;1]：

使用以下内容：

a:-1+2.0*rand[1.0];
b:-1+2.0*rand[1.0];

sq条件 ：

您确实考虑了sq >= 1条件，但是如果sq = 0，您也会遇到问题，因为您必须在最后一步中除以sq。

此外，在您的实现中，您计算了两次a，其中 1) 不是最优的，2) 与方法不一致，因为在最后一步计算sq时必须使用相同的a，后者产生非常大的数字。我从维基百科页面得到了一些灵感，他们建议如果不满足上述条件，则重新生成sq。因此，在下面的实现中，函数中的递归调用：

.stat.GetOneGaussianByBoxMuller:{ 
a:-1+2.0*rand[1.0];
b:-1+2.0*rand[1.0];
sq:(a*a)+(b*b);
if[(sq>=1)|(sq=0);
:.stat.GetOneGaussianByBoxMuller[];
];
:a*sqrt[(neg[2]*log[sq])%sq];
};

现在，您可以查看通过绘制以下数据集的直方图生成的数据：

([]val:{:.stat.GetOneGaussianByBoxMuller[]} each til 100000)

编辑：

实际上，您可以通过为每个函数调用生成 2 个随机数来获得更有效的实现，如下所示：

.stat.GetOneGaussianByBoxMuller:{ 
a:-1+2.0*rand[1.0];
b:-1+2.0*rand[1.0];
sq:(a*a)+(b*b);
if[(sq>=1)|(sq=0);
:.stat.GetOneGaussianByBoxMuller[];
];
:(a*sqrt[(neg[2]*log[sq])%sq];b*sqrt[(neg[2]*log[sq])%sq]);
};
([]val:raze {.stat.GetOneGaussianByBoxMuller[]} each til 50000)

此实现需要150 毫秒才能生成 100000 个正态分布的随机数，而上述实现需要 245 毫秒。