packages feed

statistics 0.16.2.1 → 0.16.3.0

raw patch · 18 files changed

+377/−258 lines, 18 filesdep +doctestdep +tasty-benchdep +tasty-papidep ~basedep ~mwc-randomdep ~parallelsetup-changed

Dependencies added: doctest, tasty-bench, tasty-papi

Dependency ranges changed: base, mwc-random, parallel, tasty, vector

Files

− Setup.lhs 
@@ -1,3 +0,0 @@-#!/usr/bin/env runhaskell-> import Distribution.Simple-> main = defaultMain
Statistics/Correlation.hs view
@@ -6,9 +6,11 @@ module Statistics.Correlation     ( -- * Pearson correlation       pearson+    , pearson2     , pearsonMatByRow       -- * Spearman correlation     , spearman+    , spearman2     , spearmanMatByRow     ) where @@ -25,11 +27,18 @@  -- | Pearson correlation for sample of pairs. Exactly same as -- 'Statistics.Sample.correlation'-pearson :: (G.Vector v (Double, Double), G.Vector v Double)+pearson :: (G.Vector v (Double, Double))         => v (Double, Double) -> Double pearson = correlation {-# INLINE pearson #-} +-- | Pearson correlation for sample of pairs. Exactly same as+-- 'Statistics.Sample.correlation'+pearson2 :: (G.Vector v Double)+         => v Double -> v Double -> Double+pearson2 = correlation2+{-# INLINE pearson2 #-}+ -- | Compute pairwise Pearson correlation between rows of a matrix pearsonMatByRow :: Matrix -> Matrix pearsonMatByRow m@@ -43,15 +52,13 @@ -- Spearman ---------------------------------------------------------------- --- | compute Spearman correlation between two samples+-- | Compute Spearman correlation between two samples spearman :: ( Ord a             , Ord b             , G.Vector v a             , G.Vector v b             , G.Vector v (a, b)             , G.Vector v Int-            , G.Vector v Double-            , G.Vector v (Double, Double)             , G.Vector v (Int, a)             , G.Vector v (Int, b)             )@@ -63,6 +70,27 @@   where     (x, y) = G.unzip xy {-# INLINE spearman #-}++-- | Compute Spearman correlation between two samples. Samples must+--   have same length.+spearman2 :: ( Ord a+            , Ord b+            , G.Vector v a+            , G.Vector v b+            , G.Vector v Int+            , G.Vector v (Int, a)+            , G.Vector v (Int, b)+            )+         => v a+         -> v b+         -> Double+spearman2 xs ys+  | nx /= ny  = error "Statistics.Correlation.spearman2: samples must have same length"+  | otherwise = pearson $ G.zip (rankUnsorted xs) (rankUnsorted ys)+  where+    nx = G.length xs+    ny = G.length ys+{-# INLINE spearman2 #-}  -- | compute pairwise Spearman correlation between rows of a matrix spearmanMatByRow :: Matrix -> Matrix
Statistics/Function.hs view
@@ -76,8 +76,8 @@ {-# INLINE indices #-}  -- | Zip a vector with its indices.-indexed :: (G.Vector v e, G.Vector v Int, G.Vector v (Int,e)) => v e -> v (Int,e)-indexed a = G.zip (indices a) a+indexed :: (G.Vector v e, G.Vector v (Int,e)) => v e -> v (Int,e)+indexed xs = G.imap (,) xs {-# INLINE indexed #-}  data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double
Statistics/Regression.hs view
@@ -41,8 +41,15 @@ --   element than the list of predictors; the last element is the --   /y/-intercept value. ----- * /R&#0178;/, the coefficient of determination (see 'rSquare' for+-- * /R²/, the coefficient of determination (see 'rSquare' for --   details).+--+-- >>> import qualified Data.Vector.Unboxed as VU+-- >>> :{+--  olsRegress [ VU.fromList [0,1,2,3]+--             ] (VU.fromList [1000, 1001, 1002, 1003])+-- :}+-- ([1.0000000000000218,999.9999999999999],1.0) olsRegress :: [Vector]               -- ^ Non-empty list of predictor vectors.  Must all have               -- the same length.  These will become the columns of@@ -65,7 +72,30 @@     lss@(n:ls) = map G.length preds olsRegress _ _ = error "no predictors given" --- | Compute the ordinary least-squares solution to /A x = b/.+-- | Compute the ordinary least-squares solution to overdetermined+--   linear system \(Ax = b\). In other words it finds+--+--   \[ \operatorname{argmin}|Ax-b|^2 \].+--+--   All columns of \(A\) must be linearly independent. It's not+--   checked function will return nonsensical result if resulting+--   linear system is poorly conditioned.+--+-- >>> import qualified Data.Vector.Unboxed as VU+-- >>> :{+--  ols (fromColumns [ VU.fromList [0,1,2,3]+--                   , VU.fromList [1,1,1,1]+--                   ]) (VU.fromList [1000, 1001, 1002, 1003])+-- :}+-- [1.0000000000000218,999.9999999999999]+--+-- >>> :{+--  ols (fromColumns [ VU.fromList [0,1,2,3]+--                   , VU.fromList [4,2,1,1]+--                   , VU.fromList [1,1,1,1]+--                   ]) (VU.fromList [1000, 1001, 1002, 1003])+-- :}+-- [1.0000000000005393,4.2290644612446807e-13,999.9999999999983] ols :: Matrix     -- ^ /A/ has at least as many rows as columns.     -> Vector     -- ^ /b/ has the same length as columns in /A/.     -> Vector@@ -92,7 +122,7 @@   where n = rows r         l = U.length b --- | Compute /R&#0178;/, the coefficient of determination that+-- | Compute /R²/, the coefficient of determination that -- indicates goodness-of-fit of a regression. -- -- This value will be 1 if the predictors fit perfectly, dropping to 0@@ -101,11 +131,19 @@         -> Vector               -- ^ Responders.         -> Vector               -- ^ Regression coefficients.         -> Double-rSquare pred resp coeff = 1 - r / t+rSquare pred resp coeff+  -- Data has zero variance. If fit is perfect we set R² to 1 else to+  -- 0. This is not perfect heuristic. Fit residuals may be nonzero+  -- due to rounding.+  | t == 0             = if r == 0 then 1 else 0+  -- If fit residuals are worse than average we simply set R² to 0+  | r2 >= 0 && r2 <= 1 = r2+  | otherwise          = 0   where-    r   = sum $ flip U.imap resp $ \i x -> square (x - p i)-    t   = sum $ flip U.map resp $ \x -> square (x - mean resp)-    p i = sum . flip U.imap coeff $ \j -> (* unsafeIndex pred i j)+    r2  = 1 - r / t+    r   = sum $ flip U.imap resp  $ \i x -> square (x - p i)+    t   = sum $ flip U.map  resp  $ \x   -> square (x - mean resp)+    p i = sum $ flip U.imap coeff $ \j x -> x * unsafeIndex pred i j  -- | Bootstrap a regression function.  Returns both the results of the -- regression and the requested confidence interval values.
Statistics/Sample.hs view
@@ -20,6 +20,7 @@     , range      -- * Statistics of location+    , expectation     , mean     , welfordMean     , meanWeighted@@ -54,17 +55,20 @@     -- * Joint distributions     , covariance     , correlation+    , covariance2+    , correlation2     , pair     -- * References     -- $references     ) where -import Statistics.Function (minMax)+import Statistics.Function (minMax,square) import Statistics.Sample.Internal (robustSumVar, sum) import Statistics.Types.Internal  (Sample,WeightedSample) import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U+import Numeric.Sum (kbn, Summation(zero,add))  -- Operator ^ will be overridden import Prelude hiding ((^), sum)@@ -76,9 +80,17 @@     where (lo , hi) = minMax s {-# INLINE range #-} +-- | /O(n)/ Compute expectation of function over for sample. This is+--   simply @mean . map f@ but won't create intermediate vector.+expectation :: (G.Vector v a) => (a -> Double) -> v a -> Double+expectation f xs = kbn (G.foldl' (\s -> add s . f) zero xs)+                 / fromIntegral (G.length xs)+{-# INLINE expectation #-}+ -- | /O(n)/ Arithmetic mean.  This uses Kahan-Babuška-Neumaier -- summation, so is more accurate than 'welfordMean' unless the input--- values are very large.+-- values are very large. This function is not subject to stream+-- fusion. mean :: (G.Vector v Double) => v Double -> Double mean xs = sum xs / fromIntegral (G.length xs) {-# SPECIALIZE mean :: U.Vector Double -> Double #-}@@ -122,7 +134,7 @@  -- | /O(n)/ Geometric mean of a sample containing no negative values. geometricMean :: (G.Vector v Double) => v Double -> Double-geometricMean = exp . mean . G.map log+geometricMean = exp . expectation log {-# INLINE geometricMean #-}  -- | Compute the /k/th central moment of a sample.  The central moment@@ -138,7 +150,7 @@     | a < 0  = error "Statistics.Sample.centralMoment: negative input"     | a == 0 = 1     | a == 1 = 0-    | otherwise = sum (G.map go xs) / fromIntegral (G.length xs)+    | otherwise = expectation go xs   where     go x = (x-m) ^ a     m    = mean xs@@ -354,40 +366,77 @@  -- | Covariance of sample of pairs. For empty sample it's set to --   zero-covariance :: (G.Vector v (Double,Double), G.Vector v Double)+covariance :: (G.Vector v (Double,Double))            => v (Double,Double)            -> Double covariance xy   | n == 0    = 0-  | otherwise = mean $ G.zipWith (*)-                         (G.map (\x -> x - muX) xs)-                         (G.map (\y -> y - muY) ys)+  | otherwise = expectation (\(x,y) -> (x - muX)*(y - muY)) xy   where-    n       = G.length xy-    (xs,ys) = G.unzip xy-    muX     = mean xs-    muY     = mean ys+    n   = G.length xy+    muX = expectation fst xy+    muY = expectation snd xy {-# SPECIALIZE covariance :: U.Vector (Double,Double) -> Double #-} {-# SPECIALIZE covariance :: V.Vector (Double,Double) -> Double #-}  -- | Correlation coefficient for sample of pairs. Also known as --   Pearson's correlation. For empty sample it's set to zero.-correlation :: (G.Vector v (Double,Double), G.Vector v Double)+correlation :: (G.Vector v (Double,Double))            => v (Double,Double)            -> Double correlation xy   | n == 0    = 0   | otherwise = cov / sqrt (varX * varY)   where-    n       = G.length xy-    (xs,ys) = G.unzip xy-    (muX,varX) = meanVariance xs-    (muY,varY) = meanVariance ys-    cov = mean $ G.zipWith (*)-            (G.map (\x -> x - muX) xs)-            (G.map (\y -> y - muY) ys)+    n    = G.length xy+    muX  = expectation (\(x,_) -> x) xy+    muY  = expectation (\(_,y) -> y) xy+    varX = expectation (\(x,_) -> square (x - muX))    xy+    varY = expectation (\(_,y) -> square (y - muY))    xy+    cov  = expectation (\(x,y) -> (x - muX)*(y - muY)) xy {-# SPECIALIZE correlation :: U.Vector (Double,Double) -> Double #-} {-# SPECIALIZE correlation :: V.Vector (Double,Double) -> Double #-}+++-- | Covariance of two samples. Both vectors must be of the same+--   length. If both are empty it's set to zero+covariance2 :: (G.Vector v Double)+           => v Double+           -> v Double+           -> Double+covariance2 xs ys+  | nx /= ny  = error $ "Statistics.Sample.covariance2: both samples must have same length"+  | nx == 0   = 0+  | otherwise = sum (G.zipWith (\x y -> (x - muX)*(y - muY)) xs ys)+              / fromIntegral nx+  where+    nx  = G.length xs+    ny  = G.length ys+    muX = mean xs+    muY = mean ys+{-# SPECIALIZE covariance2 :: U.Vector Double -> U.Vector Double -> Double #-}+{-# SPECIALIZE covariance2 :: V.Vector Double -> V.Vector Double -> Double #-}++-- | Correlation coefficient for two samples. Both vector must have+--   same length Also known as Pearson's correlation. For empty sample+--   it's set to zero.+correlation2 :: (G.Vector v Double)+             => v Double+             -> v Double+             -> Double+correlation2 xs ys+  | nx /= ny  = error $ "Statistics.Sample.correlation2: both samples must have same length"+  | nx == 0   = 0+  | otherwise = cov / sqrt (varX * varY)+  where+    nx         = G.length xs+    ny         = G.length ys+    (muX,varX) = meanVariance xs+    (muY,varY) = meanVariance ys+    cov = sum (G.zipWith (\x y -> (x - muX)*(y - muY)) xs ys)+        / fromIntegral nx+{-# SPECIALIZE correlation2 :: U.Vector Double -> U.Vector Double -> Double #-}+{-# SPECIALIZE correlation2 :: V.Vector Double -> V.Vector Double -> Double #-}   -- | Pair two samples. It's like 'G.zip' but requires that both
Statistics/Test/Internal.hs view
@@ -8,6 +8,7 @@ import Data.Ord import           Data.Vector.Generic           ((!)) import qualified Data.Vector.Generic         as G+import qualified Data.Vector.Unboxed         as U import qualified Data.Vector.Generic.Mutable as M import Statistics.Function @@ -27,15 +28,16 @@ --   In case of ties average of ranks of equal elements is assigned --   to each ----- >>> rank (==) (fromList [10,20,30::Int])--- > fromList [1.0,2.0,3.0]+-- >>> import qualified Data.Vector.Unboxed as VU+-- >>> rank (==) (VU.fromList [10,20,30::Int])+-- [1.0,2.0,3.0] ----- >>> rank (==) (fromList [10,10,10,30::Int])--- > fromList [2.0,2.0,2.0,4.0]-rank :: (G.Vector v a, G.Vector v Double)+-- >>> rank (==) (VU.fromList [10,10,10,30::Int])+-- [2.0,2.0,2.0,4.0]+rank :: (G.Vector v a)      => (a -> a -> Bool)        -- ^ Equivalence relation      -> v a                     -- ^ Vector to rank-     -> v Double+     -> U.Vector Double rank eq vec = G.unfoldr go (Rank 0 (-1) 1 vec)   where     go (Rank 0 _ r v)@@ -58,11 +60,10 @@ rankUnsorted :: ( Ord a                 , G.Vector v a                 , G.Vector v Int-                , G.Vector v Double                 , G.Vector v (Int, a)                 )              => v a-             -> v Double+             -> U.Vector Double rankUnsorted xs = G.create $ do     -- Put ranks into their original positions     -- NOTE: backpermute will do wrong thing
Statistics/Test/StudentT.hs view
@@ -71,7 +71,7 @@ -- | Paired two-sample t-test. Two samples are paired in a -- within-subject design. Returns @Nothing@ if sample size is not -- sufficient.-pairedTTest :: forall v. (G.Vector v (Double, Double), G.Vector v Double)+pairedTTest :: forall v. (G.Vector v (Double, Double))             => PositionTest          -- ^ one- or two-tailed test             -> v (Double, Double)    -- ^ paired samples             -> Maybe (Test StudentT)
Statistics/Types.hs view
@@ -94,7 +94,7 @@ -- second from @1 - CL@ or significance level. -- -- >>> cl95--- mkCLFromSignificance 0.05+-- mkCLFromSignificance 5.0e-2 -- -- Prior to 0.14 confidence levels were passed to function as plain -- @Doubles@. Use 'mkCL' to convert them to @CL@.@@ -134,7 +134,7 @@ --   exception if parameter is out of [0,1] range -- -- >>> mkCL 0.95    -- same as cl95--- mkCLFromSignificance 0.05+-- mkCLFromSignificance 5.0000000000000044e-2 mkCL :: (Ord a, Num a) => a -> CL a mkCL   = fromMaybe (error "Statistics.Types.mkCL: probability is out if [0,1] range")@@ -144,7 +144,7 @@ --   parameter is out of [0,1] range -- -- >>> mkCLE 0.95    -- same as cl95--- Just (mkCLFromSignificance 0.05)+-- Just (mkCLFromSignificance 5.0000000000000044e-2) mkCLE :: (Ord a, Num a) => a -> Maybe (CL a) mkCLE p   | p >= 0 && p <= 1 = Just $ CL (1 - p)@@ -155,7 +155,7 @@ --   throw exception if parameter is out of [0,1] range -- -- >>> mkCLFromSignificance 0.05    -- same as cl95--- mkCLFromSignificance 0.05+-- mkCLFromSignificance 5.0e-2 mkCLFromSignificance :: (Ord a, Num a) => a -> CL a mkCLFromSignificance = fromMaybe (error errMkCL) . mkCLFromSignificanceE @@ -163,7 +163,7 @@ --   parameter is out of [0,1] range -- -- >>> mkCLFromSignificanceE 0.05    -- same as cl95--- Just (mkCLFromSignificance 0.05)+-- Just (mkCLFromSignificance 5.0e-2) mkCLFromSignificanceE :: (Ord a, Num a) => a -> Maybe (CL a) mkCLFromSignificanceE p   | p >= 0 && p <= 1 = Just $ CL p@@ -299,6 +299,7 @@ -- >                       , confIntUDX = 6 -- >                       , confIntCL  = cl95 -- >                       }+-- >          } -- -- Prior to statistics 0.14 @Estimate@ data type used following definition: --
+ bench-papi/Bench.hs view
@@ -0,0 +1,14 @@+-- |+-- Here we reexport definitions of tasty-bench+module Bench+  ( whnf+  , nf+  , nfIO+  , whnfIO+  , bench+  , bgroup+  , defaultMain+  , benchIngredients+  ) where++import Test.Tasty.PAPI
+ bench-time/Bench.hs view
@@ -0,0 +1,14 @@+-- |+-- Here we reexport definitions of tasty-bench+module Bench+  ( whnf+  , nf+  , nfIO+  , whnfIO+  , bench+  , bgroup+  , defaultMain+  , benchIngredients+  ) where++import Test.Tasty.Bench
benchmark/bench.hs view
@@ -1,24 +1,26 @@-import Control.Monad.ST (runST)-import Criterion.Main import Data.Complex import Statistics.Sample import Statistics.Transform-import Statistics.Correlation.Pearson+import Statistics.Correlation import System.Random.MWC-import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as MVU +import Bench + -- Test sample-sample :: U.Vector Double-sample = runST $ flip uniformVector 10000 =<< create+sample :: VU.Vector Double+sample = VU.create $ do g <- create+                        MVU.replicateM 10000 (uniform g)  -- Weighted test sample-sampleW :: U.Vector (Double,Double)-sampleW = U.zip sample (U.reverse sample)+sampleW :: VU.Vector (Double,Double)+sampleW = VU.zip sample (VU.reverse sample)  -- Complex vector for FFT tests-sampleC :: U.Vector (Complex Double)-sampleC = U.zipWith (:+) sample (U.reverse sample)+sampleC :: VU.Vector (Complex Double)+sampleC = VU.zipWith (:+) sample (VU.reverse sample)   -- Simple benchmark for functions from Statistics.Sample@@ -37,9 +39,11 @@     , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample     , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW       -- Correlation-    , bench "pearson"          $ nf (\x -> pearson (U.reverse sample) x) sample-    , bench "pearson'"          $ nf (\x -> pearson' (U.reverse sample) x) sample-    , bench "pearsonFast"      $ nf (\x -> pearsonFast (U.reverse sample) x) sample+    , bench "pearson"          $ nf pearson     sampleW+    , bench "covariance"       $ nf covariance  sampleW+    , bench "correlation"      $ nf correlation sampleW+    , bench "covariance2"      $ nf (covariance2  sample) sample+    , bench "correlation2"     $ nf (correlation2 sample) sample       -- Other     , bench "stdDev"           $ nf (\x -> stdDev x)           sample     , bench "skewness"         $ nf (\x -> skewness x)         sample@@ -52,17 +56,17 @@     ]   , bgroup "FFT"     [ bgroup "fft"-      [ bench  (show n) $ whnf fft   (U.take n sampleC) | n <- fftSizes ]+      [ bench  (show n) $ whnf fft   (VU.take n sampleC) | n <- fftSizes ]     , bgroup "ifft"-      [ bench  (show n) $ whnf ifft  (U.take n sampleC) | n <- fftSizes ]+      [ bench  (show n) $ whnf ifft  (VU.take n sampleC) | n <- fftSizes ]     , bgroup "dct"-      [ bench  (show n) $ whnf dct   (U.take n sample)  | n <- fftSizes ]+      [ bench  (show n) $ whnf dct   (VU.take n sample)  | n <- fftSizes ]     , bgroup "dct_"-      [ bench  (show n) $ whnf dct_  (U.take n sampleC) | n <- fftSizes ]+      [ bench  (show n) $ whnf dct_  (VU.take n sampleC) | n <- fftSizes ]     , bgroup "idct"-      [ bench  (show n) $ whnf idct  (U.take n sample)  | n <- fftSizes ]+      [ bench  (show n) $ whnf idct  (VU.take n sample)  | n <- fftSizes ]     , bgroup "idct_"-      [ bench  (show n) $ whnf idct_ (U.take n sampleC) | n <- fftSizes ]+      [ bench  (show n) $ whnf idct_ (VU.take n sampleC) | n <- fftSizes ]     ]   ]
changelog.md view
@@ -1,3 +1,25 @@+## Changes in 0.16.3.0++ * `S.Sample.correlation`, `S.Sample.covariance`,+   `S.Correlation.pearson` do not allocate temporary arrays.++ * Variants of correlation which take two vectors as input are added:+   `S.Sample.correlation2`, `S.Sample.covariance2`, `S.Correlation.pearson2`,+   `S.Correlation.spearman2`.++ * Contexts for `S.Function.indexed`, `S.Correlation.spearman`, `S.pairedTTest`,+   `S.Sample.correlation`, `S.Sample.covariance`, reduced.++ * Computation of `rSquare` in linear regression has special case for case when+   data variation is 0.++ * Doctests added.++ * Benchmarks using `tasty-bench` and `tasty-papi` added.++ * Spurious test failures fixed.++ ## Changes in 0.16.2.1   * Unnecessary constraint dropped from `tStatisticsPaired`.
statistics.cabal view
@@ -1,5 +1,8 @@+cabal-version:  3.0+build-type:     Simple+ name:           statistics-version:        0.16.2.1+version:        0.16.3.0 synopsis:       A library of statistical types, data, and functions description:   This library provides a number of common functions and types useful@@ -22,7 +25,7 @@   * Common statistical tests for significant differences between     samples. -license:        BSD2+license:        BSD-2-Clause license-file:   LICENSE homepage:       https://github.com/haskell/statistics bug-reports:    https://github.com/haskell/statistics/issues@@ -30,32 +33,39 @@ maintainer:     Alexey Khudaykov <alexey.skladnoy@gmail.com> copyright:      2009-2014 Bryan O'Sullivan category:       Math, Statistics-build-type:     Simple-cabal-version:  >= 1.10+ extra-source-files:   README.markdown-  benchmark/bench.hs   changelog.md   examples/kde/KDE.hs   examples/kde/data/faithful.csv   examples/kde/kde.html   examples/kde/kde.tpl-  tests/Tests/Math/Tables.hs-  tests/Tests/Math/gen.py   tests/utils/Makefile   tests/utils/fftw.c  tested-with:-    GHC ==8.4.4-    GHC ==8.6.5-    GHC ==8.8.4-    GHC ==8.10.7-    GHC ==9.0.2-    GHC ==9.2.8-    GHC ==9.4.6-    GHC ==9.6.2+  GHC ==8.4.4+   || ==8.6.5+   || ==8.8.4+   || ==8.10.7+   || ==9.0.2+   || ==9.2.8+   || ==9.4.8+   || ==9.6.6+   || ==9.8.4+   || ==9.10.1 +source-repository head+  type:     git+  location: https://github.com/haskell/statistics +flag BenchPAPI+  Description: Enable building of benchmarks which use instruction counters.+               It requires libpapi and only works on Linux so it's protected by flag+  Default: False+  Manual:  True+ library   default-language: Haskell2010   exposed-modules:@@ -176,6 +186,49 @@                , vector                , vector-algorithms -source-repository head-  type:     git-  location: https://github.com/haskell/statistics+test-suite statistics-doctests+  default-language: Haskell2010+  type:             exitcode-stdio-1.0+  hs-source-dirs:   tests+  main-is:          doctest.hs+  if impl(ghcjs) || impl(ghc < 8.0)+    Buildable: False+  -- Linker on macos prints warnings to console which confuses doctests.+  -- We simply disable doctests on ma for older GHC+  -- > warning: -single_module is obsolete+  if os(darwin) && impl(ghc < 9.6)+    buildable: False+  build-depends:+            base       -any+          , statistics -any+          , doctest    >=0.15 && <0.24++-- We want to be able to build benchmarks using both tasty-bench and tasty-papi.+-- They have similar API so we just create two shim modules which reexport+-- definitions from corresponding library and pick one in cabal file.+common bench-stanza+  ghc-options:      -Wall+  default-language: Haskell2010+  build-depends: base < 5+               , vector          >= 0.12.3+               , statistics+               , mwc-random+               , tasty           >=1.3.1++benchmark statistics-bench+  import:         bench-stanza+  type:           exitcode-stdio-1.0+  hs-source-dirs: benchmark bench-time+  main-is:        bench.hs+  Other-modules:  Bench+  build-depends:  tasty-bench >= 0.3++benchmark statistics-bench-papi+  import:         bench-stanza+  type:           exitcode-stdio-1.0+  if impl(ghcjs) || !flag(BenchPAPI)+     buildable: False+  hs-source-dirs: benchmark bench-papi+  main-is:        bench.hs+  Other-modules:  Bench+  build-depends:  tasty-papi >= 0.1.2
tests/Tests/Correlation.hs view
@@ -102,11 +102,11 @@         , not (isNaN c3)         , not (isNaN c4)         ]-  ==> ( counterexample (show sample0)-      $ counterexample (show sample1)-      $ counterexample (show sample2)-      $ counterexample (show sample3)-      $ counterexample (show sample4)+  ==> ( counterexample ("S0 = " ++ show sample0)+      $ counterexample ("S1 = " ++ show sample1)+      $ counterexample ("S2 = " ++ show sample2)+      $ counterexample ("S3 = " ++ show sample3)+      $ counterexample ("S4 = " ++ show sample4)       $ counterexample (show (c1,c2,c3,c4))       $ and [ c1 == c2             , c1 == c3@@ -117,8 +117,8 @@     -- We need to stretch sample into [-10 .. 10] range to avoid     -- problems with under/overflows etc.     stretch xs-      | a == b = xs-      | otherwise = [ (x - a - 10) * 20 / (a - b) | x <- xs ]+      | a == b    = xs+      | otherwise = [ ((x - a)/(b - a) - 0.5) * 20 | x <- xs ]       where         a = minimum xs         b = maximum xs
tests/Tests/ExactDistribution.hs view
@@ -1,8 +1,9 @@-{-# LANGUAGE BangPatterns,-             FlexibleInstances,-             FlexibleContexts,-             ScopedTypeVariables-  #-}+{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE FlexibleContexts    #-}+{-# LANGUAGE FlexibleInstances   #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications    #-}+{-# LANGUAGE TypeFamilies        #-} -- | -- Module    : Tests.ExactDistribution -- Copyright : (c) 2022 Lorenz Minder@@ -64,8 +65,6 @@     , ExactHypergeomDistr(..)      -- * Linking to production distributions-    , ProductionProbFuncs(..)-    , productionProbFuncs     , ProductionLinkage      -- * Individual test routines@@ -88,6 +87,7 @@ import Test.Tasty.QuickCheck            (testProperty) import Test.QuickCheck as QC import Numeric.MathFunctions.Comparison (relativeError)+import Numeric.MathFunctions.Constants  (m_tiny)  import Statistics.Distribution import Statistics.Distribution.Binomial@@ -295,92 +295,84 @@ -- ---------------------------------------------------------------- --- | Distribution evaluation functions.------ This is used to store a-data ProductionProbFuncs = ProductionProbFuncs {-        prodProb            :: Int -> Double-    ,   prodCumulative      :: Double -> Double-    ,   prodComplCumulative :: Double -> Double-    }--productionProbFuncs :: (DiscreteDistr a) => a -> ProductionProbFuncs-productionProbFuncs d = ProductionProbFuncs {-        prodProb = probability d-    ,   prodCumulative = cumulative d-    ,   prodComplCumulative = complCumulative d-    }--class (ExactDiscreteDistr a) => ProductionLinkage a where-    productionLinkage :: a -> ProductionProbFuncs+class (ExactDiscreteDistr a, DiscreteDistr (ProdDistrib a)+      ) => ProductionLinkage a where+  type ProdDistrib a+  toProd :: a -> ProdDistrib a  instance ProductionLinkage ExactBinomialDistr where-    productionLinkage (ExactBD n p) =-        let d = binomial (fromIntegral n) (fromRational p)-        in  productionProbFuncs d+  type ProdDistrib ExactBinomialDistr = BinomialDistribution+  toProd (ExactBD n p) = binomial (fromIntegral n) (fromRational p)  instance ProductionLinkage ExactDiscreteUniformDistr where-    productionLinkage (ExactDU lower upper) =-        let d = discreteUniformAB (fromIntegral lower) (fromIntegral upper)-        in  productionProbFuncs d+  type ProdDistrib ExactDiscreteUniformDistr = DiscreteUniform+  toProd (ExactDU lower upper) = discreteUniformAB (fromIntegral lower) (fromIntegral upper)  instance ProductionLinkage ExactGeometricDistr where-    productionLinkage (ExactGeom p) =-        let d = geometric $ fromRational p-        in  productionProbFuncs d+  type ProdDistrib ExactGeometricDistr = GeometricDistribution+  toProd (ExactGeom p) = geometric $ fromRational p  instance ProductionLinkage ExactHypergeomDistr where-    productionLinkage (ExactHG nK nN n) =-        let d = hypergeometric (fromIntegral nK) (fromIntegral nN) (fromIntegral n)-        in  productionProbFuncs d+  type ProdDistrib ExactHypergeomDistr = HypergeometricDistribution+  toProd (ExactHG nK nN n) =+    hypergeometric (fromIntegral nK) (fromIntegral nN) (fromIntegral n) + ---------------------------------------------------------------- -- Tests ---------------------------------------------------------------- +-- Compare that probabilities agree. If they are denormalized just+-- return True. You can't say much about precision+probabilityAgree :: Double -> Double -> Double -> Bool+probabilityAgree tol pe pa+  | pa < 0      = False+  | pe < 0      = False+  | pe < m_tiny = True+  | otherwise   = relativeError pe pa < tol+ -- Check production probability mass function accuracy. -- -- Inputs: tolerance (max relative error) and test case-pmfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Bool-pmfMatch tol (TestCase dExact k) =-    let dProd = productionLinkage dExact-        pe = fromRational $ exactProb dExact k-        pa = prodProb dProd k'-        k' = fromIntegral k-    in  relativeError pe pa < tol+pmfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Property+pmfMatch tol (TestCase dExact k)+  = counterexample ("Exact  = " ++ show pe)+  $ counterexample ("Approx = " ++ show pa)+  $ probabilityAgree tol pe pa+  where+    pe = fromRational $ exactProb dExact k+    pa = probability (toProd dExact) (fromIntegral k)  -- Check production cumulative probability function accuracy. -- -- Inputs:  tolerance (max relative error) and test case. cdfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Bool-cdfMatch tol (TestCase dExact k) =-    let dProd = productionLinkage dExact-        pe = fromRational $ exactCumulative dExact k-        pa = prodCumulative dProd k'-        k' = fromIntegral k-    in  relativeError pe pa < tol+cdfMatch tol (TestCase dExact k)+  = probabilityAgree tol pe pa+  where+    pe = fromRational $ exactCumulative dExact k+    pa = cumulative (toProd dExact) (fromIntegral k)  -- Check production complement cumulative function accuracy. -- -- Inputs:  tolerance (max relative error) and test case. complCdfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Bool-complCdfMatch tol (TestCase dExact k) =-    let dProd = productionLinkage dExact-        pe = fromRational $ 1 - exactCumulative dExact k-        pa = prodComplCumulative dProd k'-        k' = fromIntegral k-    in  relativeError pe pa < tol+complCdfMatch tol (TestCase dExact k)+  = probabilityAgree tol pe pa+  where+    pe = fromRational $ 1 - exactCumulative dExact k+    pa = complCumulative (toProd dExact) (fromIntegral k)  -- Phantom type to encode an exact distribution. data Tag a = Tag -distTests :: (Show a, ProductionLinkage a, Arbitrary (TestCase a)) =>+distTests :: forall a. (Show a, ProductionLinkage a, Arbitrary (TestCase a)) =>     Tag a -> String -> Double -> TestTree distTests (Tag :: Tag a) name tol =-    testGroup ("Exact tests for " ++ name) [-        testProperty "PMF match" $ ((pmfMatch tol) :: TestCase a -> Bool)-    ,   testProperty "CDF match" $ ((cdfMatch tol) :: TestCase a -> Bool)-    ,   testProperty "1 - CDF match" $ ((complCdfMatch tol) :: TestCase a -> Bool)+  testGroup ("Exact tests for " ++ name)+    [ testProperty "PMF match"     $ pmfMatch      @a tol+    , testProperty "CDF match"     $ cdfMatch      @a tol+    , testProperty "1 - CDF match" $ complCdfMatch @a tol     ]  @@ -388,9 +380,8 @@  exactDistributionTests :: TestTree exactDistributionTests = testGroup "Test distributions against exact"-  [-    distTests (Tag :: Tag ExactBinomialDistr)       "Binomial"          1.0e-12-  , distTests (Tag :: Tag ExactDiscreteUniformDistr) "DiscreteUniform"  1.0e-12-  , distTests (Tag :: Tag ExactGeometricDistr)      "Geometric"         1.0e-13-  , distTests (Tag :: Tag ExactHypergeomDistr)      "Hypergeometric"    1.0e-12+  [ distTests (Tag @ExactBinomialDistr)        "Binomial"         1.0e-12+  , distTests (Tag @ExactDiscreteUniformDistr) "DiscreteUniform"  1.0e-12+  , distTests (Tag @ExactGeometricDistr)       "Geometric"        1.0e-13+  , distTests (Tag @ExactHypergeomDistr)       "Hypergeometric"   1.0e-12   ]
− tests/Tests/Math/Tables.hs 
@@ -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)-  ]
− tests/Tests/Math/gen.py 
@@ -1,51 +0,0 @@-#!/usr/bin/python-"""-"""--from mpmath import *--def printListLiteral(lines) :-    print "  [" + "\n  , ".join(lines) + "\n  ]"--################################################################-# Generate header-print "module Tests.Math.Tables where"-print--################################################################-## Generate table for logGamma-print "tableLogGamma :: [(Double,Double)]"-print "tableLogGamma ="--gammaArg = [ 1.25e-6, 6.82e-5, 2.46e-4, 8.8e-4,  3.12e-3, 2.67e-2,-             7.77e-2, 0.234,   0.86,    1.34,    1.89,    2.45,-             3.65,    4.56,    6.66,    8.25,    11.3,    25.6,-             50.4,    123.3,   487.4,   853.4,   2923.3,  8764.3,-             1.263e4, 3.45e4,  8.234e4, 2.348e5, 8.343e5, 1.23e6,-             ]-printListLiteral(-    [ '(%.15f, %.20g)' % (x, log(gamma(x))) for x in gammaArg ]-    )---################################################################-## Generate table for incompleteBeta--print "tableIncompleteBeta :: [(Double,Double,Double,Double)]"-print "tableIncompleteBeta ="--incompleteBetaArg = [-    (2,    3,    0.03),-    (2,    3,    0.23),-    (2,    3,    0.76),-    (4,    2.3,  0.89),-    (1,    1,    0.55),-    (0.3,  12.2, 0.11),-    (13.1, 9.8,  0.12),-    (13.1, 9.8,  0.42),-    (13.1, 9.8,  0.92),-    ]-printListLiteral(-    [ '(%.15f, %.15f, %.15f, %.20g)' % (p,q,x, betainc(p,q,0,x, regularized=True))-      for (p,q,x) in incompleteBetaArg-      ])
+ tests/doctest.hs view
@@ -0,0 +1,5 @@+import Test.DocTest (doctest)++main :: IO ()+main = doctest ["-XHaskell2010", "Statistics"]+