Pnorm and Dnorm on log scale

[adamlauretig]

Thanks for the package, I've found it's really useful!

I ran into overflow/underflow problems, and was able to re-implement them on the log scale, which fixed this problem. I thought I'd share it, in case this is a problem others have run into:

I did not do this using the truncated above and below, only for the separate ones, since those were the only two I'm using. Here's the c++ code:

`
using namespace Rcpp ;

/// No Truncation
double e0 (const double mean,
const double sd,
const double low,
const double high
) {
return(mean) ;
}

/// Truncated Below and Above
double e1 (const double mean,
const double sd,
const double low,
const double high
) {
double s_low = (low - mean) / sd ;
double s_high = (high - mean) / sd ;

double q1 = R::dnorm(s_low, 0.0, 1.0, false) ;
double q2 = R::dnorm(s_high, 0.0, 1.0, false) ;
double q3 = R::pnorm(s_low, 0.0, 1.0, true, false) ;
double q4 = R::pnorm(s_high, 0.0, 1.0, true, false) ;

return(mean + sd * ((q1 - q2) / (q4 - q3))) ;

}

/// Truncated Below
double e2 (const double mean,
const double sd,
const double low
) {
double s_low = (low - mean) / sd ;

double q1 = R::dnorm(s_low, 0.0, 1.0, true) ;
double q3 = R::pnorm(s_low, 0.0, 1.0, true, true) ;

return(mean + sd * exp((q1) - (1 - q3))) ;

}

/// Truncated Above
double e3 (const double mean,
const double sd,
const double high
) {
double s_high = (high - mean) / sd ;

double q2 = R::dnorm(s_high, 0.0, 1.0, false) ;
double q4 = R::pnorm(s_high, 0.0, 1.0, true, false) ;

return( mean - sd * exp(q2 - q4)) ;

}

/// Main Function
double etn2(const double mean,
const double sd,
const double low,
const double high
) {
// Init Useful Values
double out = NA_REAL ;
//

if (low == R_NegInf &&
    high == R_PosInf
    ) {
  // No Truncation
    out = e0(mean, sd, low, high) ;
} else if (low == R_NegInf) {
    // Truncated Above
    out = e3(mean, sd, high) ;
} else if (high == R_PosInf) {
    // Truncation Below
    out = e2(mean, sd, low) ;
} else {
    // Truncation Above and Below
    out = e1(mean, sd, low, high) ;
}
// Check if Mirror Problem is Numerically Better
if (!std::isfinite(out)) {
    return(-e1(-mean, sd, -high, -low)) ;
}
//
return(out) ;

}`