Случайный код генерации, переведенный с R, сбой в C++

Я работаю над кодом, реализующим алгоритм случайной генерации для выборки из хвостов нормального распределения proposed by Christian Robert. Проблема в том, что хотя код в R работал правильно, то после перевода его на C++, если он не работает. Я не вижу причин для этого, и я был бы благодарен за то, что объяснил мне, что пошло не так, и почему.Случайный код генерации, переведенный с 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 возвращает необъективные образцы. У вас есть идеи, почему?

источник

2016-10-08 Tim

@BenBolker Right ... Я посмотрел код для 'rexp.c', где' exp_rand() 'умножается на масштаб ... Не стесняйтесь публиковать его как ответ, и я его приму. В противном случае я удалю вопрос. – Tim

Как в стороне, если вы проецируете отдельные точки, вы не видите реальных разниц в скорости. Ваш вариант C++ может быть лучше для векторов этих усеченных нормалей. –

@DirkEddelbuettel Я знаю, это упрощенный пример, не предназначенный для реального использования, сделанный так, чтобы просто быть воспроизводимым. – Tim

Хотя можно использовать либо scale или rate в rexp(), параметризация по умолчанию rate - так rexp(1,astar) имеет среднее 1/astar, не astar.

Если изменить соответствующую строку кода C++ для

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

все, кажется, работает хорошо.

источник

2016-10-08 20:55:13

Случайный код генерации, переведенный с R, сбой в C++

ответ

Смежные вопросы