diff --git a/README.markdown b/README.markdown --- a/README.markdown +++ b/README.markdown @@ -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.) diff --git a/Statistics/Constants.hs b/Statistics/Constants.hs --- a/Statistics/Constants.hs +++ b/Statistics/Constants.hs @@ -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 diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs --- a/Statistics/Distribution.hs +++ b/Statistics/Distribution.hs @@ -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 diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs --- a/Statistics/Distribution/Binomial.hs +++ b/Statistics/Distribution/Binomial.hs @@ -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 { diff --git a/Statistics/Distribution/ChiSquared.hs b/Statistics/Distribution/ChiSquared.hs --- a/Statistics/Distribution/ChiSquared.hs +++ b/Statistics/Distribution/ChiSquared.hs @@ -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 diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs --- a/Statistics/Distribution/FDistribution.hs +++ b/Statistics/Distribution/FDistribution.hs @@ -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) diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs --- a/Statistics/Distribution/Gamma.hs +++ b/Statistics/Distribution/Gamma.hs @@ -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 { diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs --- a/Statistics/Distribution/Hypergeometric.hs +++ b/Statistics/Distribution/Hypergeometric.hs @@ -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 { diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs --- a/Statistics/Distribution/Normal.hs +++ b/Statistics/Distribution/Normal.hs @@ -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 diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs --- a/Statistics/Distribution/Poisson.hs +++ b/Statistics/Distribution/Poisson.hs @@ -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) diff --git a/Statistics/Distribution/Poisson/Internal.hs b/Statistics/Distribution/Poisson/Internal.hs --- a/Statistics/Distribution/Poisson/Internal.hs +++ b/Statistics/Distribution/Poisson/Internal.hs @@ -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. diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs --- a/Statistics/Distribution/StudentT.hs +++ b/Statistics/Distribution/StudentT.hs @@ -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) diff --git a/Statistics/Distribution/Uniform.hs b/Statistics/Distribution/Uniform.hs --- a/Statistics/Distribution/Uniform.hs +++ b/Statistics/Distribution/Uniform.hs @@ -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 diff --git a/Statistics/Math.hs b/Statistics/Math.hs --- a/Statistics/Math.hs +++ b/Statistics/Math.hs @@ -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 ∞ 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 --- γ(/s/,/x/). Normalization means that --- γ(/s/,∞)=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 Γ(/x/). Uses --- Algorithm AS 245 by Macleod. --- --- Gives an accuracy of 10–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 ∞ if the input is outside of the range (0 < /x/ --- ≤ 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, Γ(/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 ∞ if the input is outside of the range (0 < /x/ --- ≤ 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@ ≥ 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–256. --- --- * 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–96. --- --- * Loader, C. (2000) Fast and Accurate Computation of Binomial --- Probabilities. --- --- * 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–402. --- --- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete --- gamma integral. /Applied Statistics/ --- 37(3):466–473. --- --- * 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. --- --- * 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 --- --- --- * 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 --- diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs --- a/Statistics/Quantile.hs +++ b/Statistics/Quantile.hs @@ -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 diff --git a/Statistics/Sample/KernelDensity.hs b/Statistics/Sample/KernelDensity.hs --- a/Statistics/Sample/KernelDensity.hs +++ b/Statistics/Sample/KernelDensity.hs @@ -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) diff --git a/Statistics/Sample/KernelDensity/Simple.hs b/Statistics/Sample/KernelDensity/Simple.hs --- a/Statistics/Sample/KernelDensity/Simple.hs +++ b/Statistics/Sample/KernelDensity/Simple.hs @@ -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 diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs --- a/Statistics/Sample/Powers.hs +++ b/Statistics/Sample/Powers.hs @@ -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 diff --git a/Statistics/Test/ChiSquared.hs b/Statistics/Test/ChiSquared.hs new file mode 100644 --- /dev/null +++ b/Statistics/Test/ChiSquared.hs @@ -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 #-} diff --git a/Statistics/Test/KolmogorovSmirnov.hs b/Statistics/Test/KolmogorovSmirnov.hs new file mode 100644 --- /dev/null +++ b/Statistics/Test/KolmogorovSmirnov.hs @@ -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). diff --git a/Statistics/Test/MannWhitneyU.hs b/Statistics/Test/MannWhitneyU.hs --- a/Statistics/Test/MannWhitneyU.hs +++ b/Statistics/Test/MannWhitneyU.hs @@ -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 diff --git a/Statistics/Test/WilcoxonT.hs b/Statistics/Test/WilcoxonT.hs --- a/Statistics/Test/WilcoxonT.hs +++ b/Statistics/Test/WilcoxonT.hs @@ -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 +-- 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'. diff --git a/Statistics/Transform.hs b/Statistics/Transform.hs --- a/Statistics/Transform.hs +++ b/Statistics/Transform.hs @@ -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) diff --git a/statistics.cabal b/statistics.cabal --- a/statistics.cabal +++ b/statistics.cabal @@ -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: diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs deleted file mode 100644 --- a/tests/Tests/Distribution.hs +++ /dev/null @@ -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 diff --git a/tests/Tests/Helpers.hs b/tests/Tests/Helpers.hs deleted file mode 100644 --- a/tests/Tests/Helpers.hs +++ /dev/null @@ -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 diff --git a/tests/Tests/Math.hs b/tests/Tests/Math.hs deleted file mode 100644 --- a/tests/Tests/Math.hs +++ /dev/null @@ -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)) diff --git a/tests/Tests/Math/Tables.hs b/tests/Tests/Math/Tables.hs deleted file mode 100644 --- a/tests/Tests/Math/Tables.hs +++ /dev/null @@ -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) - ] diff --git a/tests/Tests/NonparametricTest.hs b/tests/Tests/NonparametricTest.hs deleted file mode 100644 --- a/tests/Tests/NonparametricTest.hs +++ /dev/null @@ -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 diff --git a/tests/Tests/Transform.hs b/tests/Tests/Transform.hs deleted file mode 100644 --- a/tests/Tests/Transform.hs +++ /dev/null @@ -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) diff --git a/tests/tests.hs b/tests/tests.hs deleted file mode 100644 --- a/tests/tests.hs +++ /dev/null @@ -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 - ]