subtle bug in InverseCumulativeNormal

[Kakhkhor Abdijalilov]

Bug Report:

Where:            In "normaldistribution.cpp".
What:              Undefined behavior or floating point exception.
Severity:          Subtle.
Burden:           Easy to fix. See below.


Description.
If the argument x is close to 0.0 or 1.0, the implementation sets it
exactly to 0.0 or 1.0.
This leads to an attempt to evaluate std::log(0.0) or std::log(1.0).
I propose that we cut it at x=1E-12 and x = 1 - 1E-12, which covers
all draws between
plus/minus 7 standard deviations. Any argument value beyond that range
should be considered
as either erroneous or astronomically improbable.

Regards,
Kakhkhor Abdijalilov.

=============================================================
// current implementation
Real InverseCumulativeNormal::operator()(Real x) const {
    if (x < 0.0 || x > 1.0) {
        // try to recover if due to numerical error
        if (close_enough(x, 1.0)) {
            x = 1.0;
        } else if (std::fabs(x) < QL_EPSILON) {
            x = 0.0;
        } else {
            QL_FAIL("InverseCumulativeNormal(" << x
                    << ") undefined: must be 0 < x < 1");
        }
    }
.........................................
}

// new implementation
Real InverseCumulativeNormal::operator()(Real x) const {
    if (x < 1e-12 || x > (1.0 - 1e-12)) {
        // try to recover if due to numerical error
        if (close_enough(x, 1.0)) {
            x = 1.0 - 1e-12;
        } else if (std::fabs(x) < QL_EPSILON) {
            x = 1e-12;
        } else {
            QL_FAIL("InverseCumulativeNormal(" << x
                    << ") undefined: must be 0 < x < 1");
        }
    }
.........................................
}

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Nicolai Lassesen-2]

Kakhkhor Abdijalilov <kabdijalilov <at> gmail.com> writes:

>
> Bug Report:
>
> Where:            In "normaldistribution.cpp".
> What:              Undefined behavior or floating point exception.
> Severity:          Subtle.
> Burden:           Easy to fix. See below.
>
> Description.
> If the argument x is close to 0.0 or 1.0, the implementation sets it
> exactly to 0.0 or 1.0.
> This leads to an attempt to evaluate std::log(0.0) or std::log(1.0).
> I propose that we cut it at x=1E-12 and x = 1 - 1E-12, which covers
> all draws between
> plus/minus 7 standard deviations. Any argument value beyond that range
> should be considered
> as either erroneous or astronomically improbable.
>
> Regards,
> Kakhkhor Abdijalilov.

I am currently looking into replacing the Quantlib implementation of
statistical distributions with the implementation provided by Boost. The best
solution (in my opinion anyway) would be to get totally rid of the Quantlib
distribution folder and rely totally on Boost - this would eliminate a layer of
indirection but it would break backward compatibility.

Therefore I have currently just replaced the internal parts of the Quantlib
normaldistribution.xpp with the Boost implementation. Unfortunately these makes
some testcases fail, and I am currently looking into why this can be.

If someone is interested I can post my modified normaldistribution.xpp files
here.

Br,
Nicolai


------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Luigi Ballabio]

On Thu, 2010-05-27 at 01:58 -0500, Kakhkhor Abdijalilov wrote:
> If the argument x is close to 0.0 or 1.0, the implementation sets it
> exactly to 0.0 or 1.0.
> This leads to an attempt to evaluate std::log(0.0) or std::log(1.0).

You're right.  I wonder why it does it though, since 0 and 1 are outside
the domain of the function.  If I were to return a value, I'd return the
maximum allowed double for 1 and the minimum for 0 (to get an
approximation of plus/minus infinity.)  But I'd rather throw an
exception in that case.

Thoughts?

Luigi


--

Call on God, but row away from the rocks.
-- Indian proverb



------------------------------------------------------------------------------
ThinkGeek and WIRED's GeekDad team up for the Ultimate
GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the
lucky parental unit.  See the prize list and enter to win:
http://p.sf.net/sfu/thinkgeek-promo
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev