From ade0b3d0d6794c3990a7395843980d99558e3983 Mon Sep 17 00:00:00 2001 From: Alexey Khudyakov Date: Mon, 4 Jul 2016 19:37:26 +0300 Subject: [PATCH] Avoid underflow during computation of incomplete beta I could lean to severe loss of precision --- Numeric/SpecFunctions/Internal.hs | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs index 4f8da76..f87abff 100644 --- a/Numeric/SpecFunctions/Internal.hs +++ b/Numeric/SpecFunctions/Internal.hs @@ -21,7 +21,7 @@ import Numeric.Polynomial (evaluateEvenPolynomialL,evaluateOddPolyn import Numeric.Series (sumPowerSeries,enumSequenceFrom,scanSequence,evalContFractionB) import Numeric.MathFunctions.Constants ( m_epsilon, m_NaN, m_neg_inf, m_pos_inf , m_sqrt_2_pi, m_ln_sqrt_2_pi, m_eulerMascheroni - , m_min_log + , m_min_log, m_tiny ) import Text.Printf @@ -444,12 +444,19 @@ incompleteBetaWorker beta p q x -- Constants eps = 1e-15 cx = 1 - x - -- Common multiplies for expansion + -- Common multiplies for expansion. Accurate calculation is a bit + -- tricky. Performing calculation in log-domain leads to slight + -- loss of precision for small x, while using ** prone to + -- underflows. -- - -- If beta function becomes soo small it underflows we should - -- perform calculations in log-domain - factor | beta < m_min_log = exp( p * log x + (q - 1) * log cx - beta) - | otherwise = x**p * cx**(q - 1) / exp beta + -- If either beta function of x**p·(1-x)**(q-1) underflows we + -- switch to log domain. It could waste work but there's no easy + -- switch criterion. + factor + | beta < m_min_log || prod < m_tiny = exp( p * log x + (q - 1) * log cx - beta) + | otherwise = prod / exp beta + where + prod = x**p * cx**(q - 1) -- Soper's expansion of incomplete beta function loop !psq (ns :: Int) ai term betain | done = betain' * factor / p