{-# LANGUAGE ViewPatterns #-} -- | Tests for Statistics.Math module Tests.SpecFunctions ( tests ) where import qualified Data.Vector as V import Data.Vector ((!)) import Test.QuickCheck hiding (choose) import Test.Framework import Test.Framework.Providers.QuickCheck2 import Tests.Helpers import Tests.SpecFunctions.Tables import Numeric.SpecFunctions tests :: Test tests = testGroup "Special functions" [ 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 "γ - increases" \$ \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y , testProperty "invIncompleteGamma = γ^-1" \$ invIGammaIsInverse , testProperty "invIncompleteBeta = B^-1" \$ invIBetaIsInverse -- Unit tests , testAssertion "Factorial is expected to be precise at 1e-15 level" \$ and [ eq 1e-15 (factorial (fromIntegral n)) (fromIntegral (factorial' n)) |n <- [0..170]] , testAssertion "Log factorial is expected to be precise at 1e-15 level" \$ and [ eq 1e-15 (logFactorial (fromIntegral n)) (log \$ fromIntegral \$ factorial' n) | n <- [2..170]] , testAssertion "logGamma is expected to be precise at 1e-9 level [integer points]" \$ and [ eq 1e-9 (logGamma (fromIntegral n)) (logFactorial (n-1)) | n <- [3..10000]] , testAssertion "logGamma is expected to be precise at 1e-9 level [fractional points]" \$ and [ eq 1e-9 (logGamma x) lg | (x,lg) <- tableLogGamma ] , testAssertion "logGammaL is expected to be precise at 1e-15 level" \$ and [ eq 1e-15 (logGammaL (fromIntegral n)) (logFactorial (n-1)) | n <- [3..10000]] , testAssertion "logGammaL is expected to be precise at 1e-9 level [fractional points]" \$ and [ eq 1e-10 (logGammaL x) lg | (x,lg) <- tableLogGamma ] , testAssertion "logBeta is expected to be precise at 1e-6 level" \$ and [ eq 1e-6 (logBeta p q) (logGammaL p + logGammaL q - logGammaL (p+q)) | p <- [0.1,0.2 .. 0.9] ++ [2 .. 20] , q <- [0.1,0.2 .. 0.9] ++ [2 .. 20]] -- FIXME: Why 1e-8? Is it due to poor precision of logBeta? , testAssertion "incompleteBeta is expected to be precise at 1e-8 level" \$ and [ eq 1e-8 (incompleteBeta p q x) ib | (p,q,x,ib) <- tableIncompleteBeta ] , 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]] ] ---------------------------------------------------------------- -- QC tests ---------------------------------------------------------------- -- Γ(x+1) = x·Γ(x) gammaReccurence :: (Double -> Double) -> Double -> Double -> Property gammaReccurence logG ε x = (x > 0 && x < 100) ==> (abs (g2 - g1 - log x) < ε) where g1 = logG x g2 = logG (x+1) -- γ(1,x) = 1 - exp(-x) -- Since Γ(1) = 1 normalization doesn't make any difference incompleteGammaAt1Check :: Double -> Property incompleteGammaAt1Check 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 ) \$ printTestCase ("p' = " ++ show p') \$ printTestCase ("Δp = " ++ show (p - p')) \$ abs (p - p') <= 1e-12 ) where x = invIncompleteGamma a p p' = incompleteGamma a x -- invIncompleteBeta is inverse of incompleteBeta invIBetaIsInverse :: Double -> Double -> Double -> Property invIBetaIsInverse (abs -> p) (abs -> q) (abs . snd . properFraction -> x) = p > 0 && q > 0 ==> ( printTestCase ("p = " ++ show p ) \$ printTestCase ("q = " ++ show q ) \$ printTestCase ("x = " ++ show x ) \$ printTestCase ("x' = " ++ show x') \$ printTestCase ("a = " ++ show a) \$ printTestCase ("err = " ++ (show \$ abs \$ (x - x') / x)) \$ abs (x - x') <= 1e-12 ) where x' = incompleteBeta p q a a = invIncompleteBeta p q x ---------------------------------------------------------------- -- 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))