From ffae248eb14237e4e67c1a0340dc6a249b974073 Mon Sep 17 00:00:00 2001 From: Alexey Khudyakov Date: Sun, 16 Feb 2025 18:05:01 +0300 Subject: [PATCH] Improve precision for large NDF MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit For large ndf fraction `ndf / (ndf + x*x)` becomes close to 1 and taking logarithm loses precision. It's better to rewrite log(ndf / (ndf + x²)) = -log(1 + x²/ndf) = -log1p(x²/ndf) --- Statistics/Distribution/StudentT.hs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs index c982152f..0be3debc 100644 --- a/Statistics/Distribution/StudentT.hs +++ b/Statistics/Distribution/StudentT.hs @@ -26,7 +26,7 @@ import Data.Binary (Binary(..)) import Data.Data (Data, Typeable) import GHC.Generics (Generic) import Numeric.SpecFunctions ( - logBeta, incompleteBeta, invIncompleteBeta, digamma) + logBeta, incompleteBeta, invIncompleteBeta, digamma, log1p) import qualified Statistics.Distribution as D import Statistics.Distribution.Transform (LinearTransform (..)) @@ -94,8 +94,9 @@ complCumulative (StudentT ndf) x logDensityUnscaled :: StudentT -> Double -> Double -logDensityUnscaled (StudentT ndf) x = - log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf) +logDensityUnscaled (StudentT ndf) x + = log1p (x*x/ndf) * (-(0.5 * (1 + ndf))) + - logBeta 0.5 (0.5 * ndf) quantile :: StudentT -> Double -> Double quantile (StudentT ndf) p