bug in my_etruncnorm

[stephens999]

there seems to be a bug here:

ashr:::my_etruncnorm(0,99,-2,3)
[1] 5.795534
truncnorm:::etruncnorm(0,99,-2,3)
[1] 1.795534

[stephens999]

I fixed this problem but now there still seems to be a problem with infinite limits:

ashr:::my_etruncnorm(0,99,-2,3)
[1] 1.795534
ashr:::my_etruncnorm(0,Inf,-2,3)
[1] 0
ashr:::my_etruncnorm(0,9999,-2,3)
[1] 1.795534

[stephens999]

this turned out to be a bug (or feature) of truncnorm::etruncnorm
with infinite limits when a>b

truncnorm:::etruncnorm(-1,-Inf,-2,3)
[1] NA
truncnorm:::etruncnorm(-1,-999,-2,3)
[1] -3.795471
truncnorm:::etruncnorm(-Inf,-1,-2,3)
[1] -3.795471

[willwerscheid]

etruncnorm does not do well when a and/or b are far away from zero, due to the following lines:

  /* 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;

So, for example, truncnorm::etruncnorm(-31.571, -6.379, 0, 1) gives -18.975 (the correct expected value should be -6.529). This often leads to problems for semi-nonnegative flash (and, I'd expect, for any other applications using unimix). The weird thing is that this case shouldn't actually be unstable, given their implementation (all calculations are done on log scale). We could easily implement the function ourselves -- it's not complicated. See here, ll. 47-76:
https://github.com/olafmersmann/truncnorm/blob/master/src/truncnorm.c

[stephens999]

Is this a recent update on etruncnorm?

…

On Fri, Jan 25, 2019 at 11:37 AM Jason Willwerscheid < ***@***.***> wrote: etruncnorm does not do well when a and/or b are far away from zero, due to the following lines: `static R_INLINE double e_truncnorm(double a, double b, double mean, double sd) { /* 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;So, for example,truncnorm::etruncnorm(-31.571, -6.379, 0, 1)gives -18.975 (the correct expected value should be -6.529). This often leads to problems for semi-nonnegative flash (and, I'd expect, for any other applications usingunimix`). The weird thing is that this case shouldn't actually be unstable, given their implementation (all calculations are done on log scale). We could easily implement the function ourselves -- it's not complicated. See here, ll. 47-76: https://github.com/olafmersmann/truncnorm/blob/master/src/truncnorm.c — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#78 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABt4xUCZRsihuo9hQU4-5QlBsXlG5vMjks5vG0DWgaJpZM4Q9spq> .

[willwerscheid]

No, it's been there since 2015.