c++ - 从 R 翻译的随机生成代码在 C++ 中失败

Question

我正在编写实现随机生成算法的代码，以从正态分布的尾部进行采样 proposed by Christian Robert .问题是，虽然 R 中的代码可以正常工作，但如果失败，则在将其转换为 C++ 之后。我看不出有任何原因，如果能向我解释出了什么问题以及原因，我将不胜感激。

请注意，下面的代码远谈不上优雅和高效，它经过简化以制作可重现的示例。

这是 R 中的函数:

rtnormR <- function(mean = 0, sd = 1, lower = -Inf, upper = Inf) {
  lower <- (lower - mean) / sd 
  upper <- (upper - mean) / sd

  if (lower < upper && lower >= 0) {
    while (TRUE) {
      astar <- (lower + sqrt(lower^2 + 4)) / 2
      z <- rexp(1, astar) + lower
      u <- runif(1)
      if ((u <= exp(-(z - astar)^2 / 2)) && (z <= upper)) break
    }
  } else {
    z <- NaN
  }
  z*sd + mean
}

这里是 C++ 版本:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

double rtnormCpp(double mean, double sd, double lower, double upper) {
  double z_lower = (lower - mean) / sd;
  double z_upper = (upper - mean) / sd;
  bool stop = false;
  double astar, z, u;

  if (z_lower < z_upper && z_lower >= 0) {
    while (!stop) {
      astar = (z_lower + std::sqrt(std::pow(z_lower, 2) + 4)) / 2;
      z = R::exp_rand() * astar + z_lower;
      u = R::unif_rand();
      if ((u <= std::exp(-std::pow(z-astar, 2) / 2)) && (z <= z_upper))
        stop = true;
    }
  } else {
    z = NAN;
  }
  return z*sd + mean;
}

现在比较使用这两个函数获得的样本(它们与来自 msm 库的 dtnorm 函数进行比较):

xx = seq(-6, 6, by = 0.001)
hist(replicate(5000, rtnormR(mean = 0, sd = 1, lower = 3, upper = 5)), freq= FALSE, ylab = "", xlab = "", main = "rtnormR")
lines(xx, msm::dtnorm(xx, mean = 0, sd = 1, lower = 3, upper = 5), col = "red")
hist(replicate(5000, rtnormCpp(mean = 0, sd = 1, lower = 3, upper = 5)), freq= FALSE, ylab = "", xlab = "", main = "rtnormCpp")
lines(xx, msm::dtnorm(xx, mean = 0, sd = 1, lower = 3, upper = 5), col = "red")

如您所见，rtnormCpp 返回有偏差的样本。你知道为什么吗？

Answer 1

虽然可以在 rexp() 中使用 scale 或 rate，但默认参数化是 rate -所以 rexp(1,astar) 的平均值是 1/astar，而不是 astar。

如果将相关的C++代码行改为

z = R::exp_rand() / astar + z_lower;

似乎一切正常。