etruncnorm: Expectation is not calculated correctly

[jpasquier]

Hi,

It seems that etruncnorm does not calculate the expectation correctly.

Example:

> etruncnorm(0, 1, 2, 0.1)
[1] 0.5
> etruncnorm2 <- function(a, b, m, s) {
+   a <- (a - m) / s
+   b <- (b - m) / s
+   Z <- pnorm(b) - pnorm(a)
+   m + (dnorm(a) - dnorm(b)) / Z * s
+ }
> etruncnorm2(0, 1, 2, 0.1)
[1] 0.9901907

Best regards,

Jérôme

[barbehenna]

The code causing this problem is here https://github.com/olafmersmann/truncnorm/blob/28593c0b454bdbe860bc00e8e23690c12af89927/src/truncnorm.c#L54

  /* Special case numerically instable case when (a, b) is far away from the
   * center of mass. */
  if (b < mean - 6.0 * sd || a > mean + 6.0 * sd)
    return (a + b) / 2.0;

Invoking an approximation when $(b - \mu)/\sigma < -6$ or $(a - \mu)/\sigma > 6$ seems reasonable to me, but the approximation should be something like this:

if (b < mean - 6.0 * sd)
    return b;
if (a > mean + 6.0 * sd)
    return a;