summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorBryanOSullivan <>2011-11-14 02:15:11 (GMT)
committerhdiff <hdiff@luite.com>2011-11-14 02:15:11 (GMT)
commit3e28fdda1673ca8c81d569008a361ec9b531360e (patch)
treec3216cdee95ba1ec46cf1697de88501161c4bbba
parenteebd81440ef5aab6214cc5c97b83fb0ab6565a56 (diff)
version 0.10.0.10.10.0.1
-rw-r--r--Statistics/Distribution/FDistribution.hs2
-rw-r--r--Statistics/Sample/KernelDensity.hs6
-rw-r--r--Statistics/Transform.hs26
-rw-r--r--statistics.cabal10
-rw-r--r--tests/Tests/Distribution.hs276
-rw-r--r--tests/Tests/Helpers.hs96
-rw-r--r--tests/Tests/Math.hs159
-rw-r--r--tests/Tests/Math/Tables.hs47
-rw-r--r--tests/Tests/Math/gen.py51
-rw-r--r--tests/Tests/NonparametricTest.hs148
-rw-r--r--tests/Tests/Transform.hs82
11 files changed, 892 insertions, 11 deletions
diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs
index dde5e5b..3776d34 100644
--- a/Statistics/Distribution/FDistribution.hs
+++ b/Statistics/Distribution/FDistribution.hs
@@ -25,7 +25,7 @@ import Statistics.Math (logBeta, incompleteBeta, invIncompleteBeta)
-- | Student-T distribution
data FDistribution = F { fDistributionNDF1 :: {-# UNPACK #-} !Double
, fDistributionNDF2 :: {-# UNPACK #-} !Double
- , pdfFactor :: {-# UNPACK #-} !Double
+ , _pdfFactor :: {-# UNPACK #-} !Double
}
deriving (Eq,Show,Read,Typeable)
diff --git a/Statistics/Sample/KernelDensity.hs b/Statistics/Sample/KernelDensity.hs
index 57c01bc..fa274df 100644
--- a/Statistics/Sample/KernelDensity.hs
+++ b/Statistics/Sample/KernelDensity.hs
@@ -31,7 +31,7 @@ import Statistics.Constants (m_sqrt_2_pi)
import Statistics.Function (minMax, nextHighestPowerOfTwo)
import Statistics.Math.RootFinding (fromRoot, ridders)
import Statistics.Sample.Histogram (histogram_)
-import Statistics.Transform (dct, idct)
+import Statistics.Transform (dct_, idct_)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
@@ -78,12 +78,12 @@ kde_ n0 min max xs
where
mesh = G.generate ni $ \z -> min + (d * fromIntegral z)
where d = r / (n-1)
- density = G.map (/r) . idct $ G.zipWith f a (G.enumFromTo 0 (n-1))
+ density = G.map (/r) . idct_ $ G.zipWith f a (G.enumFromTo 0 (n-1))
where f b z = b * exp (sqr z * sqr pi * t_star * (-0.5)) :+ 0
!n = fromIntegral ni
!ni = nextHighestPowerOfTwo n0
!r = max - min
- a = dct . G.map (/ G.sum h) $ h
+ a = dct_ . G.map (/ G.sum h) $ h
where h = G.map (/ (len :+ 0)) $ histogram_ ni min max xs
!len = fromIntegral (G.length xs)
!t_star = fromRoot (0.28 * len ** (-0.4)) . ridders 1e-14 (0,0.1) $ \x ->
diff --git a/Statistics/Transform.hs b/Statistics/Transform.hs
index c9b778d..4a7c4d1 100644
--- a/Statistics/Transform.hs
+++ b/Statistics/Transform.hs
@@ -17,9 +17,14 @@
module Statistics.Transform
(
+ -- * Type synonyms
CD
+ -- * Discrete cosine transform
, dct
+ , dct_
, idct
+ , idct_
+ -- * Fast Fourier transform
, fft
, ifft
) where
@@ -35,9 +40,13 @@ import qualified Data.Vector.Unboxed as U
type CD = Complex Double
--- | Discrete cosine transform.
-dct :: U.Vector CD -> U.Vector Double
-dct xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)
+-- | Discrete cosine transform (DCT-II).
+dct :: U.Vector Double -> U.Vector Double
+dct = dct_ . G.map (:+0)
+
+-- | Discrete cosine transform, with complex coefficients (DCT-II).
+dct_ :: U.Vector CD -> U.Vector Double
+dct_ xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)
where
interleaved = G.backpermute xs $ G.enumFromThenTo 0 2 (len-2) G.++
G.enumFromThenTo (len-1) (len-3) 1
@@ -46,9 +55,14 @@ dct xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)
where n = fi len
len = G.length xs
--- | Inverse discrete cosine transform.
-idct :: U.Vector CD -> U.Vector Double
-idct xs = G.generate len interleave
+-- | Inverse discrete cosine transform (DCT-III).
+idct :: U.Vector Double -> U.Vector Double
+idct = idct_ . G.map (:+0)
+
+-- | Inverse discrete cosine transform, with complex coefficients
+-- (DCT-III).
+idct_ :: U.Vector CD -> U.Vector Double
+idct_ xs = G.generate len interleave
where
interleave z | even z = vals `G.unsafeIndex` halve z
| otherwise = vals `G.unsafeIndex` (len - halve z - 1)
diff --git a/statistics.cabal b/statistics.cabal
index a9bda03..4549048 100644
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
name: statistics
-version: 0.10.0.0
+version: 0.10.0.1
synopsis: A library of statistical types, data, and functions
description:
This library provides a number of common functions and types useful
@@ -111,6 +111,7 @@ extra-source-files:
examples/kde/data/faithful.csv
examples/kde/kde.html
examples/kde/kde.tpl
+ tests/Tests/Math/gen.py
library
exposed-modules:
@@ -175,6 +176,13 @@ test-suite tests
type: exitcode-stdio-1.0
hs-source-dirs: tests
main-is: tests.hs
+ other-modules:
+ Tests.Distribution
+ Tests.Helpers
+ Tests.Math
+ Tests.Math.Tables
+ Tests.NonparametricTest
+ Tests.Transform
ghc-options:
-Wall -threaded -rtsopts
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
new file mode 100644
index 0000000..51490b9
--- /dev/null
+++ b/tests/Tests/Distribution.hs
@@ -0,0 +1,276 @@
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+-- Required for Param
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE OverlappingInstances #-}
+module Tests.Distribution (
+ distributionTests
+ ) where
+
+import Control.Applicative
+import Control.Exception
+
+import Data.List (find)
+import Data.Typeable (Typeable)
+
+import qualified Numeric.IEEE as IEEE
+
+import Test.Framework (Test,testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck as QC
+import Test.QuickCheck.Monadic as QC
+import Text.Printf
+
+import Statistics.Distribution
+import Statistics.Distribution.Binomial
+import Statistics.Distribution.ChiSquared
+import Statistics.Distribution.CauchyLorentz
+import Statistics.Distribution.Exponential
+import Statistics.Distribution.FDistribution
+import Statistics.Distribution.Gamma
+import Statistics.Distribution.Geometric
+import Statistics.Distribution.Hypergeometric
+import Statistics.Distribution.Normal
+import Statistics.Distribution.Poisson
+import Statistics.Distribution.StudentT
+import Statistics.Distribution.Uniform
+
+import Prelude hiding (catch)
+
+import Tests.Helpers
+
+
+-- | Tests for all distributions
+distributionTests :: Test
+distributionTests = testGroup "Tests for all distributions"
+ [ contDistrTests (T :: T CauchyDistribution )
+ , contDistrTests (T :: T ChiSquared )
+ , contDistrTests (T :: T ExponentialDistribution )
+ , contDistrTests (T :: T GammaDistribution )
+ , contDistrTests (T :: T NormalDistribution )
+ , contDistrTests (T :: T UniformDistribution )
+ , contDistrTests (T :: T StudentT )
+ , contDistrTests (T :: T FDistribution )
+
+ , discreteDistrTests (T :: T BinomialDistribution )
+ , discreteDistrTests (T :: T GeometricDistribution )
+ , discreteDistrTests (T :: T HypergeometricDistribution )
+ , discreteDistrTests (T :: T PoissonDistribution )
+
+ , unitTests
+ ]
+
+----------------------------------------------------------------
+-- Tests
+----------------------------------------------------------------
+
+-- Tests for continous distribution
+contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
+contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
+ cdfTests t ++
+ [ testProperty "PDF sanity" $ pdfSanityCheck t
+ , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t
+ , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
+ ]
+
+-- Tests for discrete distribution
+discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
+discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
+ cdfTests t ++
+ [ testProperty "Prob. sanity" $ probSanityCheck t
+ , testProperty "CDF is sum of prob." $ discreteCDFcorrect t
+ ]
+
+-- Tests for distributions which have CDF
+cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d) => T d -> [Test]
+cdfTests t =
+ [ testProperty "C.D.F. sanity" $ cdfSanityCheck t
+ , testProperty "CDF limit at +∞" $ cdfLimitAtPosInfinity t
+ , testProperty "CDF limit at -∞" $ cdfLimitAtNegInfinity t
+ , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing t
+ , testProperty "1-CDF is correct" $ cdfComplementIsCorrect t
+ ]
+----------------------------------------------------------------
+
+-- CDF is in [0,1] range
+cdfSanityCheck :: (Distribution d) => T d -> d -> Double -> Bool
+cdfSanityCheck _ d x = c >= 0 && c <= 1
+ where c = cumulative d x
+
+-- CDF never decreases
+cdfIsNondecreasing :: (Distribution d) => T d -> d -> Double -> Double -> Bool
+cdfIsNondecreasing _ d = monotonicallyIncreasesIEEE $ cumulative d
+
+-- CDF limit at +∞ is 1
+cdfLimitAtPosInfinity :: (Param d, Distribution d) => T d -> d -> Property
+cdfLimitAtPosInfinity _ d =
+ okForInfLimit d ==> printTestCase ("Last elements: " ++ show (drop 990 probs))
+ $ Just 1.0 == (find (>=1) probs)
+ where
+ probs = take 1000 $ map (cumulative d) $ iterate (*1.4) 1
+
+-- CDF limit at -∞ is 0
+cdfLimitAtNegInfinity :: (Param d, Distribution d) => T d -> d -> Property
+cdfLimitAtNegInfinity _ d =
+ okForInfLimit d ==> printTestCase ("Last elements: " ++ show (drop 990 probs))
+ $ case find (< IEEE.epsilon) probs of
+ Nothing -> False
+ Just p -> p >= 0
+ where
+ probs = take 1000 $ map (cumulative d) $ iterate (*1.4) (-1)
+
+-- CDF's complement is implemented correctly
+cdfComplementIsCorrect :: (Distribution d) => T d -> d -> Double -> Bool
+cdfComplementIsCorrect _ d x = (eq 1e-14) 1 (cumulative d x + complCumulative d x)
+
+
+-- PDF is positive
+pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
+pdfSanityCheck _ d x = p >= 0
+ where p = density d x
+
+-- Quantile is inverse of CDF
+quantileIsInvCDF :: (Param d, ContDistr d) => T d -> d -> Double -> Property
+quantileIsInvCDF _ d p =
+ p > 0 && p < 1 ==> ( printTestCase (printf "Quantile = %g" q )
+ $ printTestCase (printf "Probability = %g" p )
+ $ printTestCase (printf "Probability' = %g" p')
+ $ printTestCase (printf "Error = %e" (abs $ p - p'))
+ $ abs (p - p') < invQuantilePrec d
+ )
+ where
+ q = quantile d p
+ p' = cumulative d q
+
+-- Test that quantile fails if p<0 or p>1
+quantileShouldFail :: (ContDistr d) => T d -> d -> Double -> Property
+quantileShouldFail _ d p =
+ p < 0 || p > 1 ==> QC.monadicIO $ do r <- QC.run $ catch
+ (do { return $! quantile d p; return False })
+ (\(e :: SomeException) -> return True)
+ QC.assert r
+
+
+-- Probability is in [0,1] range
+probSanityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Bool
+probSanityCheck _ d x = p >= 0 && p <= 1
+ where p = probability d x
+
+-- Check that discrete CDF is correct
+discreteCDFcorrect :: (DiscreteDistr d) => T d -> d -> Int -> Int -> Property
+discreteCDFcorrect _ d a b =
+ abs (a - b) < 100 ==> abs (p1 - p2) < 3e-10
+ -- Avoid too large differeneces. Otherwise there is to much to sum
+ --
+ -- Absolute difference is used guard againist precision loss when
+ -- close values of CDF are subtracted
+ where
+ n = min a b
+ m = max a b
+ p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
+ p2 = sum $ map (probability d) [n .. m]
+
+
+
+----------------------------------------------------------------
+-- Arbitrary instances for ditributions
+----------------------------------------------------------------
+
+instance QC.Arbitrary BinomialDistribution where
+ arbitrary = binomial <$> QC.choose (1,100) <*> QC.choose (0,1)
+instance QC.Arbitrary ExponentialDistribution where
+ arbitrary = exponential <$> QC.choose (0,100)
+instance QC.Arbitrary GammaDistribution where
+ arbitrary = gammaDistr <$> QC.choose (0.1,10) <*> QC.choose (0.1,10)
+instance QC.Arbitrary GeometricDistribution where
+ arbitrary = geometric <$> QC.choose (0,1)
+instance QC.Arbitrary HypergeometricDistribution where
+ arbitrary = do l <- QC.choose (1,20)
+ m <- QC.choose (0,l)
+ k <- QC.choose (1,l)
+ return $ hypergeometric m l k
+instance QC.Arbitrary NormalDistribution where
+ arbitrary = normalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-3, 1e3)
+instance QC.Arbitrary PoissonDistribution where
+ arbitrary = poisson <$> QC.choose (0,1)
+instance QC.Arbitrary ChiSquared where
+ arbitrary = chiSquared <$> QC.choose (1,100)
+instance QC.Arbitrary UniformDistribution where
+ arbitrary = do a <- QC.arbitrary
+ b <- QC.arbitrary `suchThat` (/= a)
+ return $ uniformDistr a b
+instance QC.Arbitrary CauchyDistribution where
+ arbitrary = cauchyDistribution
+ <$> arbitrary
+ <*> ((abs <$> arbitrary) `suchThat` (> 0))
+instance QC.Arbitrary StudentT where
+ arbitrary = studentT <$> ((abs <$> arbitrary) `suchThat` (>0))
+instance QC.Arbitrary FDistribution where
+ arbitrary = fDistribution
+ <$> ((abs <$> arbitrary) `suchThat` (>0))
+ <*> ((abs <$> arbitrary) `suchThat` (>0))
+
+
+
+-- Parameters for distribution testing. Some distribution require
+-- relaxing parameters a bit
+class Param a where
+ -- Precision for quantileIsInvCDF
+ invQuantilePrec :: a -> Double
+ invQuantilePrec _ = 1e-14
+ -- Distribution is OK for testing limits
+ okForInfLimit :: a -> Bool
+ okForInfLimit _ = True
+
+
+instance Param a
+
+instance Param StudentT where
+ invQuantilePrec _ = 1e-13
+ okForInfLimit d = studentTndf d > 0.75
+
+instance Param FDistribution where
+ invQuantilePrec _ = 1e-12
+
+
+
+----------------------------------------------------------------
+-- Unit tests
+----------------------------------------------------------------
+
+unitTests :: Test
+unitTests = testGroup "Unit tests"
+ [ testAssertion "density (gammaDistr 150 1/150) 1 == 4.883311" $
+ 4.883311418525483 =~ (density (gammaDistr 150 (1/150)) 1)
+ -- Student-T
+ , testStudentPDF 0.3 1.34 0.0648215 -- PDF
+ , testStudentPDF 1 0.42 0.27058
+ , testStudentPDF 4.4 0.33 0.352994
+ , testStudentCDF 0.3 3.34 0.757146 -- CDF
+ , testStudentCDF 1 0.42 0.626569
+ , testStudentCDF 4.4 0.33 0.621739
+ -- F-distribution
+ , testFdistrPDF 1 3 3 (1/(6 * pi)) -- PDF
+ , testFdistrPDF 2 2 1.2 0.206612
+ , testFdistrPDF 10 12 8 0.000385613179281892790166
+ , testFdistrCDF 1 3 3 0.81830988618379067153 -- CDF
+ , testFdistrCDF 2 2 1.2 0.545455
+ , testFdistrCDF 10 12 8 0.99935509863451408041
+ ]
+ where
+ -- Student-T
+ testStudentPDF ndf x exact
+ = testAssertion (printf "density (studentT %f) %f ≈ %f" ndf x exact)
+ $ eq 1e-5 exact (density (studentT ndf) x)
+ testStudentCDF ndf x exact
+ = testAssertion (printf "cumulative (studentT %f) %f ≈ %f" ndf x exact)
+ $ eq 1e-5 exact (cumulative (studentT ndf) x)
+ -- F-distribution
+ testFdistrPDF n m x exact
+ = testAssertion (printf "density (fDistribution %i %i) %f ≈ %f [got %f]" n m x exact d)
+ $ eq 1e-5 exact d
+ where d = density (fDistribution n m) x
+ testFdistrCDF n m x exact
+ = testAssertion (printf "cumulative (fDistribution %i %i) %f ≈ %f [got %f]" n m x exact d)
+ $ eq 1e-5 exact d
+ where d = cumulative (fDistribution n m) x
diff --git a/tests/Tests/Helpers.hs b/tests/Tests/Helpers.hs
new file mode 100644
index 0000000..e5cd0c5
--- /dev/null
+++ b/tests/Tests/Helpers.hs
@@ -0,0 +1,96 @@
+-- | Helpers for testing
+module Tests.Helpers (
+ -- * helpers
+ T(..)
+ , typeName
+ , eq
+ , eqC
+ , (=~)
+ -- * Generic QC tests
+ , monotonicallyIncreases
+ , monotonicallyIncreasesIEEE
+ -- * HUnit helpers
+ , testAssertion
+ ) where
+
+import Data.Complex
+import Data.Typeable
+
+import qualified Numeric.IEEE as IEEE
+
+import qualified Test.HUnit as HU
+import Test.Framework
+import Test.Framework.Providers.HUnit
+
+import Statistics.Constants
+
+
+
+----------------------------------------------------------------
+-- Helpers
+----------------------------------------------------------------
+
+-- | Phantom typed value used to select right instance in QC tests
+data T a = T
+
+-- | String representation of type name
+typeName :: Typeable a => T a -> String
+typeName = show . typeOf . typeParam
+ where
+ typeParam :: T a -> a
+ typeParam _ = undefined
+
+-- | Approximate equality for 'Double'. Doesn't work well for numbers
+-- which are almost zero.
+eq :: Double -- ^ Relative error
+ -> Double -> Double -> Bool
+eq eps a b
+ | a == 0 && b == 0 = True
+ | otherwise = abs (a - b) <= eps * max (abs a) (abs b)
+
+-- | Approximate equality for 'Complex Double'
+eqC :: Double -- ^ Relative error
+ -> Complex Double
+ -> Complex Double
+ -> Bool
+eqC eps a@(ar :+ ai) b@(br :+ bi)
+ | a == 0 && b == 0 = True
+ | otherwise = abs (ar - br) <= eps * d
+ && abs (ai - bi) <= eps * d
+ where
+ d = max (realPart $ abs a) (realPart $ abs b)
+
+
+-- | Approximately equal up to 1 ulp
+(=~) :: Double -> Double -> Bool
+(=~) = eq m_epsilon
+
+
+----------------------------------------------------------------
+-- Generic QC
+----------------------------------------------------------------
+
+-- Check that function is nondecreasing
+monotonicallyIncreases :: (Ord a, Ord b) => (a -> b) -> a -> a -> Bool
+monotonicallyIncreases f x1 x2 = f (min x1 x2) <= f (max x1 x2)
+
+-- Check that function is nondecreasing taking rounding errors into
+-- account.
+--
+-- In fact funstion is allowed to decrease less than one ulp in order
+-- to guard againist problems with excess precision. On x86 FPU works
+-- with 80-bit numbers but doubles are 64-bit so rounding happens
+-- whenever values are moved from registers to memory
+monotonicallyIncreasesIEEE :: (Ord a, IEEE.IEEE b) => (a -> b) -> a -> a -> Bool
+monotonicallyIncreasesIEEE f x1 x2 =
+ y1 <= y2 || (y1 - y2) < y2 * IEEE.epsilon
+ where
+ y1 = f (min x1 x2)
+ y2 = f (max x1 x2)
+
+----------------------------------------------------------------
+-- HUnit helpers
+----------------------------------------------------------------
+
+testAssertion :: String -> Bool -> Test
+testAssertion str cont = testCase str $ HU.assertBool str cont
diff --git a/tests/Tests/Math.hs b/tests/Tests/Math.hs
new file mode 100644
index 0000000..df35fc7
--- /dev/null
+++ b/tests/Tests/Math.hs
@@ -0,0 +1,159 @@
+{-# LANGUAGE ViewPatterns #-}
+-- | Tests for Statistics.Math
+module Tests.Math (
+ mathTests
+ ) where
+
+import Data.Vector.Unboxed (fromList)
+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.Math.Tables
+import Statistics.Math
+
+
+mathTests :: Test
+mathTests = testGroup "S.Math"
+ [ 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
+ , chebyshevTests
+ -- 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
+
+-- Test that Chebyshev polynomial of low order are evaluated correctly
+chebyshevTests :: Test
+chebyshevTests = testGroup "Chebyshev polynomials"
+ [ testProperty "Chebyshev 0" $ \a0 (Ch x) ->
+ (ch0 x * a0) ≈ (chebyshev x $ fromList [a0])
+ , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->
+ (a0*ch0 x + a1*ch1 x) ≈ (chebyshev x $ fromList [a0,a1])
+ , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->
+ (a0*ch0 x + a1*ch1 x + a2*ch2 x) ≈ (chebyshev x $ fromList [a0,a1,a2])
+ , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->
+ (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3])
+ , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->
+ (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x + a4*ch4 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3,a4])
+ ]
+ where (≈) = eq 1e-12
+
+-- Chebyshev polynomials of low order
+ch0,ch1,ch2,ch3,ch4 :: Double -> Double
+ch0 _ = 1
+ch1 x = x
+ch2 x = 2*x^2 - 1
+ch3 x = 4*x^3 - 3*x
+ch4 x = 8*x^4 - 8*x^2 + 1
+
+newtype Ch = Ch Double
+ deriving Show
+instance Arbitrary Ch where
+ arbitrary = do x <- arbitrary
+ return $ Ch $ 2 * (snd . properFraction) x - 1
+
+
+
+----------------------------------------------------------------
+-- 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))
diff --git a/tests/Tests/Math/Tables.hs b/tests/Tests/Math/Tables.hs
new file mode 100644
index 0000000..89ac560
--- /dev/null
+++ b/tests/Tests/Math/Tables.hs
@@ -0,0 +1,47 @@
+module Tests.Math.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)
+ ]
diff --git a/tests/Tests/Math/gen.py b/tests/Tests/Math/gen.py
new file mode 100644
index 0000000..fde9da6
--- /dev/null
+++ b/tests/Tests/Math/gen.py
@@ -0,0 +1,51 @@
+#!/usr/bin/python
+"""
+"""
+
+from mpmath import *
+
+def printListLiteral(lines) :
+ print " [" + "\n , ".join(lines) + "\n ]"
+
+################################################################
+# Generate header
+print "module Tests.Math.Tables where"
+print
+
+################################################################
+## Generate table for logGamma
+print "tableLogGamma :: [(Double,Double)]"
+print "tableLogGamma ="
+
+gammaArg = [ 1.25e-6, 6.82e-5, 2.46e-4, 8.8e-4, 3.12e-3, 2.67e-2,
+ 7.77e-2, 0.234, 0.86, 1.34, 1.89, 2.45,
+ 3.65, 4.56, 6.66, 8.25, 11.3, 25.6,
+ 50.4, 123.3, 487.4, 853.4, 2923.3, 8764.3,
+ 1.263e4, 3.45e4, 8.234e4, 2.348e5, 8.343e5, 1.23e6,
+ ]
+printListLiteral(
+ [ '(%.15f, %.20g)' % (x, log(gamma(x))) for x in gammaArg ]
+ )
+
+
+################################################################
+## Generate table for incompleteBeta
+
+print "tableIncompleteBeta :: [(Double,Double,Double,Double)]"
+print "tableIncompleteBeta ="
+
+incompleteBetaArg = [
+ (2, 3, 0.03),
+ (2, 3, 0.23),
+ (2, 3, 0.76),
+ (4, 2.3, 0.89),
+ (1, 1, 0.55),
+ (0.3, 12.2, 0.11),
+ (13.1, 9.8, 0.12),
+ (13.1, 9.8, 0.42),
+ (13.1, 9.8, 0.92),
+ ]
+printListLiteral(
+ [ '(%.15f, %.15f, %.15f, %.20g)' % (p,q,x, betainc(p,q,0,x, regularized=True))
+ for (p,q,x) in incompleteBetaArg
+ ])
diff --git a/tests/Tests/NonparametricTest.hs b/tests/Tests/NonparametricTest.hs
new file mode 100644
index 0000000..671f97b
--- /dev/null
+++ b/tests/Tests/NonparametricTest.hs
@@ -0,0 +1,148 @@
+-- Tests for Statistics.Test.NonParametric
+module Tests.NonparametricTest (
+ nonparametricTests
+ ) where
+
+
+import qualified Data.Vector.Unboxed as U
+import Test.HUnit (Test(..),assertEqual,assertBool)
+import qualified Test.Framework as TF
+import Test.Framework.Providers.HUnit
+
+import Statistics.Test.MannWhitneyU
+import Statistics.Test.WilcoxonT
+
+
+
+
+nonparametricTests :: TF.Test
+nonparametricTests = TF.testGroup "Nonparametric tests"
+ $ hUnitTestToTests =<< concat [ mannWhitneyTests
+ , wilcoxonSumTests
+ , wilcoxonPairTests
+ ]
+
+
+----------------------------------------------------------------
+
+mannWhitneyTests :: [Test]
+mannWhitneyTests = zipWith test [(0::Int)..] testData ++
+ [TestCase $ assertEqual "Mann-Whitney U Critical Values, m=1"
+ (replicate (20*3) Nothing)
+ [mannWhitneyUCriticalValue (1,x) p | x <- [1..20], p <- [0.005,0.01,0.025]]
+ ,TestCase $ assertEqual "Mann-Whitney U Critical Values, m=2, p=0.025"
+ (replicate 7 Nothing ++ map Just [0,0,0,0,1,1,1,1,1,2,2,2,2])
+ [mannWhitneyUCriticalValue (2,x) 0.025 | x <- [1..20]]
+ ,TestCase $ assertEqual "Mann-Whitney U Critical Values, m=6, p=0.05"
+ (replicate 1 Nothing ++ map Just [0, 2,3,5,7,8,10,12,14,16,17,19,21,23,25,26,28,30,32])
+ [mannWhitneyUCriticalValue (6,x) 0.05 | x <- [1..20]]
+ ,TestCase $ assertEqual "Mann-Whitney U Critical Values, m=20, p=0.025"
+ (replicate 1 Nothing ++ map Just [2,8,14,20,27,34,41,48,55,62,69,76,83,90,98,105,112,119,127])
+ [mannWhitneyUCriticalValue (20,x) 0.025 | x <- [1..20]]
+ ]
+ where
+ test n (a, b, c, d)
+ = TestCase $ do assertEqual ("Mann-Whitney U " ++ show n) c us
+ assertEqual ("Mann-Whitney U Sig " ++ show n)
+ d $ mannWhitneyUSignificant TwoTailed (length a, length b) 0.05 us
+ where
+ us = mannWhitneyU (U.fromList a) (U.fromList b)
+
+ -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
+ testData :: [([Double], [Double], (Double, Double), Maybe TestResult)]
+ testData = [ ( [3,4,2,6,2,5]
+ , [9,7,5,10,6,8]
+ , (2, 34)
+ , Just Significant
+ )
+ , ( [540,480,600,590,605]
+ , [760,890,1105,595,940]
+ , (2, 23)
+ , Just Significant
+ )
+ , ( [19,22,16,29,24]
+ , [20,11,17,12]
+ , (17, 3)
+ , Just NotSignificant
+ )
+ , ( [126,148,85,61, 179,93, 45,189,85,93]
+ , [194,128,69,135,171,149,89,248,79,137]
+ , (35,65)
+ , Just NotSignificant
+ )
+ , ( [1..30]
+ , [1..30]
+ , (450,450)
+ , Just NotSignificant
+ )
+ , ( [1 .. 30]
+ , [11.5 .. 40 ]
+ , (190.0,710.0)
+ , Just Significant
+ )
+ ]
+
+wilcoxonSumTests :: [Test]
+wilcoxonSumTests = zipWith test [(0::Int)..] testData
+ where
+ test n (a, b, c) = TestCase $ assertEqual ("Wilcoxon Sum " ++ show n) c (wilcoxonRankSums (U.fromList a) (U.fromList b))
+
+ -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
+ testData :: [([Double], [Double], (Double, Double))]
+ testData = [ ( [8.50,9.48,8.65,8.16,8.83,7.76,8.63]
+ , [8.27,8.20,8.25,8.14,9.00,8.10,7.20,8.32,7.70]
+ , (75, 61)
+ )
+ , ( [0.45,0.50,0.61,0.63,0.75,0.85,0.93]
+ , [0.44,0.45,0.52,0.53,0.56,0.58,0.58,0.65,0.79]
+ , (71.5, 64.5)
+ )
+ ]
+
+wilcoxonPairTests :: [Test]
+wilcoxonPairTests = zipWith test [(0::Int)..] testData ++
+ -- Taken from the Mitic paper:
+ [ TestCase $ assertBool "Sig 16, 35" (to4dp 0.0467 $ wilcoxonMatchedPairSignificance 16 35)
+ , TestCase $ assertBool "Sig 16, 36" (to4dp 0.0523 $ wilcoxonMatchedPairSignificance 16 36)
+ , TestCase $ assertEqual "Wilcoxon critical values, p=0.05"
+ (replicate 4 Nothing ++ map Just [0,2,3,5,8,10,13,17,21,25,30,35,41,47,53,60,67,75,83,91,100,110,119])
+ [wilcoxonMatchedPairCriticalValue x 0.05 | x <- [1..27]]
+ , TestCase $ assertEqual "Wilcoxon critical values, p=0.025"
+ (replicate 5 Nothing ++ map Just [0,2,3,5,8,10,13,17,21,25,29,34,40,46,52,58,65,73,81,89,98,107])
+ [wilcoxonMatchedPairCriticalValue x 0.025 | x <- [1..27]]
+ , TestCase $ assertEqual "Wilcoxon critical values, p=0.01"
+ (replicate 6 Nothing ++ map Just [0,1,3,5,7,9,12,15,19,23,27,32,37,43,49,55,62,69,76,84,92])
+ [wilcoxonMatchedPairCriticalValue x 0.01 | x <- [1..27]]
+ , TestCase $ assertEqual "Wilcoxon critical values, p=0.005"
+ (replicate 7 Nothing ++ map Just [0,1,3,5,7,9,12,15,19,23,27,32,37,42,48,54,61,68,75,83])
+ [wilcoxonMatchedPairCriticalValue x 0.005 | x <- [1..27]]
+ ]
+ where
+ test n (a, b, c) = TestCase $ assertEqual ("Wilcoxon Paired " ++ show n) c (wilcoxonMatchedPairSignedRank (U.fromList a) (U.fromList b))
+
+ -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
+ testData :: [([Double], [Double], (Double, Double))]
+ testData = [ ([1..10], [1..10], (0, 0 ))
+ , ([1..5], [6..10], (0, 5*(-3)))
+ -- Worked example from the Internet:
+ , ( [125,115,130,140,140,115,140,125,140,135]
+ , [110,122,125,120,140,124,123,137,135,145]
+ , ( sum $ filter (> 0) [7,-3,1.5,9,0,-4,8,-6,1.5,-5]
+ , sum $ filter (< 0) [7,-3,1.5,9,0,-4,8,-6,1.5,-5]
+ )
+ )
+ -- Worked examples from books/papers:
+ , ( [2.4,1.9,2.3,1.9,2.4,2.5]
+ , [2.0,2.1,2.0,2.0,1.8,2.0]
+ , (18, -3)
+ )
+ , ( [130,170,125,170,130,130,145,160]
+ , [120,163,120,135,143,136,144,120]
+ , (27, -9)
+ )
+ , ( [540,580,600,680,430,740,600,690,605,520]
+ , [760,710,1105,880,500,990,1050,640,595,520]
+ , (3, -42)
+ )
+ ]
+ to4dp tgt x = x >= tgt - 0.00005 && x < tgt + 0.00005
diff --git a/tests/Tests/Transform.hs b/tests/Tests/Transform.hs
new file mode 100644
index 0000000..7021579
--- /dev/null
+++ b/tests/Tests/Transform.hs
@@ -0,0 +1,82 @@
+module Tests.Transform
+ (
+ tests
+ ) where
+
+import Data.Bits ((.&.), shiftL)
+import Data.Complex (Complex((:+)))
+import Data.Functor ((<$>))
+import Statistics.Function (within)
+import Statistics.Transform
+
+import Test.Framework (Test, testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck (Positive(..),Property,choose,vectorOf,
+ arbitrary,printTestCase)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+import Tests.Helpers
+
+
+
+tests :: Test
+tests = testGroup "fft" [
+ testProperty "t_impulse" t_impulse
+ , testProperty "t_impulse_offset" t_impulse_offset
+ , testProperty "ifft . fft = id" (t_fftInverse $ ifft . fft)
+ , testProperty "fft . ifft = id" (t_fftInverse $ fft . ifft)
+ ]
+
+-- A single real-valued impulse at the beginning of an otherwise zero
+-- vector should be replicated in every real component of the result,
+-- and all the imaginary components should be zero.
+t_impulse :: Double -> Positive Int -> Bool
+t_impulse k (Positive m) = G.all (c_near i) (fft v)
+ where v = i `G.cons` G.replicate (n-1) 0
+ i = k :+ 0
+ n = 1 `shiftL` (m .&. 6)
+
+-- If a real-valued impulse is offset from the beginning of an
+-- otherwise zero vector, the sum-of-squares of each component of the
+-- result should equal the square of the impulse.
+t_impulse_offset :: Double -> Positive Int -> Positive Int -> Bool
+t_impulse_offset k (Positive x) (Positive m) = G.all ok (fft v)
+ where v = G.concat [G.replicate xn 0, G.singleton i, G.replicate (n-xn-1) 0]
+ ok (re :+ im) = within ulps (re*re + im*im) (k*k)
+ i = k :+ 0
+ xn = x `rem` n
+ n = 1 `shiftL` (m .&. 6)
+
+-- Test that (ifft . fft ≈ id)
+--
+-- Approximate equality here is tricky. Smaller values of vector tend
+-- to have large relative error. Thus we should test that vectors as
+-- whole are approximate equal.
+t_fftInverse :: (U.Vector CD -> U.Vector CD) -> Property
+t_fftInverse roundtrip = do
+ n <- (2^) <$> choose (0,9::Int) -- Size of vector
+ x <- G.fromList <$> vectorOf n arbitrary -- Vector to transform
+ let x' = roundtrip x
+ id $ printTestCase "Original vector"
+ $ printTestCase (show x)
+ $ printTestCase "Transformed one"
+ $ printTestCase (show x)
+ $ printTestCase (show n)
+ $ vectorNorm (U.zipWith (-) x x') <= 1e-15 * vectorNorm x
+
+
+----------------------------------------------------------------
+
+-- With an error tolerance of 8 ULPs, a million QuickCheck tests are
+-- likely to all succeed. With a tolerance of 7, we fail around the
+-- half million mark.
+ulps :: Int
+ulps = 8
+
+c_near :: CD -> CD -> Bool
+c_near (a :+ b) (c :+ d) = within ulps a c && within ulps b d
+
+-- Norm of vector
+vectorNorm :: U.Vector CD -> Double
+vectorNorm = sqrt . U.sum . U.map (\(x :+ y) -> x*x + y*y)