Michael's correction

josephmure · Apr 12, 2024 · d7522c2 · d7522c2
1 parent 27ee3c2
commit d7522c2
Showing 1 changed file with 11 additions and 10 deletions.
diff --git a/lib/src/Uncertainty/Distribution/ClaytonCopula.cxx b/lib/src/Uncertainty/Distribution/ClaytonCopula.cxx
@@ -385,18 +385,19 @@ Point ClaytonCopula::computeQuantile(const Scalar prob,
                                      const Bool tail) const
 {
   if (!((prob >= 0.0) && (prob <= 1.0))) throw InvalidArgumentException(HERE) << "Error: cannot compute a quantile for a probability level outside of [0, 1]";
-  const Scalar q = tail ? 1.0 - prob : prob;
-  // Special case for boarding values
-  if (q == 0.0) return getRange().getLowerBound();
-  if (q == 1.0) return getRange().getUpperBound();
+  // Special case for boundary values
+  if (prob == 0.0) return tail ? getRange().getUpperBound() : getRange().getLowerBound();
+  if (prob == 1.0) return tail ? getRange().getLowerBound() : getRange().getUpperBound();
   // Independent case
-  if (theta_ == 0.0) return Point(2, std::sqrt(q));
+  if (theta_ == 0.0) return Point(2, std::sqrt(tail ? 1.0 - prob : prob));
   // General case
-  LOGWARN(OSS() << "In ClaytonCopula::computeQuantile: q - 1 = " << q - 1.0 << " and theta_ = " << theta_);
-  LOGWARN(OSS() << "std::pow(q, -theta_) - 1 = " << std::pow(q, -theta_) - 1.0);
-  LOGWARN(OSS() << "M_LN2 - log1p(std::pow(q, -theta_)) = " << M_LN2 - log1p(std::pow(q, -theta_)));
-  LOGWARN(OSS() << "Return value - 1 = " << std::exp((M_LN2 - log1p(std::pow(q, -theta_))) / theta_) - 1.0);
-  return Point(2, std::exp((M_LN2 - log1p(std::pow(q, -theta_))) / theta_));
+  if (tail)
+  {
+    if (prob < 1e-8) return Point(2, 1.0 - prob / 2.0);
+    return Point(2, std::exp((M_LN2 - log1p(std::exp(-theta_ * log1p(-prob)))) / theta_));
+  }
+  if (1.0 - prob < 1e-8) return Point(2, (1.0 + prob) / 2.0);
+  return Point(2, std::exp((M_LN2 - log1p(std::pow(prob, -theta_))) / theta_));
 }
 
 /* Compute the CDF of Xi | X1, ..., Xi-1. x = Xi, y = (X1,...,Xi-1) */