statistics 0.10.0.1 → 0.10.1.0

raw patch · 31 files changed

+476/−1575 lines, 31 filesdep +math-functionsdep ~mwc-randomPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: math-functions

Dependency ranges changed: mwc-random

API changes (from Hackage documentation)

Files

• README.markdown +5/−5
• Statistics/Constants.hs +6/−71
• Statistics/Distribution.hs +19/−1
• Statistics/Distribution/Binomial.hs +2/−1
• Statistics/Distribution/ChiSquared.hs +2/−2
• Statistics/Distribution/FDistribution.hs +2/−2
• Statistics/Distribution/Gamma.hs +3/−3
• Statistics/Distribution/Hypergeometric.hs +2/−2
• Statistics/Distribution/Normal.hs +5/−1
• Statistics/Distribution/Poisson.hs +1/−1
• Statistics/Distribution/Poisson/Internal.hs +3/−2
• Statistics/Distribution/StudentT.hs +2/−2
• Statistics/Distribution/Uniform.hs +4/−0
• Statistics/Math.hs +11/−626
• Statistics/Quantile.hs +1/−1
• Statistics/Sample/KernelDensity.hs +9/−10
• Statistics/Sample/KernelDensity/Simple.hs +2/−2
• Statistics/Sample/Powers.hs +6/−6
• Statistics/Test/ChiSquared.hs +39/−0
• Statistics/Test/KolmogorovSmirnov.hs +308/−0
• Statistics/Test/MannWhitneyU.hs +2/−1
• Statistics/Test/WilcoxonT.hs +4/−1
• Statistics/Transform.hs +12/−8
• statistics.cabal +26/−6
• tests/Tests/Distribution.hs +0/−276
• tests/Tests/Helpers.hs +0/−96
• tests/Tests/Math.hs +0/−159
• tests/Tests/Math/Tables.hs +0/−47
• tests/Tests/NonparametricTest.hs +0/−148
• tests/Tests/Transform.hs +0/−82
• tests/tests.hs +0/−13

@@ -33,15 +33,15 @@ # Get involved!  Please report bugs via the-[bitbucket issue tracker](http://bitbucket.org/bos/statistics/issues).+[github issue tracker](https://github.com/bos/statistics/issues). -Master [Mercurial repository](http://bitbucket.org/bos/statistics):+Master [git mirror](https://github.com/bos/statistics): -* `hg clone http://bitbucket.org/bos/statistics`+* `git clone git://github.com/bos/statistics.git` -There's also a [git mirror](http://github.com/bos/statistics):+There's also a [Mercurial mirror](https://bitbucket.org/bos/statistics): -* `git clone git://github.com/bos/statistics.git`+* `hg clone https://bitbucket.org/bos/statistics`  (You can create and contribute changes using either Mercurial or git.) 

@@ -8,78 +8,13 @@ -- Portability : portable -- -- Constant values common to much statistics code.+--+-- DEPRECATED: use module 'Numeric.MathFunctions.Constants' from+-- math-functions.  module Statistics.Constants-    (-      m_epsilon-    , m_huge-    , m_tiny-    , m_1_sqrt_2-    , m_2_sqrt_pi-    , m_ln_sqrt_2_pi-    , m_max_exp-    , m_sqrt_2-    , m_sqrt_2_pi-    , m_pos_inf-    , m_neg_inf-    , m_NaN+{-# DEPRECATED "use module Numeric.MathFunctions.Constants from math-functions" #-}+    ( module Numeric.MathFunctions.Constants     ) where --- | A very large number.-m_huge :: Double-m_huge = 1.7976931348623157e308-{-# INLINE m_huge #-}--m_tiny :: Double-m_tiny = 2.2250738585072014e-308-{-# INLINE m_tiny #-}---- | The largest 'Int' /x/ such that 2**(/x/-1) is approximately--- representable as a 'Double'.-m_max_exp :: Int-m_max_exp = 1024---- | @sqrt 2@-m_sqrt_2 :: Double-m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766-{-# INLINE m_sqrt_2 #-}---- | @sqrt (2 * pi)@-m_sqrt_2_pi :: Double-m_sqrt_2_pi = 2.5066282746310005024157652848110452530069867406099383166299-{-# INLINE m_sqrt_2_pi #-}---- | @2 / sqrt pi@-m_2_sqrt_pi :: Double-m_2_sqrt_pi = 1.1283791670955125738961589031215451716881012586579977136881-{-# INLINE m_2_sqrt_pi #-}---- | @1 / sqrt 2@-m_1_sqrt_2 :: Double-m_1_sqrt_2 = 0.7071067811865475244008443621048490392848359376884740365883-{-# INLINE m_1_sqrt_2 #-}---- | The smallest 'Double' ε such that 1 + ε ≠ 1.-m_epsilon :: Double-m_epsilon = encodeFloat (signif+1) expo - 1.0-    where (signif,expo) = decodeFloat (1.0::Double)---- | @log(sqrt((2*pi))@-m_ln_sqrt_2_pi :: Double-m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171-{-# INLINE m_ln_sqrt_2_pi #-}---- | Positive infinity.-m_pos_inf :: Double-m_pos_inf = 1/0-{-# INLINE m_pos_inf #-}---- | Negative infinity.-m_neg_inf :: Double-m_neg_inf = -1/0-{-# INLINE m_neg_inf #-}---- | Not a number.-m_NaN :: Double-m_NaN = 0/0-{-# INLINE m_NaN #-}+import Numeric.MathFunctions.Constants

@@ -21,13 +21,19 @@     , Mean(..)     , MaybeVariance(..)     , Variance(..)+      -- ** Random number generation+    , ContGen(..)+    , DiscreteGen(..)       -- * Helper functions     , findRoot     , sumProbabilities     ) where -import Control.Applicative ((<$>), Applicative(..))+import Control.Applicative     ((<$>), Applicative(..))+import Control.Monad.Primitive (PrimMonad,PrimState)+ import qualified Data.Vector.Unboxed as U+import System.Random.MWC   @@ -104,6 +110,18 @@     variance d = x * x where x = stdDev d     stdDev   :: d -> Double     stdDev = sqrt . variance+++-- | Generate discrete random variates which have given+--   distribution.+class Distribution d => ContGen d where+  genContVar :: PrimMonad m => d -> Gen (PrimState m) -> m Double++-- | Generate discrete random variates which have given+--   distribution. 'ContGen' is superclass because it's always possible+--   to generate real-valued variates from integer values+class (DiscreteDistr d, ContGen d) => DiscreteGen d where+  genDiscreteVar :: PrimMonad m => d -> Gen (PrimState m) -> m Int   

@@ -25,7 +25,8 @@  import Data.Typeable (Typeable) import qualified Statistics.Distribution as D-import Statistics.Math (choose)+import Numeric.SpecFunctions (choose)+  -- | The binomial distribution. data BinomialDistribution = BD {

@@ -18,8 +18,8 @@         , chiSquaredNDF         ) where -import Data.Typeable        (Typeable)-import Statistics.Math      (incompleteGamma,invIncompleteGamma,logGamma)+import Data.Typeable         (Typeable)+import Numeric.SpecFunctions (incompleteGamma,invIncompleteGamma,logGamma)  import qualified Statistics.Distribution as D 

@@ -17,8 +17,8 @@   ) where  import qualified Statistics.Distribution as D-import Data.Typeable   (Typeable)-import Statistics.Math (logBeta, incompleteBeta, invIncompleteBeta)+import Data.Typeable         (Typeable)+import Numeric.SpecFunctions (logBeta, incompleteBeta, invIncompleteBeta)   

@@ -25,10 +25,10 @@     ) where  import Data.Typeable (Typeable)-import Statistics.Constants (m_pos_inf, m_NaN)+import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN)+import Numeric.SpecFunctions           (incompleteGamma, invIncompleteGamma) import Statistics.Distribution.Poisson.Internal as Poisson-import Statistics.Math (incompleteGamma, invIncompleteGamma)-import qualified Statistics.Distribution as D+import qualified Statistics.Distribution        as D  -- | The gamma distribution. data GammaDistribution = GD {

@@ -27,8 +27,8 @@     , hdK     ) where -import Data.Typeable     (Typeable)-import Statistics.Math   (choose)+import Data.Typeable         (Typeable)+import Numeric.SpecFunctions (choose) import qualified Statistics.Distribution as D  data HypergeometricDistribution = HD {

@@ -22,9 +22,10 @@  import Data.Number.Erf (erfc) import Data.Typeable (Typeable)-import Statistics.Constants (m_sqrt_2, m_sqrt_2_pi)+import Numeric.MathFunctions.Constants (m_sqrt_2, m_sqrt_2_pi) import qualified Statistics.Distribution as D import qualified Statistics.Sample as S+import qualified System.Random.MWC.Distributions as MWC  -- | The normal distribution. data NormalDistribution = ND {@@ -55,6 +56,9 @@ instance D.Variance NormalDistribution where     stdDev = stdDev +instance D.ContGen NormalDistribution where+    genContVar d gen = do x <- MWC.standard gen+                          return $! stdDev d * (x - mean d)  -- | Standard normal distribution with mean equal to 0 and variance equal to 1 standard :: NormalDistribution

@@ -27,7 +27,7 @@ import Data.Typeable (Typeable) import qualified Statistics.Distribution as D import qualified Statistics.Distribution.Poisson.Internal as I-import Statistics.Math (incompleteGamma)+import Numeric.SpecFunctions (incompleteGamma)   

@@ -14,8 +14,9 @@       probability     ) where -import Statistics.Constants (m_sqrt_2_pi, m_tiny)-import Statistics.Math (bd0, logGamma, stirlingError)+import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny)+import Numeric.SpecFunctions           (logGamma, stirlingError)+import Numeric.SpecFunctions.Extra     (bd0)  -- | An unchecked, non-integer-valued version of Loader's saddle point -- algorithm.

@@ -16,8 +16,8 @@   ) where  import qualified Statistics.Distribution as D-import Data.Typeable   (Typeable)-import Statistics.Math (logBeta, incompleteBeta, invIncompleteBeta)+import Data.Typeable         (Typeable)+import Numeric.SpecFunctions (logBeta, incompleteBeta, invIncompleteBeta)   

@@ -16,6 +16,7 @@  import Data.Typeable (Typeable) import qualified Statistics.Distribution as D+import qualified System.Random.MWC       as MWC   -- | Uniform distribution@@ -60,3 +61,6 @@  instance D.MaybeVariance UniformDistribution where     maybeStdDev   = Just . D.stdDev++instance D.ContGen UniformDistribution where+    genContVar (UniformDistribution a b) gen = MWC.uniformR (a,b) gen

@@ -9,634 +9,19 @@ -- Portability : portable -- -- Mathematical functions for statistics.+--+-- DEPRECATED. Use package math-functions instead. This module is just+-- reexports functions from 'Numeric.SpecFunctions',+-- 'Numeric.SpecFunctions.Extra' and 'Numeric.Polynomial.Chebyshev'.  module Statistics.Math-    (-    -- * Functions-      choose-    -- ** Beta function-    , logBeta-    , incompleteBeta-    , incompleteBeta_-    , invIncompleteBeta-    -- ** Chebyshev polynomials-    -- $chebyshev-    , chebyshev-    , chebyshevBroucke-    -- ** Factorial-    , factorial-    , logFactorial-    -- ** Gamma function-    , logGamma-    , logGammaL-    , incompleteGamma-    , invIncompleteGamma-    -- ** Logarithm-    , log1p-    , log2-    -- ** Stirling's approximation-    , stirlingError-    , bd0-    -- * References-    -- $references+{-# DEPRECATED "Use package math-function" #-} +    ( module Numeric.Polynomial.Chebyshev+    , module Numeric.SpecFunctions+    , module Numeric.SpecFunctions.Extra     ) where -import Data.Bits ((.&.), (.|.), shiftR)-import Data.Int (Int64)-import Data.Word (Word64)-import Statistics.Constants (m_epsilon, m_sqrt_2_pi, m_ln_sqrt_2_pi, m_NaN,-                             m_neg_inf, m_pos_inf)-import Statistics.Distribution (cumulative)-import Statistics.Distribution.Normal (standard)-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Generic as G----- $chebyshev------ A Chebyshev polynomial of the first kind is defined by the--- following recurrence:------ > t 0 _ = 1--- > t 1 x = x--- > t n x = 2 * x * t (n-1) x - t (n-2) x--data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double---- | Evaluate a Chebyshev polynomial of the first kind. Uses--- Clenshaw's algorithm.-chebyshev :: (G.Vector v Double) =>-             Double      -- ^ Parameter of each function.-          -> v Double    -- ^ Coefficients of each polynomial term, in increasing order.-          -> Double-chebyshev x a = fini . G.foldr' step (C 0 0) . G.tail $ a-    where step k (C b0 b1) = C (k + x2 * b0 - b1) b0-          fini   (C b0 b1) = G.head a + x * b0 - b1-          x2               = x * 2-{-# INLINE chebyshev #-}--data B = B {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double---- | Evaluate a Chebyshev polynomial of the first kind. Uses Broucke's--- ECHEB algorithm, and his convention for coefficient handling, and so--- gives different results than 'chebyshev' for the same inputs.-chebyshevBroucke :: (G.Vector v Double) =>-             Double      -- ^ Parameter of each function.-          -> v Double    -- ^ Coefficients of each polynomial term, in increasing order.-          -> Double-chebyshevBroucke x = fini . G.foldr' step (B 0 0 0)-    where step k (B b0 b1 _) = B (k + x2 * b0 - b1) b0 b1-          fini   (B b0 _ b2) = (b0 - b2) * 0.5-          x2                 = x * 2-{-# INLINE chebyshevBroucke #-}---- | Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with--- no checking.-logChooseFast :: Double -> Double -> Double-logChooseFast n k = -log (n + 1) - logBeta (n - k + 1) (k + 1)---- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For--- values of /k/ > 30, this uses an approximation for performance--- reasons.  The approximation is accurate to 12 decimal places in the--- worst case------ Example:------ > 7 `choose` 3 == 35-choose :: Int -> Int -> Double-n `choose` k-    | k  > n         = 0-    | k' < 50        = U.foldl' go 1 . U.enumFromTo 1 $ k'-    | approx < max64 = fromIntegral . round64 $ approx-    | otherwise      = approx-  where-    k'             = min k (n-k)-    approx         = exp $ logChooseFast (fromIntegral n) (fromIntegral k')-                  -- Less numerically stable:-                  -- exp $ lg (n+1) - lg (k+1) - lg (n-k+1)-                  --   where lg = logGamma . fromIntegral-    go a i         = a * (nk + j) / j-        where j    = fromIntegral i :: Double-    nk             = fromIntegral (n - k')-    max64          = fromIntegral (maxBound :: Int64)-    round64 x      = round x :: Int64--data F = F {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64---- | Compute the factorial function /n/!.  Returns &#8734; if the--- input is above 170 (above which the result cannot be represented by--- a 64-bit 'Double').-factorial :: Int -> Double-factorial n-    | n < 0     = error "Statistics.Math.factorial: negative input"-    | n <= 1    = 1-    | n <= 14   = fini . U.foldl' goLong (F 1 1) $ ns-    | otherwise = U.foldl' goDouble 1 ns-    where goDouble t k = t * fromIntegral k-          goLong (F z x) _ = F (z * x') x'-              where x' = x + 1-          fini (F z _) = fromIntegral z-          ns = U.enumFromTo 2 n---- | Compute the natural logarithm of the factorial function.  Gives--- 16 decimal digits of precision.-logFactorial :: Int -> Double-logFactorial n-    | n <= 14   = log (factorial n)-    | otherwise = (x - 0.5) * log x - x + 9.1893853320467e-1 + z / x-    where x = fromIntegral (n + 1)-          y = 1 / (x * x)-          z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y --               2.7777777777778e-3) * y + 8.3333333333333e-2---- | Compute the normalized lower incomplete gamma function--- &#947;(/s/,/x/). Normalization means that--- &#947;(/s/,&#8734;)=1. Uses Algorithm AS 239 by Shea.-incompleteGamma :: Double       -- ^ /s/-                -> Double       -- ^ /x/-                -> Double-incompleteGamma p x-    | x < 0 || p <= 0 = m_pos_inf-    | x == 0          = 0-    | p >= 1000       = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))-    | x >= 1e8        = 1-    | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1)-                            g = a + log (pearson p 1 1)-                        in if g > limit then exp g else 0-    | otherwise       = let g = p * log x - x - logGamma p + log cf-                        in if g > limit then 1 - exp g else 1-  where-    norm = cumulative standard-    pearson !a !c !g-        | c' <= tolerance = g'-        | otherwise       = pearson a' c' g'-        where a' = a + 1-              c' = c * x / a'-              g' = g + c'-    cf = let a = 1 - p-             b = a + x + 1-             p3 = x + 1-             p4 = x * b-         in contFrac a b 0 1 x p3 p4 (p3/p4)-    contFrac !a !b !c !p1 !p2 !p3 !p4 !g-        | abs (g - rn) <= min tolerance (tolerance * rn) = g-        | otherwise = contFrac a' b' c' (f p3) (f p4) (f p5) (f p6) rn-        where a' = a + 1-              b' = b + 2-              c' = c + 1-              an = a' * c'-              p5 = b' * p3 - an * p1-              p6 = b' * p4 - an * p2-              rn = p5 / p6-              f n | abs p5 > overflow = n / overflow-                  | otherwise         = n-    limit     = -88-    tolerance = 1e-14-    overflow  = 1e37------ Adapted from Numerical Recipes §6.2.1---- | Inverse incomplete gamma function. It's approximately inverse of---   'incompleteGamma' for the same /s/. So following equality---   approximately holds:------ > invIncompleteGamma s . incompleteGamma s = id------   For @invIncompleteGamma s p@ /s/ must be positive and /p/ must be---   in [0,1] range.-invIncompleteGamma :: Double -> Double -> Double-invIncompleteGamma a p-  | a <= 0         = -      error $ "Statistics.Math.invIncompleteGamma: a must be positive. Got: " ++ show a-  | p < 0 || p > 1 = -      error $ "Statistics.Math.invIncompleteGamma: p must be in [0,1] range. Got: " ++ show p-  | p == 0         = 0-  | p == 1         = 1 / 0-  | otherwise      = loop 0 guess-  where-    -- Solve equation γ(a,x) = p using Halley method-    loop :: Int -> Double -> Double-    loop i x-      | i >= 12   = x-      | otherwise =-         let -           -- Value of γ(a,x) - p-           f    = incompleteGamma a x - p-           -- dγ(a,x)/dx-           f'   | a > 1     = afac * exp( -(x - a1) + a1 * (log x - lna1))-                | otherwise = exp( -x + a1 * log x - gln)-           u    = f / f'-           -- Halley correction to Newton-Rapson step-           corr = u * (a1 / x - 1)-           dx   = u / (1 - 0.5 * min 1.0 corr)-           -- New approximation to x-           x'   | x < dx    = 0.5 * x -- Do not go below 0-                | otherwise = x - dx-         in if abs dx < eps * x'-            then x'-            else loop (i+1) x'-    -- Calculate inital guess for root-    guess-      -- -      | a > 1   =-         let t  = sqrt $ -2 * log(if p < 0.5 then p else 1 - p)-             x1 = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t-             x2 = if p < 0.5 then -x1 else x1-         in max 1e-3 (a * (1 - 1/(9*a) - x2 / (3 * sqrt a)) ** 3)-      -- For a <= 1 use following approximations:-      --   γ(a,1) ≈ 0.253a + 0.12a²-      ---      --   γ(a,x) ≈ γ(a,1)·x^a                               x <  1-      --   γ(a,x) ≈ γ(a,1) + (1 - γ(a,1))(1 - exp(1 - x))    x >= 1-      | otherwise =-         let t = 1 - a * (0.253 + a*0.12)-         in if p < t-            then (p / t) ** (1 / a)-            else 1 - log( 1 - (p-t) / (1-t))-    -- Constants-    a1   = a - 1-    lna1 = log a1-    afac = exp( a1 * (lna1 - 1) - gln )-    gln  = logGamma a-    eps  = 1e-8------ Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html---- | Compute the logarithm of the gamma function &#915;(/x/).  Uses--- Algorithm AS 245 by Macleod.------ Gives an accuracy of 10&#8211;12 significant decimal digits, except--- for small regions around /x/ = 1 and /x/ = 2, where the function--- goes to zero.  For greater accuracy, use 'logGammaL'.------ Returns &#8734; if the input is outside of the range (0 < /x/--- &#8804; 1e305).-logGamma :: Double -> Double-logGamma x-    | x <= 0    = m_pos_inf-    | x < 1.5   = a + c *-                  ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) /-                  ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)-    | x < 4     = (x - 2) *-                  ((((r2_4 * x + r2_3) * x + r2_2) * x + r2_1) * x + r2_0) /-                  ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5)-    | x < 12    = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) /-                  ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5)-    | x > 5.1e5 = k-    | otherwise = k + x1 *-                  ((r4_2 * x2 + r4_1) * x2 + r4_0) /-                  ((x2 + r4_4) * x2 + r4_3)-  where-    (a , b , c)-        | x < 0.5   = (-y , x + 1 , x)-        | otherwise = (0  , x     , x - 1)--    y      = log x-    k      = x * (y-1) - 0.5 * y + alr2pi-    alr2pi = 0.918938533204673--    x1 = 1 / x-    x2 = x1 * x1--    r1_0 =  -2.66685511495;   r1_1 =  -24.4387534237;    r1_2 = -21.9698958928-    r1_3 =  11.1667541262;    r1_4 =    3.13060547623;   r1_5 =   0.607771387771-    r1_6 =  11.9400905721;    r1_7 =   31.4690115749;    r1_8 =  15.2346874070--    r2_0 = -78.3359299449;    r2_1 = -142.046296688;     r2_2 = 137.519416416-    r2_3 =  78.6994924154;    r2_4 =    4.16438922228;   r2_5 =  47.0668766060-    r2_6 = 313.399215894;     r2_7 =  263.505074721;     r2_8 =  43.3400022514--    r3_0 =  -2.12159572323e5; r3_1 =    2.30661510616e5; r3_2 =   2.74647644705e4-    r3_3 =  -4.02621119975e4; r3_4 =   -2.29660729780e3; r3_5 =  -1.16328495004e5-    r3_6 =  -1.46025937511e5; r3_7 =   -2.42357409629e4; r3_8 =  -5.70691009324e2--    r4_0 = 0.279195317918525;  r4_1 = 0.4917317610505968;-    r4_2 = 0.0692910599291889; r4_3 = 3.350343815022304-    r4_4 = 6.012459259764103--data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double---- | Compute the logarithm of the gamma function, &#915;(/x/).  Uses a--- Lanczos approximation.------ This function is slower than 'logGamma', but gives 14 or more--- significant decimal digits of accuracy, except around /x/ = 1 and--- /x/ = 2, where the function goes to zero.------ Returns &#8734; if the input is outside of the range (0 < /x/--- &#8804; 1e305).-logGammaL :: Double -> Double-logGammaL x-    | x <= 0    = m_pos_inf-    | otherwise = fini . U.foldl' go (L 0 (x+7)) $ a-    where fini (L l _) = log (l+a0) + log m_sqrt_2_pi - x65 + (x-0.5) * log x65-          go (L l t) k = L (l + k / t) (t-1)-          x65 = x + 6.5-          a0  = 0.9999999999995183-          a   = U.fromList [ 0.1659470187408462e-06-                           , 0.9934937113930748e-05-                           , -0.1385710331296526-                           , 12.50734324009056-                           , -176.6150291498386-                           , 771.3234287757674-                           , -1259.139216722289-                           , 676.5203681218835-                           ]---- | Compute the log gamma correction factor for @x@ &#8805; 10.  This--- correction factor is suitable for an alternate (but less--- numerically accurate) definition of 'logGamma':------ >lgg x = 0.5 * log(2*pi) + (x-0.5) * log x - x + logGammaCorrection x-logGammaCorrection :: Double -> Double-logGammaCorrection x-    | x < 10    = m_NaN-    | x < big   = chebyshevBroucke (t * t * 2 - 1) coeffs / x-    | otherwise = 1 / (x * 12)-  where-    big    = 94906265.62425156-    t      = 10 / x-    coeffs = U.fromList [-               0.1666389480451863247205729650822e+0,-              -0.1384948176067563840732986059135e-4,-               0.9810825646924729426157171547487e-8,-              -0.1809129475572494194263306266719e-10,-               0.6221098041892605227126015543416e-13,-              -0.3399615005417721944303330599666e-15,-               0.2683181998482698748957538846666e-17-             ]---- | Compute the natural logarithm of the beta function.-logBeta :: Double -> Double -> Double-logBeta a b-    | p < 0     = m_NaN-    | p == 0    = m_pos_inf-    | p >= 10   = log q * (-0.5) + m_ln_sqrt_2_pi + logGammaCorrection p + c +-                  (p - 0.5) * log ppq + q * log1p(-ppq)-    | q >= 10   = logGamma p + c + p - p * log pq + (q - 0.5) * log1p(-ppq)-    | otherwise = logGamma p + logGamma q - logGamma pq-    where-      p   = min a b-      q   = max a b-      ppq = p / pq-      pq  = p + q-      c   = logGammaCorrection q - logGammaCorrection pq---- | Regularized incomplete beta function. Uses algorithm AS63 by---   Majumder abd Bhattachrjee.-incompleteBeta :: Double -- ^ /p/ > 0-               -> Double -- ^ /q/ > 0-               -> Double -- ^ /x/, must lie in [0,1] range-               -> Double-incompleteBeta p q = incompleteBeta_ (logBeta p q) p q---- | Regularized incomplete beta function. Same as 'incompleteBeta'--- but also takes value of lo-incompleteBeta_ :: Double -- ^ logarithm of beta function-                -> Double -- ^ /p/ > 0-                -> Double -- ^ /q/ > 0-                -> Double -- ^ /x/, must lie in [0,1] range-                -> Double-incompleteBeta_ beta p q x-  | p <= 0 || q <= 0 = error "p <= 0 || q <= 0"-  | x <  0 || x >  1 = error "x <  0 || x >  1"-  | x == 0 || x == 1 = x-  | p >= (p+q) * x   = incompleteBetaWorker beta p q x-  | otherwise        = 1 - incompleteBetaWorker beta q p (1 - x)---- Worker for incomplete beta function. It is separate function to--- avoid confusion with parameter during parameter swapping-incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double-incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q) :: Int) 1 1 1-  where-    -- Constants-    eps = 1e-15-    cx  = 1 - x-    -- Loop-    loop psq ns ai term betain-      | done      = betain' * exp( p * log x + (q - 1) * log cx - beta) / p-      | otherwise = loop psq' (ns - 1) (ai + 1) term' betain'-      where-        -- New values-        term'   = term * fact / (p + ai)-        betain' = betain + term'-        fact | ns >  0   = (q - ai) * x/cx-             | ns == 0   = (q - ai) * x-             | otherwise = psq * x-        -- Iterations are complete-        done = db <= eps && db <= eps*betain' where db = abs term'-        psq' = if ns < 0 then psq + 1 else psq---- | Compute inverse of regularized incomplete beta function. Uses--- initial approximation from AS109 and Halley method to solve equation.-invIncompleteBeta :: Double     -- ^ /p/-                  -> Double     -- ^ /q/-                  -> Double     -- ^ /a/-                  -> Double-invIncompleteBeta p q a-  | p <= 0 || q <= 0 = error "p <= 0 || q <= 0"-  | a <  0 || a >  1 = error "bad a"-  | a == 0 || a == 1 = a-  | a > 0.5          = 1 - invIncompleteBetaWorker (logBeta p q) q p (1 - a)-  | otherwise        = invIncompleteBetaWorker (logBeta p q) p q a--invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double-invIncompleteBetaWorker beta p q a = loop (0::Int) guess-  where-    p1 = p - 1-    q1 = q - 1-    -- Solve equation using Halley method-    loop !i !x-      | x == 0 || x == 1             = x-      | i >= 10                      = x-      | abs dx <= 16 * m_epsilon * x = x-      | otherwise                    = loop (i+1) x'-      where-        f   = incompleteBeta_ beta p q x - a-        f'  = exp $ p1 * log x + q1 * log (1 - x) - beta-        u   = f / f'-        dx  = u / (1 - 0.5 * min 1 (u * (p1 / x - q1 / (1 - x))))-        x'  | z < 0     = x / 2-            | z > 1     = (x + 1) / 2-            | otherwise = z-            where z = x - dx-    -- Calculate initial guess-    guess -      | p > 1 && q > 1 = -          let rr = (y*y - 3) / 6-              ss = 1 / (2*p - 1)-              tt = 1 / (2*q - 1)-              hh = 2 / (ss + tt)-              ww = y * sqrt(hh + rr) / hh - (tt - ss) * (rr + 5/6 - 2 / (3 * hh))-          in p / (p + q * exp(2 * ww))-      | t'  <= 0  = 1 - exp( (log((1 - a) * q) + beta) / q )-      | t'' <= 1  = exp( (log(a * p) + beta) / p )-      | otherwise = 1 - 2 / (t'' + 1)-      where-        r   = sqrt ( - log ( a * a ) )-        y   = r - ( 2.30753 + 0.27061 * r )-                   / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )-        t   = 1 / (9 * q)-        t'  = 2 * q * (1 - t + y * sqrt t) ** 3-        t'' = (4*p + 2*q - 2) / t'-        -            ---- | Compute the natural logarithm of 1 + @x@.  This is accurate even--- for values of @x@ near zero, where use of @log(1+x)@ would lose--- precision.-log1p :: Double -> Double-log1p x-    | x == 0               = 0-    | x == -1              = m_neg_inf-    | x < -1               = m_NaN-    | x' < m_epsilon * 0.5 = x-    | (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0)-                           = x * (1 - x * 0.5)-    | x' < 0.375           = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs)-    | otherwise            = log (1 + x)-  where-    x' = abs x-    coeffs = U.fromList [-               0.10378693562743769800686267719098e+1,-              -0.13364301504908918098766041553133e+0,-               0.19408249135520563357926199374750e-1,-              -0.30107551127535777690376537776592e-2,-               0.48694614797154850090456366509137e-3,-              -0.81054881893175356066809943008622e-4,-               0.13778847799559524782938251496059e-4,-              -0.23802210894358970251369992914935e-5,-               0.41640416213865183476391859901989e-6,-              -0.73595828378075994984266837031998e-7,-               0.13117611876241674949152294345011e-7,-              -0.23546709317742425136696092330175e-8,-               0.42522773276034997775638052962567e-9,-              -0.77190894134840796826108107493300e-10,-               0.14075746481359069909215356472191e-10,-              -0.25769072058024680627537078627584e-11,-               0.47342406666294421849154395005938e-12,-              -0.87249012674742641745301263292675e-13,-               0.16124614902740551465739833119115e-13,-              -0.29875652015665773006710792416815e-14,-               0.55480701209082887983041321697279e-15,-              -0.10324619158271569595141333961932e-15-             ]---- | Calculate the error term of the Stirling approximation.  This is--- only defined for non-negative values.------ > stirlingError @n@ = @log(n!) - log(sqrt(2*pi*n)*(n/e)^n)-stirlingError :: Double -> Double-stirlingError n -  | n <= 15.0   = case properFraction (n+n) of-                    (i,0) -> sfe `U.unsafeIndex` i-                    _     -> logGamma (n+1.0) - (n+0.5) * log n + n --                             m_ln_sqrt_2_pi-  | n > 500     = (s0-s1/nn)/n-  | n > 80      = (s0-(s1-s2/nn)/nn)/n-  | n > 35      = (s0-(s1-(s2-s3/nn)/nn)/nn)/n-  | otherwise   = (s0-(s1-(s2-(s3-s4/nn)/nn)/nn)/nn)/n-  where-    nn = n*n-    s0 = 0.083333333333333333333        -- 1/12-    s1 = 0.00277777777777777777778      -- 1/360-    s2 = 0.00079365079365079365079365   -- 1/1260-    s3 = 0.000595238095238095238095238  -- 1/1680-    s4 = 0.0008417508417508417508417508 -- 1/1188-    sfe = U.fromList [ 0.0, -                0.1534264097200273452913848,   0.0810614667953272582196702,-                0.0548141210519176538961390,   0.0413406959554092940938221,-                0.03316287351993628748511048,  0.02767792568499833914878929,-                0.02374616365629749597132920,  0.02079067210376509311152277,-                0.01848845053267318523077934,  0.01664469118982119216319487,-                0.01513497322191737887351255,  0.01387612882307074799874573,-                0.01281046524292022692424986,  0.01189670994589177009505572,-                0.01110455975820691732662991,  0.010411265261972096497478567,-                0.009799416126158803298389475, 0.009255462182712732917728637,-                0.008768700134139385462952823, 0.008330563433362871256469318,-                0.007934114564314020547248100, 0.007573675487951840794972024,-                0.007244554301320383179543912, 0.006942840107209529865664152,-                0.006665247032707682442354394, 0.006408994188004207068439631,-                0.006171712263039457647532867, 0.005951370112758847735624416,-                0.005746216513010115682023589, 0.005554733551962801371038690 ]----- | Evaluate the deviance term @x log(x/np) + np - x@.-bd0 :: Double                   -- ^ @x@-    -> Double                   -- ^ @np@-    -> Double -bd0 x np -  | isInfinite x || isInfinite np || np == 0 = m_NaN-  | abs x_np >= 0.1*(x+np)                   = x * log (x/np) - x_np-  | otherwise                                = loop 1 (ej0*vv) s0-  where -    x_np = x - np-    v    = x_np / (x+np)-    s0   = x_np * v-    ej0  = 2*x*v-    vv   = v*v-    loop j ej s = case s + ej/(2*j+1) of-                    s' | s' == s   -> s'-                       | otherwise -> loop (j+1) (ej*vv) s'---- | /O(log n)/ Compute the logarithm in base 2 of the given value.-log2 :: Int -> Int-log2 v0-    | v0 <= 0   = error "Statistics.Math.log2: invalid input"-    | otherwise = go 5 0 v0-  where-    go !i !r !v | i == -1        = r-                | v .&. b i /= 0 = let si = U.unsafeIndex sv i-                                   in go (i-1) (r .|. si) (v `shiftR` si)-                | otherwise      = go (i-1) r v-    b = U.unsafeIndex bv-    !bv = U.fromList [0x2, 0xc, 0xf0, 0xff00, 0xffff0000, 0xffffffff00000000]-    !sv = U.fromList [1,2,4,8,16,32]+import Numeric.Polynomial.Chebyshev+import Numeric.SpecFunctions+import Numeric.SpecFunctions.Extra --- $references------ * Broucke, R. (1973) Algorithm 446: Ten subroutines for the---   manipulation of Chebyshev series. /Communications of the ACM/---   16(4):254&#8211;256.  <http://doi.acm.org/10.1145/362003.362037>------ * Clenshaw, C.W. (1962) Chebyshev series for mathematical---   functions. /National Physical Laboratory Mathematical Tables 5/,---   Her Majesty's Stationery Office, London.------ * Lanczos, C. (1964) A precision approximation of the gamma---   function.  /SIAM Journal on Numerical Analysis B/---   1:86&#8211;96. <http://www.jstor.org/stable/2949767>------ * Loader, C. (2000) Fast and Accurate Computation of Binomial---   Probabilities. <http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf>------ * Macleod, A.J. (1989) Algorithm AS 245: A robust and reliable---   algorithm for the logarithm of the gamma function.---   /Journal of the Royal Statistical Society, Series C (Applied Statistics)/---   38(2):397&#8211;402. <http://www.jstor.org/stable/2348078>------ * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete---   gamma integral. /Applied Statistics/---   37(3):466&#8211;473. <http://www.jstor.org/stable/2347328>------ * K. L. Majumder, G. P. Bhattacharjee (1973) Algorithm AS 63: The---   Incomplete Beta Integral. /Journal of the Royal Statistical---   Society. Series C (Applied Statistics)/ Vol. 22, No. 3 (1973),---   pp. 409-411. <http://www.jstor.org/pss/2346797>------ * K. L. Majumder, G. P. Bhattacharjee (1973) Algorithm AS 64:---   Inverse of the Incomplete Beta Function Ratio. /Journal of the---   Royal Statistical Society. Series C (Applied Statistics)/---   Vol. 22, No. 3 (1973), pp. 411-414---   <http://www.jstor.org/pss/2346798>------ * G. W. Cran, K. J. Martin and G. E. Thomas (1977) Remark AS R19---   and Algorithm AS 109: A Remark on Algorithms: AS 63: The---   Incomplete Beta Integral AS 64: Inverse of the Incomplete Beta---   Function Ratio. /Journal of the Royal Statistical Society. Series---   C (Applied Statistics)/ Vol. 26, No. 1 (1977), pp. 111-114---   <http://www.jstor.org/pss/2346887>

@@ -39,7 +39,7 @@  import Control.Exception (assert) import Data.Vector.Generic ((!))-import Statistics.Constants (m_epsilon)+import Numeric.MathFunctions.Constants (m_epsilon) import Statistics.Function (partialSort) import qualified Data.Vector.Generic as G 

@@ -25,13 +25,12 @@     -- $references     ) where -import Data.Complex (Complex(..)) import Prelude hiding (const,min,max)-import Statistics.Constants (m_sqrt_2_pi)-import Statistics.Function (minMax, nextHighestPowerOfTwo)-import Statistics.Math.RootFinding (fromRoot, ridders)-import Statistics.Sample.Histogram (histogram_)-import Statistics.Transform (dct_, idct_)+import Numeric.MathFunctions.Constants (m_sqrt_2_pi)+import Statistics.Function             (minMax, nextHighestPowerOfTwo)+import Statistics.Math.RootFinding     (fromRoot, ridders)+import Statistics.Sample.Histogram     (histogram_)+import Statistics.Transform            (dct, idct) import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U @@ -78,13 +77,13 @@   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))-      where f b z = b * exp (sqr z * sqr pi * t_star * (-0.5)) :+ 0+    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))     !n = fromIntegral ni     !ni = nextHighestPowerOfTwo n0     !r = max - min-    a = dct_ . G.map (/ G.sum h) $ h-      where h = G.map (/ (len :+ 0)) $ histogram_ ni min max xs+    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)

@@ -46,9 +46,9 @@     -- $references     ) where -import Statistics.Constants (m_1_sqrt_2, m_2_sqrt_pi)+import Numeric.MathFunctions.Constants (m_1_sqrt_2, m_2_sqrt_pi) import Statistics.Function (minMax)-import Statistics.Sample (stdDev)+import Statistics.Sample   (stdDev) import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Generic as G 

@@ -47,13 +47,13 @@     -- $references     ) where -import Data.Vector.Generic (unsafeFreeze)-import Data.Vector.Unboxed ((!))+import Data.Vector.Generic   (unsafeFreeze)+import Data.Vector.Unboxed   ((!)) import Prelude hiding (sum)-import Statistics.Function (indexed)-import Statistics.Internal (inlinePerformIO)-import Statistics.Math (choose)-import System.IO.Unsafe (unsafePerformIO)+import Statistics.Function   (indexed)+import Statistics.Internal   (inlinePerformIO)+import Numeric.SpecFunctions (choose)+import System.IO.Unsafe      (unsafePerformIO) import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed.Mutable as MU

@@ -0,0 +1,39 @@+{-# LANGUAGE FlexibleContexts #-}+-- | Pearson's chi squared test.+module Statistics.Test.ChiSquared (+    chi2test+    -- * Data types+  , TestType(..)+  , TestResult(..)+  ) where++import qualified Data.Vector.Generic as G++import Statistics.Distribution+import Statistics.Distribution.ChiSquared+import Statistics.Test.Types+++-- | Generic form of Pearson chi squared tests for binned data. Data+--   sample is supplied in form of tuples (observed quantity,+--   expected number of events). Both must be positive.+chi2test :: (G.Vector v (Int,Double), G.Vector v Double)+         => Double              -- ^ p-value+         -> Int                 -- ^ Number of additional degrees of+                                --   freedom. One degree of freedom+                                --   is due to the fact that the are+                                --   N observation in total and+                                --   accounted for automatically.+         -> v (Int,Double)      -- ^ Observation and expectation.+         -> TestResult+chi2test p ndf vec+  | ndf < 0        = error $ "Statistics.Test.ChiSquare.chi2test: negative NDF " ++ show ndf+  | n   < 0        = error $ "Statistics.Test.ChiSquare.chi2test: too short data sample"+  | p > 0 && p < 1 = significant $ complCumulative d chi2 < p+  | otherwise      = error $ "Statistics.Test.ChiSquare.chi2test: bad p-value: " ++ show p+  where+    n     = G.length vec - ndf - 1+    chi2  = G.sum $ G.map (\(o,e) -> sqr (fromIntegral o - e) / e) vec+    d     = chiSquared n+    sqr x = x * x+{-# INLINE chi2test #-}

@@ -0,0 +1,308 @@+-- |+-- Module    : Statistics.Test.KolmogorovSmirnov+-- Copyright : (c) 2011 Aleksey Khudyakov+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Kolmogov-Smirnov tests are non-parametric tests for assesing+-- whether given sample could be described by distribution or whether+-- two samples have the same distribution.+module Statistics.Test.KolmogorovSmirnov (+    -- * Kolmogorov-Smirnov test+    kolmogorovSmirnovTest+  , kolmogorovSmirnovTestCdf+  , kolmogorovSmirnovTest2+    -- * Evaluate statistics+  , kolmogorovSmirnovCdfD+  , kolmogorovSmirnovD+  , kolmogorovSmirnov2D+    -- * Probablities+  , kolmogorovSmirnovProbability+    -- * Data types+  , TestType(..)+  , TestResult(..)+    -- * References+    -- $references+  ) where++import Control.Monad+import Control.Monad.ST  (ST)++import qualified Data.Vector.Unboxed         as U+import qualified Data.Vector.Unboxed.Mutable as M++import Statistics.Distribution        (Distribution(..))+import Statistics.Types               (Sample)+import Statistics.Function            (sort)+import Statistics.Test.Types++import Text.Printf++++----------------------------------------------------------------+-- Test+----------------------------------------------------------------++-- | Check that sample could be described by+--   distribution. 'Significant' means distribution is not compatible+--   with data for given p-value.+--+--   This test uses Marsaglia-Tsang-Wang exact alogorithm for+--   calculation of p-value.+kolmogorovSmirnovTest :: Distribution d+                      => d      -- ^ Distribution+                      -> Double -- ^ p-value+                      -> Sample -- ^ Data sample+                      -> TestResult+kolmogorovSmirnovTest d = kolmogorovSmirnovTestCdf (cumulative d)+{-# INLINE kolmogorovSmirnovTest #-}++-- | Variant of 'kolmogorovSmirnovTest' which uses CFD in form of+--   function.+kolmogorovSmirnovTestCdf :: (Double -> Double) -- ^ CDF of distribution+                         -> Double             -- ^ p-value+                         -> Sample             -- ^ Data sample+                         -> TestResult+kolmogorovSmirnovTestCdf cdf p sample+  | p > 0 && p < 1 = significant $ 1 - prob < p+  | otherwise      = error "Statistics.Test.KolmogorovSmirnov.kolmogorovSmirnovTestCdf:bad p-value"+  where+    d    = kolmogorovSmirnovCdfD cdf sample+    prob = kolmogorovSmirnovProbability (U.length sample) d++-- | Two sample Kolmogorov-Smirnov test. It tests whether two data+--   samples could be described by the same distribution without+--   making any assumptions about it.+--+--   This test uses approxmate formula for computing p-value.+kolmogorovSmirnovTest2 :: Double -- ^ p-value+                       -> Sample -- ^ Sample 1+                       -> Sample -- ^ Sample 2+                       -> TestResult+kolmogorovSmirnovTest2 p xs1 xs2+  | p > 0 && p < 1 = significant $ 1 - prob( d*(en + 0.12 + 0.11/en) ) < p+  | otherwise      = error "Statistics.Test.KolmogorovSmirnov.kolmogorovSmirnovTest2:bad p-value"+  where+    d    = kolmogorovSmirnov2D xs1 xs2+    -- Effective number of data points+    n1   = fromIntegral (U.length xs1)+    n2   = fromIntegral (U.length xs2)+    en   = sqrt $ n1 * n2 / (n1 + n2)+    --+    prob z+      | z <  0    = error "kolmogorovSmirnov2D: internal error"+      | z == 0    = 1+      | z <  1.18 = let y = exp( -1.23370055013616983 / (z*z) )+                    in  2.25675833419102515 * sqrt( -log(y) ) * (y + y**9 + y**25 + y**49)+      | otherwise = let x = exp(-2 * z * z)+                    in  1 - 2*(x - x**4 + x**9)+-- FIXME: Find source for approximation for D++++----------------------------------------------------------------+-- Kolmogorov's statistic+----------------------------------------------------------------++-- | Calculate Kolmogorov's statistic /D/ for given cumulative+--   distribution function (CDF) and data sample. If sample is empty+--   returns 0.+kolmogorovSmirnovCdfD :: (Double -> Double) -- ^ CDF function+                      -> Sample             -- ^ Sample+                      -> Double+kolmogorovSmirnovCdfD cdf sample+  | U.null xs = 0+  | otherwise = U.maximum+              $ U.zipWith3 (\p a b -> abs (p-a) `max` abs (p-b))+                  ps steps (U.tail steps)+  where+    xs = sort sample+    n  = U.length xs+    --+    ps    = U.map cdf xs+    steps = U.map ((/ fromIntegral n) . fromIntegral)+          $ U.generate (n+1) id+++-- | Calculate Kolmogorov's statistic /D/ for given cumulative+--   distribution function (CDF) and data sample. If sample is empty+--   returns 0.+kolmogorovSmirnovD :: (Distribution d)+                   => d         -- ^ Distribution+                   -> Sample    -- ^ Sample+                   -> Double+kolmogorovSmirnovD d = kolmogorovSmirnovCdfD (cumulative d)+{-# INLINE kolmogorovSmirnovD #-}++-- | Calculate Kolmogorov's statistic /D/ for two data samples. If+--   either of samples is empty returns 0.+kolmogorovSmirnov2D :: Sample   -- ^ First sample+                    -> Sample   -- ^ Second sample+                    -> Double+kolmogorovSmirnov2D sample1 sample2+  | U.null sample1 || U.null sample2 = 0+  | otherwise                        = worker 0 0 0+  where+    xs1 = sort sample1+    xs2 = sort sample2+    n1  = U.length xs1+    n2  = U.length xs2+    en1 = fromIntegral n1+    en2 = fromIntegral n2+    -- Find new index+    skip x i xs = go (i+1)+      where go n | n >= U.length xs = n+                 | xs U.! n == x    = go (n+1)+                 | otherwise        = n+    -- Main loop+    worker d i1 i2+      | i1 >= n1 || i2 >= n2 = d+      | otherwise            = worker d' i1' i2'+      where+        d1  = xs1 U.! i1+        d2  = xs2 U.! i2+        i1' | d1 <= d2  = skip d1 i1 xs1+            | otherwise = i1+        i2' | d2 <= d1  = skip d2 i2 xs2+            | otherwise = i2+        d'  = max d (abs $ fromIntegral i1' / en1 - fromIntegral i2' / en2)++++-- | Calculate cumulative probability function for Kolmogorov's+--   distribution with /n/ parameters or probability of getting value+--   smaller than /d/ with n-elements sample.+--+--   It uses algorithm by Marsgalia et. al. and provide at least+--   7-digit accuracy.+kolmogorovSmirnovProbability :: Int    -- ^ Size of the sample+                             -> Double -- ^ D value+                             -> Double+kolmogorovSmirnovProbability n d+  -- Avoid potencially lengthy calculations for large N and D > 0.999+  | s > 7.24 || (s > 3.76 && n > 99) = 1 - 2 * exp( -(2.000071 + 0.331 / sqrt n' + 1.409 / n') * s)+  -- Exact computation+  | otherwise = fini $ matrixPower matrix n+  where+    s  = n' * d * d+    n' = fromIntegral n++    size = 2*k - 1+    k    = floor (n' * d) + 1+    h    = fromIntegral k - n' * d+    -- Calculate initial matrix+    matrix =+      let m = U.create $ do+            mat <- M.new (size*size)+            -- Fill matrix with 0 and 1s+            for 0 size $ \row ->+              for 0 size $ \col -> do+                let val | row + 1 >= col = 1+                        | otherwise      = 0 :: Double+                M.write mat (row * size + col) val+            -- Correct left column/bottom row+            for 0 size $ \i -> do+              let delta = h ^^ (i + 1)+              modify mat (i    * size)         (subtract delta)+              modify mat (size * size - 1 - i) (subtract delta)+            -- Correct corner element if needed+            when (2*h > 1) $ do+              modify mat ((size - 1) * size) (+ ((2*h - 1) ^ size))+            -- Divide diagonals by factorial+            let divide g num+                  | num == size = return ()+                  | otherwise   = do for num size $ \i ->+                                       modify mat (i * (size + 1) - num) (/ g)+                                     divide (g * fromIntegral (num+2)) (num+1)+            divide 2 1+            return mat+      in Matrix size m 0+    -- Last calculation+    fini m@(Matrix _ _ e) = loop 1 (matrixCenter m) e+      where+        loop i ss eQ+          | i  > n       = ss * 10 ^^ eQ+          | ss' < 1e-140 = loop (i+1) (ss' * 1e140) (eQ - 140)+          | otherwise    = loop (i+1)  ss'           eQ+          where ss' = ss * fromIntegral i / fromIntegral n+++----------------------------------------------------------------++-- Maxtrix operations.+--+-- There isn't the matrix package for haskell yet so nessesary minimum+-- is implemented here.++-- Square matrix stored in row-major order+data Matrix = Matrix+              {-# UNPACK #-} !Int -- Size of matrix+              !(U.Vector Double)  -- Matrix data+              {-# UNPACK #-} !Int -- In order to avoid overflows+                                  -- during matrix multiplication large+                                  -- exponent is stored seprately++-- Show instance useful mostly for debugging+instance Show Matrix where+  show (Matrix n vs _) = unlines $ map (unwords . map (printf "%.4f")) $ split $ U.toList vs+    where+      split [] = []+      split xs = row : split rest where (row, rest) = splitAt n xs+++-- Avoid overflow in the matrix+avoidOverflow :: Matrix -> Matrix+avoidOverflow m@(Matrix n xs e)+  | matrixCenter m > 1e140 = Matrix n (U.map (* 1e-140) xs) (e + 140)+  | otherwise              = m++-- Unsafe matrix-matrix multiplication. Matrices must be of the same+-- size. This is not checked.+matrixMultiply :: Matrix -> Matrix -> Matrix+matrixMultiply (Matrix n xs e1) (Matrix _ ys e2) =+  Matrix n (U.generate (n*n) go) (e1 + e2)+  where+    go i = U.sum $ U.zipWith (*) row col+      where+        nCol = i `rem` n+        row  = U.slice (i - nCol) n xs+        col  = U.backpermute ys $ U.enumFromStepN nCol n n++-- Raise matrix to power N. power must be positive it's not checked+matrixPower :: Matrix -> Int -> Matrix+matrixPower mat 1 = mat+matrixPower mat n = avoidOverflow res+  where+    mat2 = matrixPower mat (n `quot` 2)+    pow  = matrixMultiply mat2 mat2+    res | odd n     = matrixMultiply pow mat+        | otherwise = pow++-- Element in the center of matrix (Not corrected for exponent)+matrixCenter :: Matrix -> Double+matrixCenter (Matrix n xs _) = (U.!) xs (k*n + k) where k = n `quot` 2++-- Simple for loop+for :: Monad m => Int -> Int -> (Int -> m ()) -> m ()+for n0 n f = loop n0+  where+    loop i | i == n    = return ()+           | otherwise = f i >> loop (i+1)++-- Modify element in the vector+modify :: U.Unbox a => M.MVector s a -> Int -> (a -> a) -> ST s ()+modify arr i f = do x <- M.read arr i+                    M.write arr i (f x)+{-# INLINE modify #-}++----------------------------------------------------------------++-- $references+--+-- * G. Marsaglia, W. W. Tsang, J. Wang (2003) Evaluating Kolmogorov's+--   distribution, Journal of Statistical Software, American+--   Statistical Association, vol. 8(i18).

@@ -31,9 +31,10 @@ import Data.Ord            (comparing) import qualified Data.Vector.Unboxed as U +import Numeric.SpecFunctions          (choose)+ import Statistics.Distribution        (quantile) import Statistics.Distribution.Normal (standard)-import Statistics.Math                (choose) import Statistics.Types               (Sample) import Statistics.Function            (sortBy) import Statistics.Test.Types

@@ -10,6 +10,9 @@ -- The Wilcoxon matched-pairs signed-rank test is non-parametric test -- which could be used to whether two related samples have different -- means.+--+-- WARNING: current implementation contain critical bugs+-- <https://github.com/bos/statistics/issues/18> module Statistics.Test.WilcoxonT (     -- * Wilcoxon signed-rank matched-pair test     wilcoxonMatchedPairTest@@ -99,7 +102,7 @@ -- order to 'wilcoxonMatchedPairSignedRank', or simply swap the values in the resulting -- pair before passing them to this function. wilcoxonMatchedPairSignificant ::-     TestType            -- ^ Perform one-tailed test (see description above).+     TestType            -- ^ Perform one- or two-tailed test (see description below).   -> Int                 -- ^ The sample size from which the (T+,T-) values were derived.   -> Double              -- ^ The p-value at which to test (e.g. 0.05)   -> (Double, Double)    -- ^ The (T+, T-) values from 'wilcoxonMatchedPairSignedRank'.

@@ -29,15 +29,16 @@     , ifft     ) where -import Control.Monad (when)-import Control.Monad.ST (ST)-import Data.Bits (shiftL, shiftR)-import Data.Complex (Complex(..), conjugate, realPart)-import Statistics.Math (log2)-import qualified Data.Vector.Generic as G+import Control.Monad         (when)+import Control.Monad.ST      (ST)+import Data.Bits             (shiftL, shiftR)+import Data.Complex          (Complex(..), conjugate, realPart)+import Numeric.SpecFunctions (log2)+import qualified Data.Vector.Generic         as G import qualified Data.Vector.Generic.Mutable as M-import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed         as U + type CD = Complex Double  -- | Discrete cosine transform (DCT-II).@@ -55,7 +56,10 @@       where n = fi len     len = G.length xs --- | Inverse discrete cosine transform (DCT-III).+-- | 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) 

@@ -1,5 +1,5 @@ name:           statistics-version:        0.10.0.1+version:        0.10.1.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,23 @@   * Common statistical tests for significant differences between     samples.   .+  Changes in 0.10.1.0+  .+  * Kolmogorov-Smirnov nonparametric test added.+  .+  * Pearson's chi squared test added.+  .+  * Type class for generating random variates for given distribution+    is added.+  .+  * Modules 'Statistics.Math' and 'Statistics.Constants' are moved to+    the @math-functions@ package. They are still available but marked+    as deprecated.+  .+  Changed in 0.10.0.1+  .+  * @dct@ and @idct@ now have type @Vector Double -> Vector Double@+  .   Changes in 0.10.0.0:   .   * The type classes @Mean@ and @Variance@ are split in two. This is@@ -142,9 +159,11 @@     Statistics.Sample.KernelDensity.Simple     Statistics.Sample.Powers     Statistics.Test.NonParametric-    Statistics.Test.Types+    Statistics.Test.ChiSquared+    Statistics.Test.KolmogorovSmirnov     Statistics.Test.MannWhitneyU     Statistics.Test.WilcoxonT+    Statistics.Test.Types     Statistics.Transform     Statistics.Types   other-modules:@@ -156,10 +175,11 @@     base < 5,     deepseq >= 1.1.0.2,     erf,-    monad-par >= 0.1.0.1,-    mwc-random >= 0.8.0.5,-    primitive >= 0.3,-    vector >= 0.7.1,+    monad-par         >= 0.1.0.1,+    mwc-random        >= 0.11.0.0,+    math-functions    >= 0.1.1,+    primitive         >= 0.3,+    vector            >= 0.7.1,     vector-algorithms >= 0.4   if impl(ghc >= 6.10)     build-depends:

@@ -1,276 +0,0 @@-{-# 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

@@ -1,96 +0,0 @@--- | 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

@@ -1,159 +0,0 @@-{-# 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))

@@ -1,47 +0,0 @@-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)-  ]

@@ -1,148 +0,0 @@--- 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

@@ -1,82 +0,0 @@-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)

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