Случайный код генерации, переведенный из R, не выполняется в Stack Overflow

Я работаю над кодом, реализующим алгоритм случайной генерации для выборки из хвостов нормального распределения предложенный Кристианом Робертом. Проблема в том, что хотя код в 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;
}

Теперь сравните образцы, полученные с использованием обеих функций (они сравниваются с dtnorm функция от MSM библиотека):

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 возвращает смещенные образцы. У тебя есть идеи почему?