Я работаю над кодом, реализующим алгоритм случайной генерации для выборки из хвостов нормального распределения предложенный Кристианом Робертом. Проблема в том, что хотя код в 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
возвращает смещенные образцы. У тебя есть идеи почему?
Хотя можно использовать либо scale
или же rate
в rexp()
параметризация по умолчанию rate
— так rexp(1,astar)
имеет среднее значение 1/astar
не astar
,
Если вы измените соответствующую строку кода C ++ на
z = R::exp_rand() / astar + z_lower;
Кажется, все работает нормально.
Других решений пока нет …