statistics 0.10.1.0 → 0.10.2.0

raw patch · 27 files changed

+1117/−66 lines, 27 filesdep ~math-functionsPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: math-functions

API changes (from Hackage documentation)

Files

• README.markdown +0/−15
• Statistics/Distribution.hs +8/−1
• Statistics/Distribution/Beta.hs +98/−0
• Statistics/Distribution/Binomial.hs +3/−3
• Statistics/Distribution/CauchyLorentz.hs +3/−0
• Statistics/Distribution/ChiSquared.hs +6/−1
• Statistics/Distribution/Exponential.hs +6/−2
• Statistics/Distribution/FDistribution.hs +4/−1
• Statistics/Distribution/Gamma.hs +15/−3
• Statistics/Distribution/Geometric.hs +2/−3
• Statistics/Distribution/Normal.hs +2/−2
• Statistics/Distribution/Poisson.hs +3/−3
• Statistics/Distribution/StudentT.hs +3/−0
• Statistics/Distribution/Uniform.hs +14/−7
• Statistics/Sample/Histogram.hs +2/−1
• Statistics/Sample/KernelDensity.hs +11/−10
• Statistics/Test/WilcoxonT.hs +5/−2
• Statistics/Transform.hs +22/−9
• statistics.cabal +26/−3
• tests/Tests/Distribution.hs +284/−0
• tests/Tests/Function.hs +26/−0
• tests/Tests/Helpers.hs +100/−0
• tests/Tests/KDE.hs +45/−0
• tests/Tests/NonparametricTest.hs +233/−0
• tests/Tests/NonparametricTest/Table.hs +36/−0
• tests/Tests/Transform.hs +145/−0
• tests/tests.hs +15/−0

@@ -14,21 +14,6 @@ libraries and applications that use this library using a high level of optimisation. -Suggested GHC options:--    -O -funbox-strict-fields--To illustrate, here are the times (in seconds) to generate and sum 250-million random Word32 values, on a laptop with a 2.4GHz Core2 Duo-P8600 processor, running Fedora 11 and GHC 6.10.3:--    no flags   200+-    -O           1.249-    -O -fvia-C   0.991--As the numbers above suggest, compiling without optimisation will-yield unacceptable performance.-  # Get involved! 

@@ -24,6 +24,7 @@       -- ** Random number generation     , ContGen(..)     , DiscreteGen(..)+    , genContinous       -- * Helper functions     , findRoot     , sumProbabilities@@ -123,7 +124,13 @@ class (DiscreteDistr d, ContGen d) => DiscreteGen d where   genDiscreteVar :: PrimMonad m => d -> Gen (PrimState m) -> m Int -+-- | Generate variates from continous distribution using inverse+--   transform rule.+genContinous :: (ContDistr d, PrimMonad m) => d -> Gen (PrimState m) -> m Double+genContinous d gen = do+  x <- uniform gen+  return $! quantile d x+{-# INLINE genContinous #-}  data P = P {-# UNPACK #-} !Double {-# UNPACK #-} !Double 

@@ -0,0 +1,98 @@+{-# LANGUAGE DeriveDataTypeable #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Statistics.Distribution.Beta+-- Copyright   :  (C) 2012 Edward Kmett,+-- License     :  BSD-style (see the file LICENSE)+--+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>+-- Stability   :  provisional+-- Portability :  DeriveDataTypeable+--+----------------------------------------------------------------------------+module Statistics.Distribution.Beta+  ( BetaDistribution+    -- * Constructor+  , betaDistr+  , improperBetaDistr+    -- * Accessors+  , bdAlpha+  , bdBeta+  ) where++import Numeric.SpecFunctions           (incompleteBeta, invIncompleteBeta, logBeta)+import Numeric.MathFunctions.Constants (m_NaN)+import qualified Statistics.Distribution as D+import Data.Typeable++-- | The beta distribution+data BetaDistribution = BD+ { bdAlpha :: {-# UNPACK #-} !Double+   -- ^ Alpha shape parameter+ , bdBeta  :: {-# UNPACK #-} !Double+   -- ^ Beta shape parameter+ } deriving (Eq,Read,Show,Typeable)++-- | Create beta distribution. Both shape parameters must be positive.+betaDistr :: Double             -- ^ Shape parameter alpha+          -> Double             -- ^ Shape parameter beta+          -> BetaDistribution+betaDistr a b+  | a > 0 && b > 0 = improperBetaDistr a b+  | otherwise      =+      error $  "Statistics.Distribution.Beta.betaDistr: "+            ++ "shape parameters must be positive. Got a = "+            ++ show a+            ++ " b = "+            ++ show b+{-# INLINE betaDistr #-}++-- | Create beta distribution. This construtor doesn't check parameters.+improperBetaDistr :: Double             -- ^ Shape parameter alpha+                  -> Double             -- ^ Shape parameter beta+                  -> BetaDistribution+improperBetaDistr = BD+{-# INLINE improperBetaDistr #-}++instance D.Distribution BetaDistribution where+  cumulative (BD a b) x+    | x <= 0    = 0+    | x >= 1    = 1+    | otherwise = incompleteBeta a b x+  {-# INLINE cumulative #-}++instance D.Mean BetaDistribution where+  mean (BD a b) = a / (a + b)+  {-# INLINE mean #-}++instance D.MaybeMean BetaDistribution where+  maybeMean = Just . D.mean+  {-# INLINE maybeMean #-}++instance D.Variance BetaDistribution where+  variance (BD a b) = a*b / (apb*apb*(apb+1))+    where apb = a + b+  {-# INLINE variance #-}++instance D.MaybeVariance BetaDistribution where+  maybeVariance = Just . D.variance+  {-# INLINE maybeVariance #-}++instance D.ContDistr BetaDistribution where+  density (BD a b) x+   | a <= 0 || b <= 0 = m_NaN+   | x <= 0 = 0+   | x >= 1 = 0+   | otherwise = exp $ (a-1)*log x + (b-1)*log (1-x) - logBeta a b+  {-# INLINE density #-}++  quantile (BD a b) p+    | p == 0         = 0+    | p == 1         = 1+    | p > 0 && p < 1 = invIncompleteBeta a b p+    | otherwise      =+        error $ "Statistics.Distribution.Gamma.quantile: p must be in [0,1] range. Got: "++show p+  {-# INLINE quantile #-}++instance D.ContGen BetaDistribution where+  genContVar = D.genContinous

@@ -85,13 +85,13 @@ {-# INLINE variance #-}  -- | Construct binomial distribution. Number of trials must be---   positive and probability must be in [0,1] range+--   non-negative and probability must be in [0,1] range binomial :: Int                 -- ^ Number of trials.          -> Double              -- ^ Probability.          -> BinomialDistribution binomial n p -  | n <= 0         = -    error $ msg ++ "number of trials must be positive. Got " ++ show n+  | n < 0          =+    error $ msg ++ "number of trials must be non-negative. Got " ++ show n   | p < 0 || p > 1 =      error $ msg++"probability must be in [0,1] range. Got " ++ show p   | otherwise      = BD n p

@@ -63,3 +63,6 @@     | p == 1         =  1 / 0     | otherwise      =       error $ "Statistics.Distribution.CauchyLorentz..quantile: p must be in [0,1] range. Got: "++show p++instance D.ContGen CauchyDistribution where+  genContVar = D.genContinous

@@ -21,7 +21,8 @@ import Data.Typeable         (Typeable) import Numeric.SpecFunctions (incompleteGamma,invIncompleteGamma,logGamma) -import qualified Statistics.Distribution as D+import qualified Statistics.Distribution         as D+import qualified System.Random.MWC.Distributions as MWC   -- | Chi-squared distribution@@ -63,6 +64,10 @@ instance D.MaybeVariance ChiSquared where     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance++instance D.ContGen ChiSquared where+    genContVar (ChiSquared n) = MWC.chiSquare n+  cumulative :: ChiSquared -> Double -> Double cumulative chi x

@@ -24,8 +24,9 @@     ) where  import Data.Typeable (Typeable)-import qualified Statistics.Distribution as D-import qualified Statistics.Sample as S+import qualified Statistics.Distribution         as D+import qualified Statistics.Sample               as S+import qualified System.Random.MWC.Distributions as MWC import Statistics.Types (Sample)  newtype ExponentialDistribution = ED {@@ -54,6 +55,9 @@ instance D.MaybeVariance ExponentialDistribution where     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance++instance D.ContGen ExponentialDistribution where+  genContVar = MWC.exponential . edLambda  cumulative :: ExponentialDistribution -> Double -> Double cumulative (ED l) x | x <= 0    = 0

@@ -22,7 +22,7 @@   --- | Student-T distribution+-- | F distribution data FDistribution = F { fDistributionNDF1 :: {-# UNPACK #-} !Double                        , fDistributionNDF2 :: {-# UNPACK #-} !Double                        , _pdfFactor        :: {-# UNPACK #-} !Double@@ -74,6 +74,9 @@   maybeStdDev (F n m _)      | m > 4     = Just $ 2 * sqr m * (m + n - 2) / (n * sqr (m - 2) * (m - 4))     | otherwise = Nothing++instance D.ContGen FDistribution where+  genContVar = D.genContinous  sqr :: Double -> Double sqr x = x * x

@@ -19,6 +19,7 @@       GammaDistribution     -- * Constructors     , gammaDistr+    , improperGammaDistr     -- * Accessors     , gdShape     , gdScale@@ -27,8 +28,9 @@ import Data.Typeable (Typeable) import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN) import Numeric.SpecFunctions           (incompleteGamma, invIncompleteGamma)-import Statistics.Distribution.Poisson.Internal as Poisson-import qualified Statistics.Distribution        as D+import Statistics.Distribution.Poisson.Internal  as Poisson+import qualified Statistics.Distribution         as D+import qualified System.Random.MWC.Distributions as MWC  -- | The gamma distribution. data GammaDistribution = GD {@@ -44,10 +46,18 @@ gammaDistr k theta   | k     <= 0 = error $ msg ++ "shape must be positive. Got " ++ show k   | theta <= 0 = error $ msg ++ "scale must be positive. Got " ++ show theta-  | otherwise  = GD k theta+  | otherwise  = improperGammaDistr k theta     where msg = "Statistics.Distribution.Gamma.gammaDistr: " {-# INLINE gammaDistr #-} +-- | Create gamma distribution. This constructor do not check whether+--   parameters are valid+improperGammaDistr :: Double            -- ^ Shape parameter. /k/+                   -> Double            -- ^ Scale parameter, &#977;.+                   -> GammaDistribution+improperGammaDistr = GD+{-# INLINE improperGammaDistr #-}+ instance D.Distribution GammaDistribution where     cumulative = cumulative @@ -70,6 +80,8 @@     maybeStdDev   = Just . D.stdDev     maybeVariance = Just . D.variance +instance D.ContGen GammaDistribution where+    genContVar (GD a l) = MWC.gamma a l   density :: GammaDistribution -> Double -> Double

@@ -59,10 +59,9 @@ geometric :: Double                -- ^ Success rate           -> GeometricDistribution geometric x-  | x < 0 || x > 1 = +  | x >= 0 && x <= 1 = GD x+  | otherwise        =     error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x-  | otherwise      =-    GD x {-# INLINE geometric #-}  probability :: GeometricDistribution -> Int -> Double

@@ -57,8 +57,8 @@     stdDev = stdDev  instance D.ContGen NormalDistribution where-    genContVar d gen = do x <- MWC.standard gen-                          return $! stdDev d * (x - mean d)+    genContVar d = MWC.normal (mean d) (stdDev d)+    {-# INLINE genContVar #-}  -- | Standard normal distribution with mean equal to 0 and variance equal to 1 standard :: NormalDistribution

@@ -63,9 +63,9 @@ -- | Create Poisson distribution. poisson :: Double -> PoissonDistribution poisson l-  | l <= 0    = error $ "Statistics.Distribution.Poisson.poisson:\-                        \ lambda must be positive. Got " ++ show l-  | otherwise = PD l+  | l >=  0   = PD l+  | otherwise = error $ "Statistics.Distribution.Poisson.poisson:\+                        \ lambda must be non-negative. Got " ++ show l {-# INLINE poisson #-}  -- $references

@@ -69,3 +69,6 @@ instance D.MaybeVariance StudentT where   maybeStdDev (StudentT ndf) | ndf > 2   = Just $ ndf / (ndf - 2)                              | otherwise = Nothing++instance D.ContGen StudentT where+  genContVar = D.genContinous

@@ -9,19 +9,26 @@ -- Portability : portable -- -- Variate distributed uniformly in the interval.-module Statistics.Distribution.Uniform (-    UniformDistribution-  , uniformDistr-  ) where+module Statistics.Distribution.Uniform+    (+      UniformDistribution+    -- * Constructors+    , uniformDistr+    -- ** Accessors+    , uniformA+    , uniformB+    ) where  import Data.Typeable (Typeable) import qualified Statistics.Distribution as D import qualified System.Random.MWC       as MWC  --- | Uniform distribution-data UniformDistribution = UniformDistribution {-# UNPACK #-} !Double {-# UNPACK #-} !Double-                           deriving (Eq,Show,Read,Typeable)+-- | Uniform distribution from A to B+data UniformDistribution = UniformDistribution {+      uniformA :: {-# UNPACK #-} !Double -- ^ Low boundary of distribution+    , uniformB :: {-# UNPACK #-} !Double -- ^ Upper boundary of distribution+    } deriving (Eq, Read, Show, Typeable)  -- | Create uniform distribution. uniformDistr :: Double -> Double -> UniformDistribution

@@ -19,6 +19,7 @@     , range     ) where +import Numeric.MathFunctions.Constants (m_epsilon) import Statistics.Function (minMax) import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as GM@@ -72,7 +73,7 @@          GM.write bins b . (+1) =<< GM.read bins b          go (i+1)        len = G.length xs-       d = (hi - lo) / fromIntegral numBins+       d = ((hi - lo) * (1 + realToFrac m_epsilon)) / fromIntegral numBins {-# INLINE histogram_ #-}  -- | /O(n)/ Compute decent defaults for the lower and upper bounds of

@@ -52,7 +52,8 @@ kde n0 xs = kde_ n0 (lo - range / 10) (hi + range / 10) xs   where     (lo,hi) = minMax xs-    range = hi - lo+    range   | U.length xs <= 1 = 1       -- Unreasonable guess+            | otherwise        = hi - lo  -- | Gaussian kernel density estimator for one-dimensional data, using -- the method of Botev et al.@@ -72,19 +73,20 @@      -- ^ Upper bound (@max@) of the mesh range.      -> U.Vector Double -> (U.Vector Double, U.Vector Double) kde_ n0 min max xs+  | U.null xs = error "Statistics.KernelDensity.kde: empty sample"   | n0 < 1    = error "Statistics.KernelDensity.kde: invalid number of points"   | otherwise = (mesh, density)   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 (/(2 * 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))-    !n = fromIntegral ni+    !n  = fromIntegral ni     !ni = nextHighestPowerOfTwo n0-    !r = max - min-    a = dct . G.map (/ G.sum h) $ h-      where h = G.map (/ len) $ histogram_ ni min max xs-    !len = fromIntegral (G.length xs)+    !r  = max - min+    a   = dct . G.map (/ G.sum h) $ h+        where h = G.map (/ len) $ 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 ->               x - (len * (2 * sqrt pi) * go 6 (f 7 x)) ** (-0.4)       where@@ -94,10 +96,9 @@                 iv = G.map sqr $ G.enumFromTo 1 (n-1)         go s !h | s == 1    = h                 | otherwise = go (s-1) (f s time)-          where time = (2 * const * k0 / len / h) ** (2 / (3 + 2 * s))+          where time  = (2 * const * k0 / len / h) ** (2 / (3 + 2 * s))                 const = (1 + 0.5 ** (s+0.5)) / 3-                k0 = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi-    _bandwidth = sqrt t_star * r+                k0    = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi     sqr x = x * x  -- $references

@@ -74,7 +74,7 @@ -- all the coefficients by r down the list. -- -- This list will be processed lazily from the head.-coefficients :: Int -> [Int]+coefficients :: Int -> [Integer] coefficients 1 = [1, 1] -- 1 + x coefficients r = let coeffs = coefficients (r-1)                      (firstR, rest) = splitAt r coeffs@@ -86,7 +86,10 @@  -- This list will be processed lazily from the head. summedCoefficients :: Int -> [Double]-summedCoefficients = map fromIntegral . scanl1 (+) . coefficients+summedCoefficients n+  | n < 1     = error "Statistics.Test.WilcoxonT.summedCoefficients: nonpositive sample size"+  | n > 1023  = error "Statistics.Test.WilcoxonT.summedCoefficients: sample is too large (see bug #18)"+  | otherwise = map fromIntegral $ scanl1 (+) $ coefficients n  -- | Tests whether a given result from a Wilcoxon signed-rank matched-pairs test -- is significant at the given level.

@@ -43,37 +43,50 @@  -- | Discrete cosine transform (DCT-II). dct :: U.Vector Double -> U.Vector Double-dct = dct_ . G.map (:+0)+dct = dctWorker . G.map (:+0) --- | Discrete cosine transform, with complex coefficients (DCT-II).+-- | Discrete cosine transform (DCT-II). Only real part of vector is+--   transformed, imaginary part is ignored. dct_ :: U.Vector CD -> U.Vector Double-dct_ xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)+dct_ = dctWorker . G.map (\(i :+ _) -> i :+ 0)++dctWorker :: U.Vector CD -> U.Vector Double+dctWorker 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-    weights = G.cons 1 . G.generate (len-1) $ \x ->+    weights = G.cons 2 . G.generate (len-1) $ \x ->               2 * exp ((0:+(-1))*fi (x+1)*pi/(2*n))       where n = fi len     len = G.length xs ++ -- | Inverse discrete cosine transform (DCT-III). It's inverse of -- 'dct' only up to scale parameter: -- -- > (idct . dct) x = (* lenngth x) idct :: U.Vector Double -> U.Vector Double-idct = idct_ . G.map (:+0)+idct = idctWorker . G.map (:+0) --- | Inverse discrete cosine transform, with complex coefficients--- (DCT-III).+-- | Inverse discrete cosine transform (DCT-III). Only real part of vector is+--   transformed, imaginary part is ignored. idct_ :: U.Vector CD -> U.Vector Double-idct_ xs = G.generate len interleave+idct_ = idctWorker . G.map (\(i :+ _) -> i :+ 0)++idctWorker :: U.Vector CD -> U.Vector Double+idctWorker xs = G.generate len interleave   where     interleave z | even z    = vals `G.unsafeIndex` halve z                  | otherwise = vals `G.unsafeIndex` (len - halve z - 1)     vals = G.map realPart . ifft $ G.zipWith (*) weights xs-    weights = G.generate len $ \x -> n * exp ((0:+1)*fi x*pi/(2*n))+    weights +      = G.cons n+      $ G.generate (len - 1) $ \x -> 2 * n * exp ((0:+1) * fi (x+1) * pi/(2*n))       where n = fi len     len = G.length xs+  -- | Inverse fast Fourier transform. ifft :: U.Vector CD -> U.Vector CD

@@ -1,5 +1,5 @@ name:           statistics-version:        0.10.1.0+version:        0.10.2.0 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful@@ -22,6 +22,26 @@   * Common statistical tests for significant differences between     samples.   .+  Changes in 0.10.2.0+  .+  * Bugs in DCT and IDCT are fixed.+  .+  * Accesors for uniform distribution are added.+  .+  * 'ContGen' instances for all continous distribtuions are added.+  .+  * Beta distribution is added.+  .+  * Constructor for improper gamma distribtuion is added.+  .+  * Binomial distribution allows zero trials.+  .+  * Poisson distribution now accept zero parameter.+  .+  * Integer overflow in caculation of Wilcoxon-T test is fixed.+  .+  * Bug in 'ContGen' instance for normal distribution is fixed.+  .   Changes in 0.10.1.0   .   * Kolmogorov-Smirnov nonparametric test added.@@ -135,6 +155,7 @@     Statistics.Autocorrelation     Statistics.Constants     Statistics.Distribution+    Statistics.Distribution.Beta     Statistics.Distribution.Binomial     Statistics.Distribution.CauchyLorentz     Statistics.Distribution.ChiSquared@@ -199,10 +220,11 @@   other-modules:     Tests.Distribution     Tests.Helpers-    Tests.Math-    Tests.Math.Tables+    Tests.Function     Tests.NonparametricTest+    Tests.NonparametricTest.Table     Tests.Transform+    Tests.KDE    ghc-options:     -Wall -threaded -rtsopts@@ -215,6 +237,7 @@     test-framework,     test-framework-quickcheck2,     test-framework-hunit,+    math-functions,     statistics,     primitive,     vector,

@@ -0,0 +1,284 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# LANGUAGE ScopedTypeVariables  #-}+-- Required for Param+{-# LANGUAGE FlexibleInstances    #-}+{-# LANGUAGE OverlappingInstances #-}+{-# LANGUAGE ViewPatterns #-}+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.Beta+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 BetaDistribution        )+  , 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 (snd . properFraction -> 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+  = printTestCase (printf "CDF = %g" p1)+  $ printTestCase (printf "Sum = %g" p2)+  $ printTestCase (printf "Δ   = %g" (abs (p1 - p2)))+  $ 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  = n + (abs (a - b) `mod` 100)+    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 BetaDistribution where+  arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,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

@@ -0,0 +1,26 @@+module Tests.Function ( tests ) where++import qualified Data.Vector.Unboxed as U+import           Data.Vector.Unboxed   ((!))++import Test.QuickCheck+import Test.Framework+import Test.Framework.Providers.QuickCheck2++import Statistics.Function++++tests :: Test+tests = testGroup "S.Function"+  [ testProperty "Sort is sort" p_sort+  ]+++p_sort :: [Double] -> Property+p_sort xs =+  not (null xs) ==> U.all (uncurry (<=)) (U.zip v $ U.tail v)+    where+      v = sort $ U.fromList xs++      

@@ -0,0 +1,100 @@+-- | Helpers for testing+module Tests.Helpers (+    -- * helpers+    T(..)+  , typeName+  , eq+  , eqC+  , (=~)+    -- * Generic QC tests+  , monotonicallyIncreases+  , monotonicallyIncreasesIEEE+    -- * HUnit helpers+  , testAssertion+  , testEquality+  ) 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 Numeric.MathFunctions.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++testEquality :: (Show a, Eq a) => String -> a -> a -> Test+testEquality msg a b = testCase msg $ HU.assertEqual msg a b

@@ -0,0 +1,45 @@+-- | Tests for Kernel density estimates.+module Tests.KDE ( +  tests +  )where++import Data.Vector.Unboxed ((!))+import qualified Data.Vector.Unboxed as U++import Test.Framework                       (Test, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Test.QuickCheck                      +import Text.Printf++import Statistics.Sample.KernelDensity++++tests :: Test+tests = testGroup "KDE"+  [ testProperty "integral(PDF) == 1" t_densityIsPDF+  ]++t_densityIsPDF :: [Double] -> Property+t_densityIsPDF vec +  = not (null vec) ==> test+  where+    (xs,ys)  = kde 4096 (U.fromList vec)+    step     = (xs ! 1) - (xs ! 0)+    integral = integratePDF step ys+    --+    test = printTestCase (printf "Integral %f" integral)+         $ abs (1 - integral) <= 1e-3++          ++integratePDF :: Double -> U.Vector Double -> Double+integratePDF step vec +  = step * U.sum (U.zipWith (*) vec weights)+  where+    n       = U.length vec+    weights = U.generate n go+      where+        go i | i == 0    = 0.5+             | i == n-1  = 0.5+             | otherwise = 1

@@ -0,0 +1,233 @@+-- Tests for Statistics.Test.NonParametric+module Tests.NonparametricTest (+  nonparametricTests+  ) where+++import qualified Data.Vector.Unboxed as U+import Test.HUnit                     (assertEqual)+import Test.Framework+import Test.Framework.Providers.HUnit++import Statistics.Test.MannWhitneyU+import Statistics.Test.WilcoxonT++import Tests.Helpers+import Tests.NonparametricTest.Table++import Statistics.Test.KolmogorovSmirnov+import Statistics.Distribution.Normal    (standard)++++nonparametricTests :: Test+nonparametricTests = testGroup "Nonparametric tests"+                   $ concat [ mannWhitneyTests+                            , wilcoxonSumTests+                            , wilcoxonPairTests+                            , kolmogorovSmirnovDTest+                            ]++----------------------------------------------------------------++mannWhitneyTests :: [Test]+mannWhitneyTests = zipWith test [(0::Int)..] testData +++  [ testEquality "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]]+  , testEquality "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]]+  , testEquality "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]]+  , testEquality "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 "Mann-Whitney" $ do+          assertEqual ("Mann-Whitney U "     ++ show n) c us+          assertEqual ("Mann-Whitney U Sig " ++ show n) d ss+      where+        us = mannWhitneyU (U.fromList a) (U.fromList b)+        ss = mannWhitneyUSignificant TwoTailed (length a, length b) 0.05 us+    -- 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 "Wilcoxon Sum"+                     $ 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:+  [ testAssertion "Sig 16, 35" (to4dp 0.0467 $ wilcoxonMatchedPairSignificance 16 35)+  , testAssertion "Sig 16, 36" (to4dp 0.0523 $ wilcoxonMatchedPairSignificance 16 36)+  , testEquality   "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]]+  , testEquality "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]]+  , testEquality "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]]+  , testEquality "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) = testEquality ("Wilcoxon Paired " ++ show n) c res+      where res = (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++++----------------------------------------------------------------+-- K-S test+----------------------------------------------------------------+++kolmogorovSmirnovDTest :: [Test]+kolmogorovSmirnovDTest =+  [ testAssertion "K-S D statistics" $+    and [ eq 1e-6 (kolmogorovSmirnovD standard (toU sample)) reference+        | (reference,sample) <- tableKSD+        ]+  , testAssertion "K-S 2-sample statistics" $+    and [ eq 1e-6 (kolmogorovSmirnov2D (toU xs) (toU ys)) reference+        | (reference,xs,ys) <- tableKS2D+        ]+  , testAssertion "K-S probability" $+    and [ eq 1e-14 (kolmogorovSmirnovProbability n d) p+        | (d,n,p) <- testData+        ]+  ]+  where+    toU = U.fromList+    -- Test data for the calculation of cumulative probability +    -- P(D[n] < d).+    -- +    -- Test data is:+    --    (D[n], n, p)+    -- Table is generated using sample program from paper+    testData :: [(Double,Int,Double)]+    testData = +      [ (0.09           ,    3, 0                   )+      , (0.2            ,    3, 0.00177777777777778 )+      , (0.301          ,    3, 0.116357025777778   )+      , (0.392          ,    3, 0.383127210666667   )+      , (0.5003         ,    3, 0.667366306558667   )+      , (0.604          ,    3, 0.861569877333333   )+      , (0.699          ,    3, 0.945458198         )+      , (0.802          ,    3, 0.984475216         )+      , (0.9            ,    3, 0.998               )+      , (0.09           ,    5, 0                   )+      , (0.2            ,    5, 0.0384              )+      , (0.301          ,    5, 0.33993786080016    )+      , (0.392          ,    5, 0.66931908083712    )+      , (0.5003         ,    5, 0.888397260183794   )+      , (0.604          ,    5, 0.971609957879808   )+      , (0.699          ,    5, 0.994331075994008   )+      , (0.802          ,    5, 0.999391366368064   )+      , (0.9            ,    5, 0.99998             )+      , (0.09           ,    8, 3.37615237575e-06   )+      , (0.2            ,    8, 0.151622071801758   )+      , (0.301          ,    8, 0.613891042670582   )+      , (0.392          ,    8, 0.871491561427005   )+      , (0.5003         ,    8, 0.977534089199071   )+      , (0.604          ,    8, 0.997473116268255   )+      , (0.699          ,    8, 0.999806082005123   )+      , (0.802          ,    8, 0.999995133786947   )+      , (0.9            ,    8, 0.99999998          )+      , (0.09           ,   10, 3.89639433093119e-05)+      , (0.2            ,   10, 0.25128096          )+      , (0.301          ,   10, 0.732913126355935   )+      , (0.392          ,   10, 0.932185254518767   )+      , (0.5003         ,   10, 0.992276179340446   )+      , (0.604          ,   10, 0.999495533516769   )+      , (0.699          ,   10, 0.999979691783985   )+      , (0.802          ,   10, 0.999999801409237   )+      , (0.09           ,   20, 0.00794502217168886 )+      , (0.2            ,   20, 0.647279826376584   )+      , (0.301          ,   20, 0.958017466965765   )+      , (0.392          ,   20, 0.997206424843499   )+      , (0.5003         ,   20, 0.999962641414228   )+      , (0.09           ,   30, 0.0498147538075168  )+      , (0.2            ,   30, 0.842030838984526   )+      , (0.301          ,   30, 0.993403560017612   )+      , (0.392          ,   30, 0.99988478803318    )+      , (0.09           ,  100, 0.629367974413669   )+      ]

@@ -0,0 +1,36 @@+module Tests.NonparametricTest.Table where++-- Table for Kolmogorov-Smirnov statistics for standard normal+-- distribution. Generated using R.+--+-- First element of tuple is D second is sample for which it was+-- calculated. +tableKSD :: [(Double,[Double])]+tableKSD = +  [ (0.2012078,[1.360645,-0.3151904,-1.245443,0.1741977,-0.1421206,-1.798246,1.171594,-1.335844,-5.050093e-2,1.030063,-1.849005,0.6491455,-0.7028004])+  , (0.2569956,[0.3884734,-1.227821,-0.4166262,0.429118,-0.9280124,0.8025867,-0.6703089,-0.2124872,0.1224496,0.1087734,-4.285284e-2,-1.039936,-0.7071956])+  , (0.1960356,[-1.814745,-0.6327167,0.7082493,0.6264716,1.02061,-0.4094635,0.821026,-0.4255047,-0.4820728,-0.2239833,0.648517,1.114283,0.3610216])+  , (0.2095746,[0.187011,0.1805498,0.4448389,0.6065506,0.2308673,0.5292549,-1.489902,-1.455191,0.5449396,-0.1436403,-0.7977073,-0.2693545,0.8260888,-1.474473,-2.158696e-2,-0.1455387])+  , (0.1922603,[0.5772317,-1.255561,1.605823,0.4923361,0.2470848,1.176101,-0.3767689,-0.6896885,0.4509345,-0.5048447,0.9436534,1.025599,0.2998393,-3.415219e-2,1.264315,-1.44433,-1.646449e-2])+  , (0.2173401,[1.812807,-0.8687497,-0.5710508,1.003647,1.142621,0.6546577,-6.083323e-3,1.628574e-2,1.067499,-1.953143,-0.6060077,1.90859,-0.7480553,0.6715162,-0.928759,1.862,1.604621,-0.2171044,-0.1835918])+  , (0.2510541,[-0.4769572,1.062319,0.9952284,1.198086,1.015589,-0.4154523,-0.6711762,1.202902,0.2217098,5.381759e-2,0.6679715,0.2551287,-0.1371492])+  , (0.1996022,[1.158607,-0.7354863,1.526559,-0.7855418,-2.82999,-0.6045106,-0.1830228,0.3306812,-0.819657,-1.223715,0.2536423,-0.4155781,1.447042])+  , (0.2284761,[1.239965,0.8187093,0.5199788,1.172072,0.748259,1.869376e-2,0.1625921,-1.712065,0.7043582,-1.702702,-0.4792806,-0.1023351,0.1187189])+  , (0.2337866,[0.9417261,-0.1024297,-0.7354359,1.099991,0.801984,-0.3745397,-1.749564,1.795771,1.099963,-0.605557,-2.035897,1.893603,-0.3468928,-0.2593938,2.100988,0.9665698,0.8757091,0.7696328,0.8730729,-0.3990352,2.04361,-0.4617864,-0.155021,2.15774,0.2687795,-0.9853512,-0.3264898,1.260026,4.267695,-0.5571145,0.6307067,-0.1691405,-1.730686])+  , (0.3389167,[2.025542,-1.542641,-1.090238,3.99027,9.949129e-2,-0.8974433,-2.508418,6.390346,-2.675515,1.154459,1.688072,2.220727,-0.4743102])+  , (0.4920231,[0.5192906,-3.260813,-1.245185,1.693084,3.561318,4.058924,2.27063,0.9446943,4.794159,-3.423733,0.8240817,0.644059,0.900175,1.932513,1.024586,2.82823,2.072192,-0.353231,-0.4319673,1.505952,1.0199,4.555054,2.364929,5.531467,3.279415,3.19821,2.726925,1.680027,-0.9041334,-0.8246765,-1.343979,8.454955,1.354581])+  , (0.6727408,[-6.705672,-3.193988,-4.612611,-3.207994,-5.070172,-6.141169,-0.397149,-4.093359,-1.204801,-3.986585,-2.724662,0.9868107,-6.295266,-5.95839,-6.35114,-1.679555,-2.635889,-4.050329,1.557428,-2.548465,-0.9073924,-1.502018,-4.535688,-4.158818,-8.833434,-5.944697,-1.569672,-4.70399,-7.832059,-4.093708,-8.393417,-2.085432,-7.06495,-0.4230419,-3.046822,-3.23895,-0.9265873,-9.227822,3.293713,-5.593577,-5.942398,-4.358421,2.660044,-4.301572,-1.258879,0.1499903,3.572833,-3.19844,0.8652432,-0.3025793,-1.576673,-7.666265,-6.751463,-1.398944,-2.690656,-1.429654,7.508364e-2,0.7998344,-3.562074,-1.021431,1.342968,2.110244,-7.561497,-2.372083,-3.649193,-5.7723,-1.068083,0.7537809,-4.569546,-1.198005,-5.638384,-1.227226,-1.195852,-1.118175,-9.130527,0.9675821,-2.497391,0.5988562,-1.965783,-4.25741,-6.547006,-1.459294,-2.380556,-3.977307,-7.809006,-4.276819,-4.028746,-9.055546,-3.599239,-1.470512,-8.253329,-1.351687,-4.269324,-6.140353,-6.30808,-1.834091,-3.135146,-9.391791,3.117815,-5.554733,-2.556769,-3.287376,-2.064013,-5.741995,-5.047918,-4.808841,-1.488526,-0.2351115,-5.760833,-2.722929,-7.012353,2.281171,-3.890514,-1.516824,-1.41011,-1.828457,-5.561244,-3.472142,-10.16919,-0.4369042,-5.698953,-4.587462,-4.897086])+  ]++-- Table for 2-sample Kolmogorov-Smirnov statistics. Generated using R+--+-- First element is D, second and third are samples+tableKS2D :: [(Double,[Double],[Double])]+tableKS2D =+  [ (0.2820513,[-0.4212928,2.146532,0.7585263,-0.5086105,-0.7725486,6.235548e-2,-0.1849861,0.861972,-0.1958534,-3.379697e-2,-1.316854,0.6701269],[0.4957582,0.4241167,0.9822869,0.4504248,-0.1749617,1.178098,-1.117222,-0.859273,0.3073736,0.4344583,-0.4761338,-1.332374,1.487291])+  , (0.2820513,[-0.712252,0.7990333,-0.7968473,1.443609,1.163096,-1.349071,-0.1553941,-2.003104,-0.3400618,-0.7019282,0.183293,-0.2352167],[-0.4622455,-0.8132221,0.1161614,-1.472115e-2,1.001454,-6.557789e-2,-0.2531216,-1.032432,0.4105478,1.749614,0.9722899,5.850942e-2,-0.3352746])+  , (0.2564103,[0.3509882,-0.2982833,1.314731,1.264223,-0.8156374,0.3734029,-3.288915e-2,0.6766016,0.9786335,0.1079949,-0.4211722,1.58656],[0.8024675,7.464538e-2,0.2739861,-2.334255e-2,0.5611802,0.6683374,0.4358206,0.349843,1.207834,1.402578,-0.4049183,0.4286042,1.665129])+  , (0.1833333,[1.376196,9.926384e-2,2.199292,-2.04993,0.5585353,-0.4812132,0.1041527,2.084774,0.71194,-1.398245,-4.458574e-2,1.484945,-1.473182,1.020076,-0.7019646,0.2182066,-1.702963,-0.3522622,-0.8129267,-0.6338972],[-1.020371,0.3323861,1.513288,0.1958708,-1.0723,5.323446e-2,-0.9993713,-0.7046356,-0.6781067,-0.4471603,1.512042,-0.2650665,-4.765228e-2,-1.501205,1.228664,0.5332935,-0.2960315,-0.1509683])+  , (0.5666667,[0.7145305,0.1255674,2.001531,0.1419216],[2.113474,-0.3352839,-0.4962429,-1.386079,0.6404667,-0.7145304,0.1084008,-0.9821421,-2.270472,-1.003846,-0.5644588,2.699695,-1.296494,-0.1538839,1.319094,-1.127544,0.3568889,0.2004726,-1.313291,0.3581084,0.3313498,0.9336278,0.9850203,-1.309506,1.170459,-0.7517466,-1.771269,0.7156381,-1.129691,0.877729])+  , (0.5,[0.6950626,0.1643805,-0.3102472,0.4810762,0.1844602,1.338836,-0.8083386,-0.5482141,0.9532421,-0.2644837],[7.527945,-1.95654,1.513725,-1.318431,2.453895,0.2078194,0.7371092,2.834245,-2.134794,3.938259])+  ] 

@@ -0,0 +1,145 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts  #-}+{-# LANGUAGE ViewPatterns      #-}+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,Arbitrary(..),Gen,choose,vectorOf,+                                             printTestCase, quickCheck)++import Text.Printf++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)+        , testProperty "idct . dct = id [up to scale]"+            (t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))+        , testProperty "dct . idct = id [up to scale]"+            (t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))+          -- Exact small size DCT+          -- 2+        , testDCT [1,0] $ map (*2) [1, cos (pi/4)   ]+        , testDCT [0,1] $ map (*2) [1, cos (3*pi/4) ]+          -- 4+        , testDCT [1,0,0,0] $ map (*2) [1, cos(  pi/8), cos( 2*pi/8), cos( 3*pi/8)]+        , testDCT [0,1,0,0] $ map (*2) [1, cos(3*pi/8), cos( 6*pi/8), cos( 9*pi/8)]+        , testDCT [0,0,1,0] $ map (*2) [1, cos(5*pi/8), cos(10*pi/8), cos(15*pi/8)]+        , testDCT [0,0,0,1] $ map (*2) [1, cos(7*pi/8), cos(14*pi/8), cos(21*pi/8)]+          -- Exact small size IDCT+          -- 2+        , testIDCT [1,0]            [1,         1          ]+        , testIDCT [0,1] $ map (*2) [cos(pi/4), cos(3*pi/4)]+          -- 4+        , testIDCT [1,0,0,0]            [1,            1,            1,            1            ]+        , testIDCT [0,1,0,0] $ map (*2) [cos(   pi/8), cos( 3*pi/8), cos( 5*pi/8), cos( 7*pi/8) ]+        , testIDCT [0,0,1,0] $ map (*2) [cos( 2*pi/8), cos( 6*pi/8), cos(10*pi/8), cos(14*pi/8) ]+        , testIDCT [0,0,0,1] $ map (*2) [cos( 3*pi/8), cos( 9*pi/8), cos(15*pi/8), cos(21*pi/8) ]+        ]++-- 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 :: (HasNorm (U.Vector a), U.Unbox a, Num a, Show a, Arbitrary a)+             => (U.Vector a -> U.Vector a) -> Property+t_fftInverse roundtrip = do+  x <- genFftVector+  let n  = G.length x+      x' = roundtrip x+      d  = G.zipWith (-) x x'+      nd = vectorNorm d+      nx = vectorNorm x+  id $ printTestCase "Original vector"+     $ printTestCase (show x )+     $ printTestCase "Transformed one"+     $ printTestCase (show x')+     $ printTestCase (printf "Length = %i" n)+     $ printTestCase (printf "|x - x'| / |x| = %.6g" (nd / nx))+     $ nd <= 3e-14 * nx++-- Test discrete cosine transform+testDCT :: [Double] -> [Double] -> Test+testDCT (U.fromList -> vec) (U.fromList -> res)+  = testAssertion ("DCT test for " ++ show vec)+  $ vecEqual 3e-14 (dct vec) res++-- Test inverse discrete cosine transform+testIDCT :: [Double] -> [Double] -> Test+testIDCT (U.fromList -> vec) (U.fromList -> res)+  = testAssertion ("IDCT test for " ++ show vec)+  $ vecEqual 3e-14 (idct vec) res++++----------------------------------------------------------------++-- 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++-- Arbitrary vector for FFT od DCT+genFftVector :: (U.Unbox a, Arbitrary a) => Gen (U.Vector a)+genFftVector = do+  n <- (2^)  <$> choose (1,9::Int)    -- Size of vector+  G.fromList <$> vectorOf n arbitrary -- Vector to transform++-- Ad-hoc type class for calculation of vector norm+class HasNorm a where+  vectorNorm :: a -> Double++instance HasNorm (U.Vector Double) where+  vectorNorm = sqrt . U.sum . U.map (\x -> x*x)++instance HasNorm (U.Vector CD) where+  vectorNorm = sqrt . U.sum . U.map (\(x :+ y) -> x*x + y*y)++-- Approximate equality for vectors+vecEqual :: Double -> U.Vector Double -> U.Vector Double -> Bool+vecEqual ε v u+  = vectorNorm (U.zipWith (-) v u) < ε * vectorNorm v

@@ -0,0 +1,15 @@+import Test.Framework       (defaultMain)++import Tests.Distribution+import Tests.NonparametricTest+import qualified Tests.Transform+import qualified Tests.Function+import qualified Tests.KDE++main :: IO ()+main = defaultMain [ distributionTests +                   , nonparametricTests+                   , Tests.Transform.tests+                   , Tests.Function.tests+                   , Tests.KDE.tests+                   ]