From b26cd6a4fd57e7e904fcce59654da9be79f9e14e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Mon, 19 Aug 2019 14:11:04 +0200 Subject: [PATCH] [apps/proba] Fix Chi2 inverse DCF interval computation --- apps/probability/law/chi_squared_law.cpp | 20 +++++++++++++++++--- apps/probability/law/chi_squared_law.h | 2 +- apps/probability/law/law.cpp | 2 +- 3 files changed, 19 insertions(+), 5 deletions(-) diff --git a/apps/probability/law/chi_squared_law.cpp b/apps/probability/law/chi_squared_law.cpp index 75736377e..d2bcbc641 100644 --- a/apps/probability/law/chi_squared_law.cpp +++ b/apps/probability/law/chi_squared_law.cpp @@ -4,6 +4,8 @@ namespace Probability { +static inline double maxDouble(double x, double y) { return x > y ? x : y; } + float ChiSquaredLaw::xMin() const { return -k_displayLeftMarginRatio * xMax(); } @@ -50,9 +52,21 @@ double ChiSquaredLaw::cumulativeDistributiveFunctionAtAbscissa(double x) const { } double ChiSquaredLaw::cumulativeDistributiveInverseForProbability(double * probability) { - /* We cannot put xMin because xMin is < 0 for fisplay purposes, and negative - * values are not accepted. */ - return cumulativeDistributiveInverseForProbabilityUsingBrentRoots(probability, DBL_EPSILON, xMax()); + /* We have to compute the values of the interval in chich to look for x. + * We cannot put xMin because xMin is < 0 for display purposes, and negative + * values are not accepted. + * The maximum of the interval: we want + * cumulativeDistributiveFunctionAtAbscissa(b) > proba + * => regularizedGamma(halfk, b/2, k_regularizedGammaPrecision) >proba + * => int(exp(-t)*t^(k/2 - 1), t, 0, b/2)/gamma(k/2) > proba + * => int(exp(-t)*t^(k/2 - 1), t, 0, b/2) > proba * gamma(k/2) + * true if: (b/2 - 0) * exp(-(k/2 - 1))*(k/2 - 1)^(k/2 - 1) > proba * gamma(k/2) + * => b/2 * exp(-(k/2 - 1) (1 + ln(k/2 - 1))) > proba * gamma(k/2) + * => b > 2 * proba * gamma(k/2) / exp(-(k/2 - 1) (1 + ln(k/2 - 1))) */ + const double ceilKOver2 = std::ceil(m_parameter1/2); + const double kOver2Minus1 = std::ceil(m_parameter1/2); + double xmax = 2 * *probability * std::exp(std::lgamma(ceilKOver2)) / exp(-(kOver2Minus1)*(1 + std::log(kOver2Minus1))); + return cumulativeDistributiveInverseForProbabilityUsingBrentRoots(probability, DBL_EPSILON, maxDouble(xMax(), xmax)); } float ChiSquaredLaw::coefficient() const { diff --git a/apps/probability/law/chi_squared_law.h b/apps/probability/law/chi_squared_law.h index 5ee3c2d3c..0edd0c666 100644 --- a/apps/probability/law/chi_squared_law.h +++ b/apps/probability/law/chi_squared_law.h @@ -9,7 +9,7 @@ namespace Probability { class ChiSquaredLaw : public OneParameterLaw { public: static constexpr int k_maxRegularizedGammaIterations = 1000; - static constexpr double k_regularizedGammaPrecision = FLT_EPSILON; + static constexpr double k_regularizedGammaPrecision = DBL_EPSILON; ChiSquaredLaw() : OneParameterLaw(1.0f) {} I18n::Message title() override { return I18n::Message::ChiSquaredLaw; } diff --git a/apps/probability/law/law.cpp b/apps/probability/law/law.cpp index 073fb59ec..71440aba5 100644 --- a/apps/probability/law/law.cpp +++ b/apps/probability/law/law.cpp @@ -149,7 +149,7 @@ double Law::cumulativeDistributiveInverseForProbabilityUsingBrentRoots(double * this, probability, nullptr); - assert(std::fabs(result.value()) < FLT_EPSILON*10); // TODO FLT_EPSILON is too strict + assert(std::fabs(result.value()) < FLT_EPSILON*100); // TODO FLT_EPSILON is too strict return result.abscissa(); }