mirror of
https://github.com/UpsilonNumworks/Upsilon.git
synced 2026-01-19 16:57:31 +01:00
018dd91ca1
Plots are still rendered in float but computations are now in double Change-Id: I7e0a38effb780861b1443ee92a097cd319de3bc8
85 lines
3.0 KiB
C++
85 lines
3.0 KiB
C++
#include "chi_squared_distribution.h"
|
|
#include "regularized_gamma.h"
|
|
#include <cmath>
|
|
#include <algorithm>
|
|
|
|
namespace Probability {
|
|
|
|
float ChiSquaredDistribution::xMin() const {
|
|
return -k_displayLeftMarginRatio * xMax();
|
|
}
|
|
|
|
float ChiSquaredDistribution::xMax() const {
|
|
assert(m_parameter1 != 0.0);
|
|
return (m_parameter1 + 5.0f * std::sqrt(m_parameter1)) * (1.0f + k_displayRightMarginRatio);
|
|
}
|
|
|
|
float ChiSquaredDistribution::yMax() const {
|
|
float result;
|
|
if (m_parameter1/2.0f <= 1.0f + FLT_EPSILON) {
|
|
result = 0.5f;
|
|
} else {
|
|
result = evaluateAtAbscissa(m_parameter1 - 1.0f) * 1.2f;
|
|
}
|
|
return result * (1.0f + k_displayTopMarginRatio);
|
|
}
|
|
|
|
float ChiSquaredDistribution::evaluateAtAbscissa(float x) const {
|
|
if (x < 0.0f) {
|
|
return NAN;
|
|
}
|
|
const float halfk = m_parameter1/2.0;
|
|
const float halfX = x/2.0f;
|
|
return std::exp(-lgamma(halfk) - halfX + (halfk-1.0f) * std::log(halfX)) / 2.0f;
|
|
}
|
|
|
|
bool ChiSquaredDistribution::authorizedValueAtIndex(float x, int index) const {
|
|
assert(index == 0);
|
|
return x > 0.0f && x == (float)((int)x) && x <= k_maxK;
|
|
}
|
|
|
|
double ChiSquaredDistribution::cumulativeDistributiveFunctionAtAbscissa(double x) const {
|
|
if (x < DBL_EPSILON) {
|
|
return 0.0;
|
|
}
|
|
double result = 0.0;
|
|
if (regularizedGamma(m_parameter1/2.0, x/2.0, k_regularizedGammaPrecision, k_maxRegularizedGammaIterations, &result)) {
|
|
return result;
|
|
}
|
|
return NAN;
|
|
}
|
|
|
|
double ChiSquaredDistribution::cumulativeDistributiveInverseForProbability(double * probability) {
|
|
/* We have to compute the values of the interval in which 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: for (k/2 - 1) > 0 which is k > 2
|
|
* (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)))
|
|
* else
|
|
* 2/k * (b/2)^(k/2) > proba * gamma(k/2)
|
|
* => b^(k/2) > k/2 * 2^(k/2) * proba * gamma(k/2)
|
|
* => exp(k/2 * ln(b)) > k/2 * 2^(k/2) * proba * gamma(k/2)
|
|
* => b > exp(2/k * ln(k/2 * 2^(k/2) * proba * gamma(k/2))) */
|
|
if (*probability < DBL_EPSILON) {
|
|
return 0.0;
|
|
}
|
|
const double k = m_parameter1;
|
|
const double ceilKOver2 = std::ceil(k/2.0);
|
|
const double kOver2Minus1 = k/2.0 - 1.0;
|
|
double xmax = m_parameter1 > 2.0 ?
|
|
2.0 * *probability * std::exp(std::lgamma(ceilKOver2)) / (exp(-kOver2Minus1) * std::pow(kOver2Minus1, kOver2Minus1)) :
|
|
30.0; // Ad hoc value
|
|
xmax = std::isnan(xmax) ? 1000000000.0 : xmax;
|
|
return cumulativeDistributiveInverseForProbabilityUsingIncreasingFunctionRoot(probability, FLT_EPSILON, std::max<double>(xMax(), xmax));
|
|
}
|
|
|
|
}
|