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

How do I most effectively prevent my normally-distributed random variable from being zero?

本文关键字：随机 变量 正态分布 何最 有效地      更新时间：2023-10-16

我正在写一个蒙特卡罗算法，其中有一点我需要除以一个随机变量。更准确地说：随机变量被用作差商的步长，所以我实际上首先用这个变量乘以一些东西，然后再把它从这个表达式的一些局部线性函数中除以。像

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。但在数字实现中，我每次≈2³²的函数调用都会遇到一个h==0（无论梅森龙卷风的周期比宇宙长，它仍然输出普通的long s！）。

避免这种麻烦很简单，现在我正在做

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

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

我正在评估的函数实际上不仅仅是ℝ->ℝ与f一样ℝᵐxℝⁿ -> ℝ其中我计算了ℝᵐ变量，同时对ℝⁿ变量。整个函数叠加了不可预测（但"相干"）的噪声，有时还有特定（但未知）的突出频率，当我尝试使用固定值的h时，这就是我遇到麻烦的原因。

你的方式似乎足够优雅，也许有点不同：

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

两个正态分布随机变量的比值为柯西分布。柯西分布是一种具有无穷方差的恶劣分布。确实很恶心。柯西分布会把蒙特卡洛实验搞得一团糟。

在计算两个随机变量的比率的许多情况下，分母是不正常的。人们经常使用正态分布来近似这个非正态分布的随机变量，因为

  • 正态分布通常很容易处理
  • 通常具有很好的数学性质
  • 正常假设似乎或多或少是正确的，并且
  • 真正的分配是一只熊

假设你除以距离。根据定义，距离是半正定的，并且通常作为随机变量是正定的。因此，即时距离永远不可能是正态分布的。尽管如此，在平均值远大于标准差的情况下，人们通常假设距离的正态分布。当做出这种正常的假设时，你需要防范那些非真实的价值。一个简单的解决方案是截断法线。

如果要保持正态分布，必须排除0或将0分配给以前未出现的新值。由于在计算机科学的有限范围内，第二种方法很可能是不可能的，所以第一种方法是我们唯一的选择。

函数（f（x+h）-f（x））/h的极限为h->0，因此如果遇到h==0，则应使用该极限。极限是f'（x），所以如果你知道导数，你可以使用它

如果你实际上正在做的是创建一些离散点，尽管这些点近似于正态分布，并且这对你的分布来说已经足够好了，那么以一种没有一个点实际值为0的方式创建它。

根据您试图计算的内容，也许可以使用以下方法：

double h = r();
double a;
if (h != 0)
    a = ( f(x+h) - f(x) ) / h;
else
    a = 0;

如果f是一个线性函数，它应该（我认为？）在h=0时保持连续。

您可能还想考虑通过零异常捕获除法，以避免分支的成本请注意，这可能会也可能不会对性能产生不利影响-双向基准测试！

在Linux上，您需要使用-fnon-call-exceptions构建包含潜在的除以零的文件，并安装SIGFPE处理程序：

struct fp_exception { };
void sigfpe(int) {
  signal(SIGFPE, sigfpe);
  throw fp_exception();
}
void setup() {
  signal(SIGFPE, sigfpe);
}
// Later...
    try {
        run_one_monte_carlo_trial();
    } catch (fp_exception &) {
        // skip this trial
    }

在Windows上，使用SEH:

__try 
{ 
    run_one_monte_carlo_trial();
} 
__except(GetExceptionCode() == EXCEPTION_INT_DIVIDE_BY_ZERO ? 
         EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH)
{ 
    // skip this trial
}

这具有潜在地对快速路径影响较小的优点。没有分支，尽管可能会对异常处理程序记录进行一些调整。在Linux上，由于编译器为-fnon-call-exceptions生成了更保守的代码，因此可能会对性能造成较小的影响。如果在-fnon-call-exceptions下编译的代码不分配任何自动（堆栈）C++对象，那么这就不太可能成为问题。同样值得注意的是，这使得除以零确实的情况非常昂贵。

相关文章：