1{-# LANGUAGE BangPatterns #-}
2{-# LANGUAGE CPP          #-}
3{-# LANGUAGE ViewPatterns #-}
4-- | Tests for Statistics.Math
5module Tests.SpecFunctions (
6  tests
7  ) where
8
9import Control.Monad
10import Data.List
11import Data.Maybe
12import qualified Data.Vector as V
13import           Data.Vector   ((!))
14import qualified Data.Vector.Unboxed as U
15
16import Test.QuickCheck  hiding (choose,within)
17import Test.Tasty
18import Test.Tasty.QuickCheck   (testProperty)
19import Test.Tasty.HUnit
20
21import Tests.Helpers
22import Tests.SpecFunctions.Tables
23import Numeric.SpecFunctions
24import Numeric.SpecFunctions.Internal   (factorialTable)
25import Numeric.MathFunctions.Comparison (within,relativeError,ulpDistance)
26import Numeric.MathFunctions.Constants  (m_epsilon,m_tiny)
27
28erfTol,erfcTol,erfcLargeTol :: Int
29#if USE_SYSTEM_ERF && !defined(__GHCJS__)
30erfTol       = 1
31erfcTol      = 2
32erfcLargeTol = 2
33#else
34erfTol       = 4
35erfcTol      = 4
36erfcLargeTol = 64
37#endif
38
39isGHCJS :: Bool
40#if defined(__GHCJS__)
41isGHCJS = True
42#else
43isGHCJS = False
44#endif
45
46tests :: TestTree
47tests = testGroup "Special functions"
48  [ testGroup "erf"
49    [ -- implementation from numerical recipes loses presision for
50      -- large arguments
51      testCase "erfc table" $
52        forTable "tests/tables/erfc.dat" $ \[x, exact] ->
53          checkTabularPure erfcTol (show x) exact (erfc x)
54    , testCase "erfc table [large]" $
55        forTable "tests/tables/erfc-large.dat" $ \[x, exact] ->
56          checkTabularPure erfcLargeTol (show x) exact (erfc x)
57      --
58    , testCase "erf table" $
59        forTable "tests/tables/erf.dat" $ \[x, exact] -> do
60          checkTabularPure erfTol (show x) exact (erf x)
61    , testProperty "id = erfc . invErfc" invErfcIsInverse
62    , testProperty "id = invErfc . erfc" invErfcIsInverse2
63    , testProperty "invErf  = erf^-1"    invErfIsInverse
64    ]
65  --
66  , testGroup "log1p & Co"
67    [ testCase "expm1 table" $
68        forTable "tests/tables/expm1.dat" $ \[x, exact] ->
69          checkTabularPure 2 (show x) exact (expm1 x)
70    , testCase "log1p table" $
71        forTable "tests/tables/log1p.dat" $ \[x, exact] ->
72          checkTabularPure 1 (show x) exact (log1p x)
73    ]
74  ----------------
75  , testGroup "gamma function"
76    [ testCase "logGamma table [fractional points" $
77        forTable "tests/tables/loggamma.dat" $ \[x, exact] -> do
78          checkTabularPure 2 (show x) exact (logGamma x)
79    , testProperty "Gamma(x+1) = x*Gamma(x)" $ gammaReccurence
80    , testCase     "logGamma is expected to be precise at 1e-15 level" $
81        forM_ [3..10000::Int] $ \n -> do
82          let exact = logFactorial (n-1)
83              val   = logGamma (fromIntegral n)
84          checkTabular 8 (show n) exact val
85    ]
86  ----------------
87  , testGroup "incomplete gamma"
88    [ testCase "incompleteGamma table" $
89        forTable "tests/tables/igamma.dat" $ \[a,x,exact] -> do
90          let err | a < 10    = 16
91                  | a <= 101  = if isGHCJS then 64 else 32
92                  | a == 201  = 200
93                  | otherwise = 32
94          checkTabularPure err (show (a,x)) exact (incompleteGamma a x)
95    , testProperty "incomplete gamma - increases" $
96        \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y
97    , testProperty "0 <= gamma <= 1"               incompleteGammaInRange
98    , testProperty "gamma(1,x) = 1 - exp(-x)"      incompleteGammaAt1Check
99    , testProperty "invIncompleteGamma = gamma^-1" invIGammaIsInverse
100    ]
101  ----------------
102  , testGroup "beta function"
103    [ testCase "logBeta table" $
104        forTable "tests/tables/logbeta.dat" $ \[p,q,exact] ->
105          let errEst
106                -- For Stirling approx. errors are very good
107                | b > 10          = 2
108                -- Partial Stirling approx
109                | a > 10 = case () of
110                    _| b >= 1    -> 4
111                     | otherwise -> 2 * est
112                -- sum of logGamma
113                | otherwise = case () of
114                    _| a <= 1 && b <= 1 -> 8
115                     | a >= 1 && b >= 1 -> 8
116                     | otherwise        -> 2 * est
117                where
118                  a = max p q
119                  b = min p q
120                  --
121                  est = ceiling
122                      $ abs (logGamma a) + abs (logGamma b) + abs (logGamma (a + b))
123                      / abs (logBeta a b)
124          in checkTabularPure errEst (show (p,q)) exact (logBeta p q)
125    , testCase "logBeta factorial" betaFactorial
126    , testProperty "beta(1,p) = 1/p"   beta1p
127    -- , testProperty "beta recurrence"   betaRecurrence
128    ]
129  ----------------
130  , testGroup "incomplete beta"
131    [ testCase "incompleteBeta table" $
132        forM_ tableIncompleteBeta $ \(p,q,x,exact) ->
133          checkTabular 64 (show (x,p,q)) (incompleteBeta p q x) exact
134    , testCase "incompleteBeta table with p > 3000 and q > 3000" $
135        forM_ tableIncompleteBetaP3000 $ \(x,p,q,exact) ->
136          checkTabular 7000 (show (x,p,q)) (incompleteBeta p q x) exact
137    --
138    , testProperty "0 <= I[B] <= 1" incompleteBetaInRange
139    , testProperty "ibeta symmetry" incompleteBetaSymmetry
140    -- XXX FIXME DISABLED due to failures
141    -- , testProperty "invIncompleteBeta  = B^-1" $ invIBetaIsInverse
142    ]
143  ----------------
144  , testGroup "digamma"
145    [ testAssertion "digamma is expected to be precise at 1e-14 [integers]"
146        $ digammaTestIntegers 1e-14
147      -- Relative precision is lost when digamma(x) ≈ 0
148    , testCase "digamma is expected to be precise at 1e-12" $
149      forTable "tests/tables/digamma.dat" $ \[x, exact] ->
150        checkTabularPure 2048
151          (show x) (digamma x) exact
152    ]
153  ----------------
154  , testGroup "factorial"
155    [ testCase "Factorial table" $
156      forM_ [0 .. 170] $ \n -> do
157        checkTabular 1
158          (show n)
159          (fromIntegral (factorial' n))
160          (factorial (fromIntegral n :: Int))
161      --
162    , testCase "Log factorial from integer" $
163      forM_ [2 .. 170] $ \n -> do
164        checkTabular 1
165          (show n)
166          (log $ fromIntegral $ factorial' n)
167          (logFactorial (fromIntegral n :: Int))
168    , testAssertion "Factorial table is OK"
169    $ U.length factorialTable == 171
170    , testCase "Log factorial table" $
171      forTable "tests/tables/factorial.dat" $ \[i,exact] ->
172        checkTabularPure 3
173          (show i) (logFactorial (round i :: Int)) exact
174    ]
175  ----------------
176  , testGroup "combinatorics"
177    [ testCase "choose table" $
178      forM_ [0 .. 1000] $ \n ->
179        forM_ [0 .. n]  $ \k -> do
180          checkTabular 2048
181            (show (n,k))
182            (fromIntegral $ choose' n k)
183            (choose (fromInteger n) (fromInteger k))
184    --
185    , testCase "logChoose == log . choose" $
186      forM_ [0 .. 1000] $ \n ->
187        forM_ [0 .. n]  $ \k -> do
188          checkTabular 2
189            (show (n,k))
190            (log $ choose n k)
191            (logChoose n k)
192    ]
193    ----------------------------------------------------------------
194    -- Self tests
195  , testGroup "self-test"
196    [ testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0
197    ]
198  ]
199
200----------------------------------------------------------------
201-- efr tests
202----------------------------------------------------------------
203
204roundtrip_erfc_invErfc,
205  roundtrip_invErfc_erfc,
206  roundtrip_erf_invErf
207  :: (Double,Double)
208#if USE_SYSTEM_ERF && !defined(__GHCJS__)
209roundtrip_erfc_invErfc = (2,2)
210roundtrip_invErfc_erfc = (2,2)
211roundtrip_erf_invErf   = (1,1)
212#else
213roundtrip_erfc_invErfc = (2,8)
214roundtrip_invErfc_erfc = (8,4)
215roundtrip_erf_invErf   = (128,128)
216#endif
217
218-- id ≈ erfc . invErfc
219invErfcIsInverse :: Double -> Property
220invErfcIsInverse ((*2) . range01 -> x)
221  = (not $ isInfinite x) ==>
222  ( counterexample ("x        = " ++ show x )
223  $ counterexample ("y        = " ++ show y )
224  $ counterexample ("x'       = " ++ show x')
225  $ counterexample ("calc.err = " ++ show (delta, delta-e'))
226  $ counterexample ("ulps     = " ++ show (ulpDistance x x'))
227  $ ulpDistance x x' <= round delta
228  )
229  where
230    (e,e') = roundtrip_erfc_invErfc
231    delta  = e' + e * abs ( y / x  *  2 / sqrt pi * exp( -y*y ))
232    y  = invErfc x
233    x' = erfc y
234
235-- id ≈ invErfc . erfc
236invErfcIsInverse2 :: Double -> Property
237invErfcIsInverse2 x
238  = (not $ isInfinite x') ==>
239    (y > m_tiny)          ==>
240    (x /= 0) ==>
241    counterexample ("x        = " ++ show x )
242  $ counterexample ("y        = " ++ show y )
243  $ counterexample ("x'       = " ++ show x')
244  $ counterexample ("calc.err = " ++ show delta)
245  $ counterexample ("ulps     = " ++ show (ulpDistance x x'))
246  $ ulpDistance x x' <= delta
247  where
248    (e,e') = roundtrip_invErfc_erfc
249    delta  = round
250           $ e' + e * abs (y / x  /  (2 / sqrt pi * exp( -x*x )))
251    y  = erfc x
252    x' = invErfc y
253
254-- id ≈ erf . invErf
255invErfIsInverse :: Double -> Property
256invErfIsInverse a
257  = (x /= 0) ==>
258    counterexample ("x        = " ++ show x )
259  $ counterexample ("y        = " ++ show y )
260  $ counterexample ("x'       = " ++ show x')
261  $ counterexample ("calc.err = " ++ show delta)
262  $ counterexample ("ulps     = " ++ show (ulpDistance x x'))
263  $ ulpDistance x x' <= delta
264  where
265    (e,e') = roundtrip_erf_invErf
266    delta  = round
267           $ e + e' * abs (y / x  *  2 / sqrt pi * exp ( -y * y ))
268    x  | a < 0     = - range01 a
269       | otherwise =   range01 a
270    y  = invErf x
271    x' = erf y
272
273----------------------------------------------------------------
274-- QC tests
275----------------------------------------------------------------
276
277-- B(p,q) = (x - 1)!(y-1)! / (x + y - 1)!
278betaFactorial :: IO ()
279betaFactorial = do
280  forM_ prod $ \(p,q,facP,facQ,facProd) -> do
281    let exact = fromIntegral (facQ * facP)
282              / fromIntegral facProd
283    checkTabular 16 (show (p,q))
284      (logBeta (fromIntegral p) (fromIntegral q))
285      (log exact)
286  where
287    prod    = [ (p,q,facP,facQ, factorial' (p + q - 1))
288              | (p,facP) <- facList
289              , (q,facQ) <- facList
290              , p + q < 170
291              , not (p == 1 && q== 1)
292              ]
293    facList = [(p,factorial' (p-1)) | p <- [1 .. 170]]
294
295-- B(1,p) = 1/p
296beta1p :: Double -> Property
297beta1p (abs -> p)
298  = p > 2 ==>
299    counterexample ("p    = " ++ show p)
300  $ counterexample ("logB = " ++ show lb)
301  $ counterexample ("err  = " ++ show d)
302  $ d <= 24
303  where
304    lb = logBeta 1 p
305    d  = ulpDistance lb (- log p)
306
307{-
308-- B(p+1,q) = B(p,q) · p/(p+q)
309betaRecurrence :: Double -> Double -> Property
310betaRecurrence (abs -> p) (abs -> q)
311  = p > 0  &&  q > 0  ==>
312    counterexample ("p          = " ++ show p)
313  $ counterexample ("q          = " ++ show q)
314  $ counterexample ("log B(p,q) = " ++ show (logBeta p q))
315  $ counterexample ("log B(p+1,q) = " ++ show (logBeta (p+1) q))
316  $ counterexample ("err        = " ++ show d)
317  $ d <= 128
318  where
319    logB  = logBeta p q + log (p / (p + q))
320    logB' = logBeta (p + 1) q
321    d     = ulpDistance logB logB'
322-}
323
324-- Γ(x+1) = x·Γ(x)
325gammaReccurence :: Double -> Property
326gammaReccurence x
327  = x > 0  ==>  err < errEst
328    where
329      g1     = logGamma x
330      g2     = logGamma (x+1)
331      err    = abs (g2 - g1 - log x)
332      -- logGamma apparently is not as precise for small x. See #59 for details
333      errEst = max 1e-14
334             $ 2 * m_epsilon * sum (map abs [ g1 , g2 , log x ])
335
336-- γ(s,x) is in [0,1] range
337incompleteGammaInRange :: Double -> Double -> Property
338incompleteGammaInRange (abs -> s) (abs -> x) =
339  x >= 0 && s > 0  ==> let i = incompleteGamma s x in i >= 0 && i <= 1
340
341-- γ(1,x) = 1 - exp(-x)
342-- Since Γ(1) = 1 normalization doesn't make any difference
343incompleteGammaAt1Check :: Double -> Bool
344incompleteGammaAt1Check (abs -> x) =
345  ulpDistance (incompleteGamma 1 x) (-expm1(-x)) < 16
346
347-- invIncompleteGamma is inverse of incompleteGamma
348invIGammaIsInverse :: Double -> Double -> Property
349invIGammaIsInverse (abs -> a) (range01 -> p) =
350  a > m_tiny && p > m_tiny && p < 1 && x > m_tiny  ==>
351    ( counterexample ("a    = " ++ show a )
352    $ counterexample ("p    = " ++ show p )
353    $ counterexample ("x    = " ++ show x )
354    $ counterexample ("p'   = " ++ show p')
355    $ counterexample ("err  = " ++ show (ulpDistance p p'))
356    $ counterexample ("est  = " ++ show est)
357    $ ulpDistance p p' <= est
358    )
359  where
360    x  = invIncompleteGamma a p
361    f' = exp ( log x * (a-1) - x - logGamma a)
362    p' = incompleteGamma    a x
363    -- FIXME: Test should be rechecked when #42 is fixed
364    (e,e') = (32,32)
365    est    = round
366           $ e' + e * abs (x / p * f')
367
368-- I(x;p,q) is in [0,1] range
369incompleteBetaInRange :: Double -> Double -> Double -> Property
370incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) =
371  p > 0 && q > 0  ==> let i = incompleteBeta p q x in i >= 0 && i <= 1
372
373-- I(0.5; p,p) = 0.5
374incompleteBetaSymmetry :: Double -> Property
375incompleteBetaSymmetry (abs -> p) =
376    p > 0 ==>
377    counterexample ("p   = " ++ show p)
378  $ counterexample ("ib  = " ++ show ib)
379  $ counterexample ("err = " ++ show d)
380  $ counterexample ("est = " ++ show est)
381  $ d <= est
382  where
383    est | p < 1     = 80
384        | p < 10    = 200
385        | otherwise = round $ 6 * p
386    d  = ulpDistance ib 0.5
387    ib = incompleteBeta p p 0.5
388
389-- invIncompleteBeta is inverse of incompleteBeta
390invIBetaIsInverse :: Double -> Double -> Double -> Property
391invIBetaIsInverse (abs -> p) (abs -> q) (range01 -> x) =
392  p > 0 && q > 0  ==> ( counterexample ("p   = " ++ show p )
393                      $ counterexample ("q   = " ++ show q )
394                      $ counterexample ("x   = " ++ show x )
395                      $ counterexample ("x'  = " ++ show x')
396                      $ counterexample ("a   = " ++ show a)
397                      $ counterexample ("err = " ++ (show $ abs $ (x - x') / x))
398                      $ abs (x - x') <= 1e-12
399                      )
400  where
401    x' = incompleteBeta    p q a
402    a  = invIncompleteBeta p q x
403
404-- Table for digamma function:
405--
406-- Uses equality ψ(n) = H_{n-1} - γ where
407--   H_{n} = Σ 1/k, k = [1 .. n]     - harmonic number
408--   γ     = 0.57721566490153286060  - Euler-Mascheroni number
409digammaTestIntegers :: Double -> Bool
410digammaTestIntegers eps
411  = all (uncurry $ eq eps) $ take 3000 digammaInt
412  where
413    ok approx exact = approx
414    -- Harmonic numbers starting from 0
415    harmN = scanl (\a n -> a + 1/n) 0 [1::Rational .. ]
416    gam   = 0.57721566490153286060
417    -- Digamma values
418    digammaInt = zipWith (\i h -> (digamma i, realToFrac h - gam)) [1..] harmN
419
420
421----------------------------------------------------------------
422-- Unit tests
423----------------------------------------------------------------
424
425-- Lookup table for fact factorial calculation. It has fixed size
426-- which is bad but it's OK for this particular case
427factorial_table :: V.Vector Integer
428factorial_table = V.generate 2000 (\n -> product [1..fromIntegral n])
429
430-- Exact implementation of factorial
431factorial' :: Integer -> Integer
432factorial' n = factorial_table ! fromIntegral n
433
434-- Exact albeit slow implementation of choose
435choose' :: Integer -> Integer -> Integer
436choose' n k = factorial' n `div` (factorial' k * factorial' (n-k))
437
438-- Truncate double to [0,1]
439range01 :: Double -> Double
440range01 = abs . (snd :: (Integer, Double) -> Double) . properFraction
441
442
443readTable :: FilePath -> IO [[Double]]
444readTable
445  = fmap (fmap (fmap read . words) . lines)
446  . readFile
447
448forTable :: FilePath -> ([Double] -> Maybe String) -> IO ()
449forTable path fun = do
450  rows <- readTable path
451  case mapMaybe fun rows of
452    []   -> return ()
453    errs -> assertFailure $ intercalate "---\n" errs
454
455checkTabular :: Int -> String -> Double -> Double -> IO ()
456checkTabular prec x exact val =
457  case checkTabularPure prec x exact val of
458    Nothing -> return ()
459    Just s  -> assertFailure s
460
461checkTabularPure :: Int -> String -> Double -> Double -> Maybe String
462checkTabularPure prec x exact val
463  | within prec exact val = Nothing
464  | otherwise             = Just $ unlines
465      [ " x         = " ++ x
466      , " expected  = " ++ show exact
467      , " got       = " ++ show val
468      , " ulps diff = " ++ show (ulpDistance exact val)
469      , " err.est.  = " ++ show prec
470      ]
471