UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-01-19 00:37:25 +01:00
Code Issues Packages Projects Releases Wiki Activity

[apps/proba] Fix Chi2 inverse DCF interval computation

Browse Source
This commit is contained in:
Léa Saviot
2019-08-19 14:11:04 +02:00
parent 07e639aa80
commit b26cd6a4fd
3 changed files with 19 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
apps/probability/law/chi_squared_law.cpp
View File

@@ -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 {

2
apps/probability/law/chi_squared_law.h
View File

@@ -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; }

2
apps/probability/law/law.cpp
View File

@@ -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();
}