从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,lower = 3,upper = 5)),freq= FALSE,ylab = "",xlab = "",main = "rtnormR")
lines(xx,msm::dtnorm(xx,mean = 0,upper = 5),col = "red")
hist(replicate(5000,rtnormCpp(mean = 0,main = "rtnormCpp")
lines(xx,col = "red")

如您所见,rtnormCpp返回有偏差的样本.你有什么想法吗？

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