c++ - 如何最有效地防止正态分布的随机变量为零？

Question

我正在编写一个蒙特卡洛算法，其中某一时刻我需要除以一个随机变量。更准确地说:随机变量用作差商的步长，因此我实际上首先将某个值乘以该变量，然后再次将其除以该表达式的某个局部线性函数。喜欢

double f(double);

std::tr1::variate_generator<std::tr1::mt19937, std::tr1::normal_distribution<> >
  r( std::tr1::mt19937(time(NULL)),
     std::tr1::normal_distribution<>(0) );

double h = r();
double a = ( f(x+h) - f(x) ) / h;

这在大多数情况下工作正常，但在 h=0 时失败。从数学上讲，这不是一个问题，因为在正态分布随机变量的任何有限(或者实际上是可数)选择中，所有这些变量都将非零且概率为 1。但是在数字实现中我会遇到 h= =0 每 ≈2³² 函数调用一次(不管梅森扭曲器的周期比宇宙还长，它仍然输出普通的 longs!)。

避免这个麻烦很简单，目前我正在做

double h = r();
while (h==0) h=r();

但我不认为这特别优雅。有没有更好的办法？

---

我正在评估的函数实际上不仅仅是像 f 那样简单的 ℝ->ℝ，而是一个 ℝᵐxℝⁿ -> ℝ，我在其中计算 ℝᵐ 变量的梯度，同时对 ℝⁿ 变量进行数值积分.整个函数叠加了不可预测(但“相干”)的噪声，有时具有特定(但未知)的突出频率，这就是我在尝试使用 h 的固定值时遇到麻烦的原因。

Answer 1

你的方式看起来很优雅，也许有点不同:

do {
    h = r();
} while (h == 0.0);