从R转换的随机生成代码在C中失败

我正在研究实现随机生成算法的代码,用于从正常分布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返回有偏差的样本.你有什么想法吗？

解决方法:

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

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

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

一切似乎都很好.

标签：c,r,rcpp
来源： https://codeday.me/bug/20190727/1553784.html