Bug in Phi(x) for x < -8.5 or so

[ksvanhorn]

[This is for Stan 2.0.1] I have found that Phi(x) rounds down to 0 for x < -8.5 or so. Specifically, Phi(-8) does NOT round to zero, but Phi(x) for every x <= -8.52352 I have tried has rounded to zero. For sufficiently negative x this is unavoidable, but the correct value for Phi(x) is representable as a nonzero double all the way down to x = -37.

[bob-carpenter]

Thanks. Given that the implementation is just

 template <typename T>
 inline typename boost::math::tools::promote_args<T>::type
 Phi(const T x) {
   return 0.5 * (1.0 + boost::math::erf(INV_SQRT_2 * x));
 }

this signals a problem in how we're using boost::math::erf

http://www.boost.org/doc/libs/1_54_0/libs/math/doc/html/math_toolkit/sf_erf/error_function.html

I'll reproduce with a test in a branch and report back.

We'll need to test both erf() and erfc() for similar behavior.

• Bob

On 11/8/13, 4:52 PM, Kevin S. Van Horn wrote:

[This is for Stan 2.0.1] I have found that Phi(x) rounds down to 0 for x < -8.5 or so. Specifically, Phi(-8) does NOT
round to zero, but Phi(x) for every x <= -8.52352 I have tried has rounded to zero. For sufficiently negative x this is
unavoidable, but the correct value for Phi(x) is representable as a nonzero double all the way down to x = -37.

—
Reply to this email directly or view it on GitHub #376.

[bob-carpenter]

We've been pointed before at this, which appears to have a simple evaluation strategy for Phi() with a low error rate:

Marsaglia, George. 2004. Evaluating the Normal Distribution. J. Stat. Soft.