diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs --- a/Statistics/Distribution/FDistribution.hs +++ b/Statistics/Distribution/FDistribution.hs @@ -25,7 +25,7 @@ -- | 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 --- a/Statistics/Sample/KernelDensity.hs +++ b/Statistics/Sample/KernelDensity.hs @@ -31,7 +31,7 @@ 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 @@ 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 --- 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 @@ 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 @@ 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 --- 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 @@ 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 @@ 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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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 --- /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)