packages feed

statistics 0.10.0.0 → 0.10.0.1

raw patch · 11 files changed

+892/−11 lines, 11 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

+ Statistics.Transform: dct_ :: Vector CD -> Vector Double
+ Statistics.Transform: idct_ :: Vector CD -> Vector Double
- Statistics.Transform: dct :: Vector CD -> Vector Double
+ Statistics.Transform: dct :: Vector Double -> Vector Double
- Statistics.Transform: idct :: Vector CD -> Vector Double
+ Statistics.Transform: idct :: Vector Double -> Vector Double

Files

Statistics/Distribution/FDistribution.hs view
@@ -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)
Statistics/Sample/KernelDensity.hs view
@@ -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 ->
Statistics/Transform.hs view
@@ -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)
statistics.cabal view
@@ -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
+ tests/Tests/Distribution.hs view
@@ -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
+ tests/Tests/Helpers.hs view
@@ -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
+ tests/Tests/Math.hs view
@@ -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))
+ tests/Tests/Math/Tables.hs view
@@ -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)+  ]
+ tests/Tests/Math/gen.py view
@@ -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+      ])
+ tests/Tests/NonparametricTest.hs view
@@ -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
+ tests/Tests/Transform.hs view
@@ -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)