math-functions-0.3.4.2: tests/Tests/SpecFunctions.hs
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE ViewPatterns #-}
-- | Tests for Statistics.Math
module Tests.SpecFunctions (
tests
) where
import Control.Monad
import Data.List
import Data.Maybe
import qualified Data.Vector as V
import Data.Vector ((!))
import qualified Data.Vector.Unboxed as U
import Test.QuickCheck hiding (choose,within)
import Test.Tasty
import Test.Tasty.QuickCheck (testProperty)
import Test.Tasty.HUnit
import Tests.Helpers
import Tests.SpecFunctions.Tables
import Numeric.SpecFunctions
import Numeric.SpecFunctions.Internal (factorialTable)
import Numeric.MathFunctions.Comparison (within,relativeError,ulpDistance)
import Numeric.MathFunctions.Constants (m_epsilon,m_tiny)
erfTol,erfcTol,erfcLargeTol :: Int
-- Pure haskell implementation is not very good
#if !defined(USE_SYSTEM_ERF) || defined(__GHCJS__)
erfTol = 4
erfcTol = 4
erfcLargeTol = 64
-- Macos's erf is slightly worse that GNU one
#elif defined(darwin_HOST_OS)
erfTol = 2
erfcTol = 2
erfcLargeTol = 2
#else
erfTol = 1
erfcTol = 2
erfcLargeTol = 2
#endif
isGHCJS :: Bool
#if defined(__GHCJS__)
isGHCJS = True
#else
isGHCJS = False
#endif
isWindows :: Bool
#if defined(mingw32_HOST_OS)
isWindows = True
#else
isWindows = False
#endif
tests :: TestTree
tests = testGroup "Special functions"
[ testGroup "erf"
[ -- implementation from numerical recipes loses presision for
-- large arguments
testCase "erfc table" $
forTable "tests/tables/erfc.dat" $ \[x, exact] ->
checkTabularPure erfcTol (show x) exact (erfc x)
, testCase "erfc table [large]" $
forTable "tests/tables/erfc-large.dat" $ \[x, exact] ->
checkTabularPure erfcLargeTol (show x) exact (erfc x)
--
, testCase "erf table" $
forTable "tests/tables/erf.dat" $ \[x, exact] -> do
checkTabularPure erfTol (show x) exact (erf x)
, testProperty "id = erfc . invErfc" invErfcIsInverse
, testProperty "id = invErfc . erfc" invErfcIsInverse2
, testProperty "invErf = erf^-1" invErfIsInverse
]
--
, testGroup "log1p & Co"
[ testCase "expm1 table" $
forTable "tests/tables/expm1.dat" $ \[x, exact] ->
checkTabularPure 2 (show x) exact (expm1 x)
, testCase "log1p table" $
forTable "tests/tables/log1p.dat" $ \[x, exact] ->
checkTabularPure 1 (show x) exact (log1p x)
]
----------------
, testGroup "gamma function"
[ testCase "logGamma table [fractional points" $
forTable "tests/tables/loggamma.dat" $ \[x, exact] -> do
checkTabularPure 2 (show x) exact (logGamma x)
, testProperty "Gamma(x+1) = x*Gamma(x)" $ gammaReccurence
, testCase "logGamma is expected to be precise at 1e-15 level" $
forM_ [3..10000::Int] $ \n -> do
let exact = logFactorial (n-1)
val = logGamma (fromIntegral n)
checkTabular 8 (show n) exact val
]
----------------
, testGroup "incomplete gamma"
[ testCase "incompleteGamma table" $
forTable "tests/tables/igamma.dat" $ \[a,x,exact] -> do
let err | a < 10 = 16
| a <= 101 = case () of
_| isGHCJS -> 64
| isWindows -> 64
| otherwise -> 32
| a == 201 = 200
| otherwise = 32
checkTabularPure err (show (a,x)) exact (incompleteGamma a x)
, testProperty "incomplete gamma - increases" $
\(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y
, testProperty "0 <= gamma <= 1" incompleteGammaInRange
, testProperty "gamma(1,x) = 1 - exp(-x)" incompleteGammaAt1Check
, testProperty "invIncompleteGamma = gamma^-1" invIGammaIsInverse
]
----------------
, testGroup "beta function"
[ testCase "logBeta table" $
forTable "tests/tables/logbeta.dat" $ \[p,q,exact] ->
let errEst
-- For Stirling approx. errors are very good
| b > 10 = 2
-- Partial Stirling approx
| a > 10 = case () of
_| b >= 1 -> 4
| otherwise -> 2 * est
-- sum of logGamma
| otherwise = case () of
_| a <= 1 && b <= 1 -> 8
| a >= 1 && b >= 1 -> 8
| otherwise -> 2 * est
where
a = max p q
b = min p q
--
est = ceiling
$ abs (logGamma a) + abs (logGamma b) + abs (logGamma (a + b))
/ abs (logBeta a b)
in checkTabularPure errEst (show (p,q)) exact (logBeta p q)
, testCase "logBeta factorial" betaFactorial
, testProperty "beta(1,p) = 1/p" beta1p
-- , testProperty "beta recurrence" betaRecurrence
]
----------------
, testGroup "incomplete beta"
[ testCase "incompleteBeta table" $
forM_ tableIncompleteBeta $ \(p,q,x,exact) ->
checkTabular 64 (show (x,p,q)) (incompleteBeta p q x) exact
, testCase "incompleteBeta table with p > 3000 and q > 3000" $
forM_ tableIncompleteBetaP3000 $ \(x,p,q,exact) ->
checkTabular 7000 (show (x,p,q)) (incompleteBeta p q x) exact
--
, testProperty "0 <= I[B] <= 1" incompleteBetaInRange
, testProperty "ibeta symmetry" incompleteBetaSymmetry
, testCase "Regression #68" $ do
let a = 1
b = 0.3
p = 0.3
x = invIncompleteBeta a b p
assertBool "Inversion OK" $ incompleteBeta a b x `ulpDistance` p < 4
-- XXX FIXME DISABLED due to failures
-- , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse
]
----------------
, testGroup "digamma"
[ testAssertion "digamma is expected to be precise at 1e-14 [integers]"
$ digammaTestIntegers 1e-14
-- Relative precision is lost when digamma(x) ≈ 0
, testCase "digamma is expected to be precise at 1e-12" $
forTable "tests/tables/digamma.dat" $ \[x, exact] ->
checkTabularPure 2048
(show x) (digamma x) exact
]
----------------
, testGroup "factorial"
[ testCase "Factorial table" $
forM_ [0 .. 170] $ \n -> do
checkTabular 1
(show n)
(fromIntegral (factorial' n))
(factorial (fromIntegral n :: Int))
--
, testCase "Log factorial from integer" $
forM_ [2 .. 170] $ \n -> do
checkTabular 1
(show n)
(log $ fromIntegral $ factorial' n)
(logFactorial (fromIntegral n :: Int))
, testAssertion "Factorial table is OK"
$ U.length factorialTable == 171
, testCase "Log factorial table" $
forTable "tests/tables/factorial.dat" $ \[i,exact] ->
checkTabularPure 3
(show i) (logFactorial (round i :: Int)) exact
]
----------------
, testGroup "combinatorics"
[ testCase "choose table" $
forM_ [0 .. 1000] $ \n ->
forM_ [0 .. n] $ \k -> do
checkTabular (if isWindows then 3072 else 2048)
(show (n,k))
(fromIntegral $ choose' n k)
(choose (fromInteger n) (fromInteger k))
--
, testCase "logChoose == log . choose" $
forM_ [0 .. 1000] $ \n ->
forM_ [0 .. n] $ \k -> do
checkTabular 2
(show (n,k))
(log $ choose n k)
(logChoose n k)
]
----------------------------------------------------------------
-- Self tests
, testGroup "self-test"
[ testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0
]
]
----------------------------------------------------------------
-- efr tests
----------------------------------------------------------------
roundtrip_erfc_invErfc,
roundtrip_invErfc_erfc,
roundtrip_erf_invErf
:: (Double,Double)
#if !defined(USE_SYSTEM_ERF) || defined(__GHCJS__)
roundtrip_erfc_invErfc = (2,8)
roundtrip_invErfc_erfc = (8,4)
roundtrip_erf_invErf = (128,128)
#elif defined(darwin_HOST_OS)
roundtrip_erfc_invErfc = (4,4)
roundtrip_invErfc_erfc = (4,4)
roundtrip_erf_invErf = (2,2)
#else
roundtrip_erfc_invErfc = (2,2)
roundtrip_invErfc_erfc = (2,2)
roundtrip_erf_invErf = (1,1)
#endif
-- id ≈ erfc . invErfc
invErfcIsInverse :: Double -> Property
invErfcIsInverse ((*2) . range01 -> x)
= (not $ isInfinite x) ==>
( counterexample ("x = " ++ show x )
$ counterexample ("y = " ++ show y )
$ counterexample ("x' = " ++ show x')
$ counterexample ("calc.err = " ++ show (delta, delta-e'))
$ counterexample ("ulps = " ++ show (ulpDistance x x'))
$ ulpDistance x x' <= round delta
)
where
(e,e') = roundtrip_erfc_invErfc
delta = e' + e * abs ( y / x * 2 / sqrt pi * exp( -y*y ))
y = invErfc x
x' = erfc y
-- id ≈ invErfc . erfc
invErfcIsInverse2 :: Double -> Property
invErfcIsInverse2 x
= (not $ isInfinite x') ==>
(y > m_tiny) ==>
(x /= 0) ==>
counterexample ("x = " ++ show x )
$ counterexample ("y = " ++ show y )
$ counterexample ("x' = " ++ show x')
$ counterexample ("calc.err = " ++ show delta)
$ counterexample ("ulps = " ++ show (ulpDistance x x'))
$ ulpDistance x x' <= delta
where
(e,e') = roundtrip_invErfc_erfc
delta = round
$ e' + e * abs (y / x / (2 / sqrt pi * exp( -x*x )))
y = erfc x
x' = invErfc y
-- id ≈ erf . invErf
invErfIsInverse :: Double -> Property
invErfIsInverse a
= (x /= 0) ==>
counterexample ("x = " ++ show x )
$ counterexample ("y = " ++ show y )
$ counterexample ("x' = " ++ show x')
$ counterexample ("calc.err = " ++ show delta)
$ counterexample ("ulps = " ++ show (ulpDistance x x'))
$ ulpDistance x x' <= delta
where
(e,e') = roundtrip_erf_invErf
delta = round
$ e + e' * abs (y / x * 2 / sqrt pi * exp ( -y * y ))
x | a < 0 = - range01 a
| otherwise = range01 a
y = invErf x
x' = erf y
----------------------------------------------------------------
-- QC tests
----------------------------------------------------------------
-- B(p,q) = (x - 1)!(y-1)! / (x + y - 1)!
betaFactorial :: IO ()
betaFactorial = do
forM_ prod $ \(p,q,facP,facQ,facProd) -> do
let exact = fromIntegral (facQ * facP)
/ fromIntegral facProd
checkTabular 16 (show (p,q))
(logBeta (fromIntegral p) (fromIntegral q))
(log exact)
where
prod = [ (p,q,facP,facQ, factorial' (p + q - 1))
| (p,facP) <- facList
, (q,facQ) <- facList
, p + q < 170
, not (p == 1 && q== 1)
]
facList = [(p,factorial' (p-1)) | p <- [1 .. 170]]
-- B(1,p) = 1/p
beta1p :: Double -> Property
beta1p (abs -> p)
= p > 2 ==>
counterexample ("p = " ++ show p)
$ counterexample ("logB = " ++ show lb)
$ counterexample ("err = " ++ show d)
$ d <= 24
where
lb = logBeta 1 p
d = ulpDistance lb (- log p)
{-
-- B(p+1,q) = B(p,q) · p/(p+q)
betaRecurrence :: Double -> Double -> Property
betaRecurrence (abs -> p) (abs -> q)
= p > 0 && q > 0 ==>
counterexample ("p = " ++ show p)
$ counterexample ("q = " ++ show q)
$ counterexample ("log B(p,q) = " ++ show (logBeta p q))
$ counterexample ("log B(p+1,q) = " ++ show (logBeta (p+1) q))
$ counterexample ("err = " ++ show d)
$ d <= 128
where
logB = logBeta p q + log (p / (p + q))
logB' = logBeta (p + 1) q
d = ulpDistance logB logB'
-}
-- Γ(x+1) = x·Γ(x)
gammaReccurence :: Double -> Property
gammaReccurence x
= x > 0 ==> err < errEst
where
g1 = logGamma x
g2 = logGamma (x+1)
err = abs (g2 - g1 - log x)
-- logGamma apparently is not as precise for small x. See #59 for details
errEst = max 1e-14
$ 2 * m_epsilon * sum (map abs [ g1 , g2 , log x ])
-- γ(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 -> Bool
incompleteGammaAt1Check (abs -> x) =
ulpDistance (incompleteGamma 1 x) (-expm1(-x)) < 16
-- invIncompleteGamma is inverse of incompleteGamma
invIGammaIsInverse :: Double -> Double -> Property
invIGammaIsInverse (abs -> a) (range01 -> p) =
a > m_tiny && p > m_tiny && p < 1 && x > m_tiny ==>
( counterexample ("a = " ++ show a )
$ counterexample ("p = " ++ show p )
$ counterexample ("x = " ++ show x )
$ counterexample ("p' = " ++ show p')
$ counterexample ("err = " ++ show (ulpDistance p p'))
$ counterexample ("est = " ++ show est)
$ ulpDistance p p' <= est
)
where
x = invIncompleteGamma a p
f' = exp ( log x * (a-1) - x - logGamma a)
p' = incompleteGamma a x
-- FIXME: Test should be rechecked when #42 is fixed
(e,e') = (32,32)
est = round
$ e' + e * abs (x / p * f')
-- I(x;p,q) 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
-- I(0.5; p,p) = 0.5
incompleteBetaSymmetry :: Double -> Property
incompleteBetaSymmetry (abs -> p) =
p > 0 ==>
counterexample ("p = " ++ show p)
$ counterexample ("ib = " ++ show ib)
$ counterexample ("err = " ++ show d)
$ counterexample ("est = " ++ show est)
$ d <= est
where
est | p < 1 = 80
| p < 10 = 200
| otherwise = round $ 6 * p
d = ulpDistance ib 0.5
ib = incompleteBeta p p 0.5
-- invIncompleteBeta is inverse of incompleteBeta
invIBetaIsInverse :: Double -> Double -> Double -> Property
invIBetaIsInverse (abs -> p) (abs -> q) (range01 -> x) =
p > 0 && q > 0 ==> ( counterexample ("p = " ++ show p )
$ counterexample ("q = " ++ show q )
$ counterexample ("x = " ++ show x )
$ counterexample ("x' = " ++ show x')
$ counterexample ("a = " ++ show a)
$ counterexample ("err = " ++ (show $ abs $ (x - x') / x))
$ abs (x - x') <= 1e-12
)
where
x' = incompleteBeta p q a
a = invIncompleteBeta p q x
-- Table for digamma function:
--
-- Uses equality ψ(n) = H_{n-1} - γ where
-- H_{n} = Σ 1/k, k = [1 .. n] - harmonic number
-- γ = 0.57721566490153286060 - Euler-Mascheroni number
digammaTestIntegers :: Double -> Bool
digammaTestIntegers eps
= all (uncurry $ eq eps) $ take 3000 digammaInt
where
ok approx exact = approx
-- Harmonic numbers starting from 0
harmN = scanl (\a n -> a + 1/n) 0 [1::Rational .. ]
gam = 0.57721566490153286060
-- Digamma values
digammaInt = zipWith (\i h -> (digamma i, realToFrac h - gam)) [1..] harmN
----------------------------------------------------------------
-- Unit tests
----------------------------------------------------------------
-- Lookup table for fact factorial calculation. It has fixed size
-- which is bad but it's OK for this particular case
factorial_table :: V.Vector Integer
factorial_table = V.generate 2000 (\n -> product [1..fromIntegral n])
-- Exact implementation of factorial
factorial' :: Integer -> Integer
factorial' n = factorial_table ! fromIntegral n
-- 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 :: (Integer, Double) -> Double) . properFraction
readTable :: FilePath -> IO [[Double]]
readTable
= fmap (fmap (fmap read . words) . lines)
. readFile
forTable :: FilePath -> ([Double] -> Maybe String) -> IO ()
forTable path fun = do
rows <- readTable path
case mapMaybe fun rows of
[] -> return ()
errs -> assertFailure $ intercalate "---\n" errs
checkTabular :: Int -> String -> Double -> Double -> IO ()
checkTabular prec x exact val =
case checkTabularPure prec x exact val of
Nothing -> return ()
Just s -> assertFailure s
checkTabularPure :: Int -> String -> Double -> Double -> Maybe String
checkTabularPure prec x exact val
| within prec exact val = Nothing
| otherwise = Just $ unlines
[ " x = " ++ x
, " expected = " ++ show exact
, " got = " ++ show val
, " ulps diff = " ++ show (ulpDistance exact val)
, " err.est. = " ++ show prec
]