2

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

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³² 函数调用（不管 mersenne twister 的周期比宇宙长，它仍然输出普通long的 s！）。

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

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

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

我正在评估的函数实际上不仅仅是一个简单的 ℝ->ℝ f，而是一个 ℝᵐxℝⁿ -> ℝ，我在其中计算 ℝᵐ 变量的梯度，同时对 ℝⁿ 变量进行数值积分。整个函数叠加了不可预测的（但“连贯的”）噪声，有时还有特定的（但未知的）突出频率，当我尝试使用h.

4

5 回答 5

3

你的方式看起来足够优雅，也许有点不同：

do {
    h = r();
} while (h == 0.0);
于 2011-06-18T22:19:34.887 回答
3

两个正态分布的随机变量之比就是柯西分布。柯西分布是具有无限方差的那些讨厌的分布之一。确实非常恶心。柯西分布会使您的蒙特卡洛实验一团糟。

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

  • 正态分布通常很容易使用，
  • 通常具有如此好的数学性质，
  • 正常的假设似乎或多或少是正确的，并且
  • 真正的分布是熊。

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

于 2011-06-18T22:35:25.597 回答
0

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

于 2011-06-18T22:13:36.613 回答
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++ 对象，这不太可能成为问题。还值得注意的是，这使得除以零的情况非常昂贵。

于 2011-06-18T22:14:35.523 回答
0

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

如果您实际上正在做的是创建离散点的数量，尽管这近似于正态分布，并且这对于您的分布来说已经足够了，那么创建它的方式是使它们中的任何一个都不会真正具有值 0。

于 2011-06-18T22:29:06.870 回答