diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -1,7 +1,7 @@ -{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} -- | -- Module : Numeric.SpecFunctions --- Copyright : (c) 2009, 2011 Bryan O'Sullivan +-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan -- License : BSD3 -- -- Maintainer : bos@serpentine.com @@ -71,7 +71,7 @@ ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5) | x < 12 = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) / ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5) - | x > 5.1e5 = k + | x > 3e6 = k | otherwise = k + x1 * ((r4_2 * x2 + r4_1) * x2 + r4_0) / ((x2 + r4_4) * x2 + r4_3) @@ -168,17 +168,18 @@ -> Double -- ^ /x/ -> Double incompleteGamma p x - | x < 0 || p <= 0 = m_pos_inf - | x == 0 = 0 - | p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1)) - | x >= 1e8 = 1 - | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1) - g = a + log (pearson p 1 1) - in if g > limit then exp g else 0 - | otherwise = let g = p * log x - x - logGamma p + log cf - in if g > limit then 1 - exp g else 1 + | isNaN p || isNaN x = m_NaN + | x < 0 || p <= 0 = m_pos_inf + | x == 0 = 0 + | p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1)) + | x >= 1e8 = 1 + | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1) + g = a + log (pearson p 1 1) + in if g > limit then exp g else 0 + | otherwise = let g = p * log x - x - logGamma p + log cf + in if g > limit then 1 - exp g else 1 where - norm a = erfc (- a / m_sqrt_2) + norm a = 0.5 * erfc (- a / m_sqrt_2) pearson !a !c !g | c' <= tolerance = g' | otherwise = pearson a' c' g' @@ -231,24 +232,25 @@ -- Solve equation γ(a,x) = p using Halley method loop :: Int -> Double -> Double loop i x - | i >= 12 = x - | otherwise = - let - -- Value of γ(a,x) - p - f = incompleteGamma a x - p - -- dγ(a,x)/dx - f' | a > 1 = afac * exp( -(x - a1) + a1 * (log x - lna1)) - | otherwise = exp( -x + a1 * log x - gln) - u = f / f' - -- Halley correction to Newton-Rapson step - corr = u * (a1 / x - 1) - dx = u / (1 - 0.5 * min 1.0 corr) - -- New approximation to x - x' | x < dx = 0.5 * x -- Do not go below 0 - | otherwise = x - dx - in if abs dx < eps * x' - then x' - else loop (i+1) x' + | i >= 12 = x' + -- For small s derivative becomes approximately 1/x*exp(-x) and + -- skyrockets for small x. If it happens correct answer is 0. + | isInfinite f' = 0 + | abs dx < eps * x' = x' + | otherwise = loop (i + 1) x' + where + -- Value of γ(a,x) - p + f = incompleteGamma a x - p + -- dγ(a,x)/dx + f' | a > 1 = afac * exp( -(x - a1) + a1 * (log x - lna1)) + | otherwise = exp( -x + a1 * log x - gln) + u = f / f' + -- Halley correction to Newton-Rapson step + corr = u * (a1 / x - 1) + dx = u / (1 - 0.5 * min 1.0 corr) + -- New approximation to x + x' | x < dx = 0.5 * x -- Do not go below 0 + | otherwise = x - dx -- Calculate inital guess for root guess -- @@ -297,7 +299,7 @@ c = logGammaCorrection q - logGammaCorrection pq -- | Regularized incomplete beta function. Uses algorithm AS63 by --- Majumder abd Bhattachrjee. +-- Majumder and Bhattachrjee. incompleteBeta :: Double -- ^ /p/ > 0 -> Double -- ^ /q/ > 0 -> Double -- ^ /x/, must lie in [0,1] range @@ -312,22 +314,22 @@ -> Double -- ^ /x/, must lie in [0,1] range -> Double incompleteBeta_ beta p q x - | p <= 0 || q <= 0 = error "p <= 0 || q <= 0" - | x < 0 || x > 1 = error "x < 0 || x > 1" - | x == 0 || x == 1 = x + | p <= 0 || q <= 0 = error "p <= 0 || q <= 0" + | x < 0 || x > 1 || isNaN x = error "x out of [0,1] range" + | x == 0 || x == 1 = x | p >= (p+q) * x = incompleteBetaWorker beta p q x | otherwise = 1 - incompleteBetaWorker beta q p (1 - x) -- Worker for incomplete beta function. It is separate function to -- avoid confusion with parameter during parameter swapping incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double -incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q) :: Int) 1 1 1 +incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1 where -- Constants eps = 1e-15 cx = 1 - x -- Loop - loop psq ns ai term betain + loop !psq (ns :: Int) ai term betain | done = betain' * exp( p * log x + (q - 1) * log cx - beta) / p | otherwise = loop psq' (ns - 1) (ai + 1) term' betain' where diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -1,5 +1,5 @@ name: math-functions -version: 0.1.1.1 +version: 0.1.1.2 cabal-version: >= 1.8 license: BSD3 license-file: LICENSE @@ -31,10 +31,12 @@ Numeric.MathFunctions.Constants test-suite tests + buildable: False type: exitcode-stdio-1.0 hs-source-dirs: tests main-is: tests.hs other-modules: + Tests.Helpers Tests.Chebyshev Tests.SpecFunctions Tests.SpecFunctions.Tables diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -21,9 +21,11 @@ [ testProperty "Γ(x+1) = x·Γ(x) logGamma" $ gammaReccurence logGamma 3e-8 , testProperty "Γ(x+1) = x·Γ(x) logGammaL" $ gammaReccurence logGammaL 2e-13 , testProperty "γ(1,x) = 1 - exp(-x)" $ incompleteGammaAt1Check + , testProperty "0 <= γ <= 1" $ incompleteGammaInRange , testProperty "γ - increases" $ \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y , testProperty "invIncompleteGamma = γ^-1" $ invIGammaIsInverse + , testProperty "0 <= I[B] <= 1" $ incompleteBetaInRange , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse -- Unit tests , testAssertion "Factorial is expected to be precise at 1e-15 level" @@ -57,6 +59,9 @@ , testAssertion "choose is expected to precise at 1e-12 level" $ and [ eq 1e-12 (choose (fromIntegral n) (fromIntegral k)) (fromIntegral $ choose' n k) | n <- [0..300], k <- [0..n]] + ---------------------------------------------------------------- + -- Self tests + , testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0 ] ---------------------------------------------------------------- @@ -71,19 +76,25 @@ g1 = logG x g2 = logG (x+1) +-- γ(s,x) is in [0,1] range +incompleteGammaInRange :: Double -> Double -> Property +incompleteGammaInRange (abs -> s) (abs -> x) = + x >= 0 && s > 0 ==> let i = incompleteGamma s x in i >= 0 && i <= 1 -- γ(1,x) = 1 - exp(-x) -- Since Γ(1) = 1 normalization doesn't make any difference incompleteGammaAt1Check :: Double -> Property -incompleteGammaAt1Check x = +incompleteGammaAt1Check (abs -> x) = x > 0 ==> (incompleteGamma 1 x + exp(-x)) ≈ 1 where (≈) = eq 1e-13 -- invIncompleteGamma is inverse of incompleteGamma invIGammaIsInverse :: Double -> Double -> Property -invIGammaIsInverse (abs -> a) (abs . snd . properFraction -> p) = - a > 0 && p > 0 && p < 1 ==> ( printTestCase ("x = " ++ show x ) +invIGammaIsInverse (abs -> a) (range01 -> p) = + a > 0 && p > 0 && p < 1 ==> ( printTestCase ("a = " ++ show a ) + $ printTestCase ("p = " ++ show p ) + $ printTestCase ("x = " ++ show x ) $ printTestCase ("p' = " ++ show p') $ printTestCase ("Δp = " ++ show (p - p')) $ abs (p - p') <= 1e-12 @@ -92,6 +103,11 @@ x = invIncompleteGamma a p p' = incompleteGamma a x +-- B(s,x) is in [0,1] range +incompleteBetaInRange :: Double -> Double -> Double -> Property +incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) = + p > 0 && q > 0 ==> let i = incompleteBeta p q x in i >= 0 && i <= 1 + -- invIncompleteBeta is inverse of incompleteBeta invIBetaIsInverse :: Double -> Double -> Double -> Property invIBetaIsInverse (abs -> p) (abs -> q) (abs . snd . properFraction -> x) = @@ -125,3 +141,7 @@ -- Exact albeit slow implementation of choose choose' :: Integer -> Integer -> Integer choose' n k = factorial' n `div` (factorial' k * factorial' (n-k)) + +-- Truncate double to [0,1] +range01 :: Double -> Double +range01 = abs . snd . properFraction diff --git a/tests/Tests/SpecFunctions/Tables.hs b/tests/Tests/SpecFunctions/Tables.hs new file mode 100644 --- /dev/null +++ b/tests/Tests/SpecFunctions/Tables.hs @@ -0,0 +1,47 @@ +module Tests.SpecFunctions.Tables where + +tableLogGamma :: [(Double,Double)] +tableLogGamma = + [(0.000001250000000, 13.592366285131769033) + , (0.000068200000000, 9.5930266308318756785) + , (0.000246000000000, 8.3100370767447966358) + , (0.000880000000000, 7.03508133735248542) + , (0.003120000000000, 5.768129358365567505) + , (0.026700000000000, 3.6082588918892977148) + , (0.077700000000000, 2.5148371858768232556) + , (0.234000000000000, 1.3579557559432759994) + , (0.860000000000000, 0.098146578027685615897) + , (1.340000000000000, -0.11404757557207759189) + , (1.890000000000000, -0.0425116422978701336) + , (2.450000000000000, 0.25014296569217625565) + , (3.650000000000000, 1.3701041997380685178) + , (4.560000000000000, 2.5375143317949580002) + , (6.660000000000000, 5.9515377269550207018) + , (8.250000000000000, 9.0331869196051233217) + , (11.300000000000001, 15.814180681373947834) + , (25.600000000000001, 56.711261598328121636) + , (50.399999999999999, 146.12815158702164808) + , (123.299999999999997, 468.85500075897556371) + , (487.399999999999977, 2526.9846647543727158) + , (853.399999999999977, 4903.9359135978220365) + , (2923.300000000000182, 20402.93198938705973) + , (8764.299999999999272, 70798.268343590112636) + , (12630.000000000000000, 106641.77264982508495) + , (34500.000000000000000, 325976.34838781820145) + , (82340.000000000000000, 849629.79603036714252) + , (234800.000000000000000, 2668846.4390507959761) + , (834300.000000000000000, 10540830.912557534873) + , (1230000.000000000000000, 16017699.322315014899) + ] +tableIncompleteBeta :: [(Double,Double,Double,Double)] +tableIncompleteBeta = + [(2.000000000000000, 3.000000000000000, 0.030000000000000, 0.0051864299999999996862) + , (2.000000000000000, 3.000000000000000, 0.230000000000000, 0.22845923000000001313) + , (2.000000000000000, 3.000000000000000, 0.760000000000000, 0.95465728000000005249) + , (4.000000000000000, 2.300000000000000, 0.890000000000000, 0.93829812158347802864) + , (1.000000000000000, 1.000000000000000, 0.550000000000000, 0.55000000000000004441) + , (0.300000000000000, 12.199999999999999, 0.110000000000000, 0.95063000053947077639) + , (13.100000000000000, 9.800000000000001, 0.120000000000000, 1.3483109941962659385e-07) + , (13.100000000000000, 9.800000000000001, 0.420000000000000, 0.071321857831804780226) + , (13.100000000000000, 9.800000000000001, 0.920000000000000, 0.99999578339197081611) + ]