diff --git a/include/diagnostics/raftery.hpp b/include/diagnostics/raftery.hpp index de7928202..1fa087a5d 100644 --- a/include/diagnostics/raftery.hpp +++ b/include/diagnostics/raftery.hpp @@ -29,10 +29,9 @@ #define DIAGNOSTICS_RAFTERY_HPP template -NT fix(NT x) +NT round_to_zero(NT x) { - if (x > 0.0) return std::floor(x); - return std::ceil(x); + return (x > 0.0) ? std::floor(x) : std::ceil(x); } #include "raftery_subroutines/empquant.hpp" @@ -121,13 +120,13 @@ MT perform_raftery(MT const& samples, NT const& q, NT const& r, NT const& s) NT psum = alpha + beta; NT tmp1 = std::log(psum * epss / std::max(alpha, beta)) / std::log(std::abs(1.0 - psum)); - NT nburn = fix((tmp1 + 1.0) * NT(kthin)); + NT nburn = round_to_zero((tmp1 + 1.0) * NT(kthin)); NT phi = ppnd((s + 1.0) / 2.0); NT tmp2 = (2.0 - psum) * alpha * beta * (phi * phi) / (psum * psum * psum * (r * r)); - NT nprec = fix(tmp2 + 1.0) * kthin; - NT nmin = fix(((1.0 - q) * q * (phi * phi) / (r * r)) + 1.0); + NT nprec = round_to_zero(tmp2 + 1.0) * kthin; + NT nmin = round_to_zero(((1.0 - q) * q * (phi * phi) / (r * r)) + 1.0); NT irl = (nburn + nprec) / nmin; - NT kind = std::max(fix(irl + 1.0), NT(kmind)); + NT kind = std::max(round_to_zero(irl + 1.0), NT(kmind)); results(i, 0) = NT(kthin); results(i, 1) = NT(nburn); diff --git a/include/diagnostics/raftery_subroutines/empquant.hpp b/include/diagnostics/raftery_subroutines/empquant.hpp index cd305372e..6604f0ad7 100644 --- a/include/diagnostics/raftery_subroutines/empquant.hpp +++ b/include/diagnostics/raftery_subroutines/empquant.hpp @@ -17,10 +17,10 @@ template NT empquant(VT const& sorted_samples, NT const& q) { unsigned int n = sorted_samples.rows(); - + NT order = (n - 1) * q + 1.0; NT fract = order - NT(int(order)); - int low = std::max(fix(order), 1.0); + int low = std::max(round_to_zero(order), 1.0); int high = std::min(low + 1.0, NT(n)); NT y = (1.0 - fract) * sorted_samples(low - 1) + fract * sorted_samples(high - 1);