statistics 0.15.2.0 → 0.16.5.0
raw patch · 61 files changed
Files
- README.markdown +3/−10
- Setup.lhs +0/−3
- Statistics/Correlation.hs +34/−6
- Statistics/Correlation/Kendall.hs +1/−6
- Statistics/Distribution.hs +27/−33
- Statistics/Distribution/Binomial.hs +23/−5
- Statistics/Distribution/CauchyLorentz.hs +31/−11
- Statistics/Distribution/DiscreteUniform.hs +8/−1
- Statistics/Distribution/Exponential.hs +6/−9
- Statistics/Distribution/FDistribution.hs +27/−6
- Statistics/Distribution/Gamma.hs +6/−1
- Statistics/Distribution/Geometric.hs +16/−12
- Statistics/Distribution/Hypergeometric.hs +17/−4
- Statistics/Distribution/Laplace.hs +3/−3
- Statistics/Distribution/Lognormal.hs +172/−0
- Statistics/Distribution/NegativeBinomial.hs +188/−0
- Statistics/Distribution/Normal.hs +18/−10
- Statistics/Distribution/Poisson.hs +12/−1
- Statistics/Distribution/Poisson/Internal.hs +3/−3
- Statistics/Distribution/StudentT.hs +4/−3
- Statistics/Distribution/Transform.hs +0/−2
- Statistics/Distribution/Uniform.hs +6/−6
- Statistics/Distribution/Weibull.hs +224/−0
- Statistics/Function.hs +2/−2
- Statistics/Internal.hs +0/−2
- Statistics/Quantile.hs +1/−11
- Statistics/Regression.hs +45/−8
- Statistics/Resampling.hs +0/−4
- Statistics/Resampling/Bootstrap.hs +2/−13
- Statistics/Sample.hs +75/−24
- Statistics/Sample/Histogram.hs +5/−3
- Statistics/Sample/Internal.hs +7/−2
- Statistics/Test/Bartlett.hs +99/−0
- Statistics/Test/ChiSquared.hs +14/−8
- Statistics/Test/Internal.hs +9/−8
- Statistics/Test/KolmogorovSmirnov.hs +1/−1
- Statistics/Test/KruskalWallis.hs +0/−1
- Statistics/Test/Levene.hs +153/−0
- Statistics/Test/MannWhitneyU.hs +0/−1
- Statistics/Test/StudentT.hs +5/−5
- Statistics/Test/WilcoxonT.hs +0/−1
- Statistics/Types.hs +6/−14
- bench-papi/Bench.hs +14/−0
- bench-time/Bench.hs +14/−0
- benchmark/Main.hs +77/−0
- benchmark/bench.hs +0/−71
- changelog.md +98/−12
- statistics.cabal +95/−33
- tests/Tests/Correlation.hs +16/−12
- tests/Tests/Distribution.hs +77/−32
- tests/Tests/ExactDistribution.hs +387/−0
- tests/Tests/Helpers.hs +2/−2
- tests/Tests/Math/Tables.hs +0/−47
- tests/Tests/Math/gen.py +0/−51
- tests/Tests/Matrix.hs +14/−4
- tests/Tests/Matrix/Types.hs +17/−4
- tests/Tests/Orphanage.hs +14/−4
- tests/Tests/Parametric.hs +128/−6
- tests/Tests/Serialization.hs +6/−0
- tests/Tests/Transform.hs +0/−1
- tests/doctest.hs +5/−0
README.markdown view
@@ -18,18 +18,11 @@ # Get involved! Please report bugs via the-[github issue tracker](https://github.com/bos/statistics/issues).--Master [git mirror](https://github.com/bos/statistics):--* `git clone git://github.com/bos/statistics.git`--There's also a [Mercurial mirror](https://bitbucket.org/bos/statistics):--* `hg clone https://bitbucket.org/bos/statistics`+[github issue tracker](https://github.com/haskell/statistics/issues). -(You can create and contribute changes using either Mercurial or git.)+Master [git mirror](https://github.com/haskell/statistics): +* `git clone git://github.com/haskell/statistics.git` # Authors
− 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,12 +27,19 @@ -- | 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 #-} --- | Compute pairwise pearson correlation between rows of a matrix+-- | 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 = generateSym (rows 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) )@@ -64,7 +71,28 @@ (x, y) = G.unzip xy {-# INLINE spearman #-} --- | compute pairwise spearman correlation between rows of a matrix+-- | 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 spearmanMatByRow = pearsonMatByRow . fromRows . fmap rankUnsorted . toRows
Statistics/Correlation/Kendall.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns, CPP, FlexibleContexts #-}+{-# LANGUAGE BangPatterns, FlexibleContexts #-} -- | -- Module : Statistics.Correlation.Kendall --@@ -130,11 +130,6 @@ _ -> do GM.unsafeWrite src iIns eLow wroteLow low (iLow+1) high iHigh eHigh (iIns+1) {-# INLINE merge #-}--#if !MIN_VERSION_base(4,6,0)-modifySTRef' :: STRef s a -> (a -> a) -> ST s ()-modifySTRef' = modifySTRef-#endif -- $references --
Statistics/Distribution.hs view
@@ -29,18 +29,15 @@ , ContGen(..) , DiscreteGen(..) , genContinuous- , genContinous -- * Helper functions , findRoot , sumProbabilities ) where -import Control.Applicative ((<$>), Applicative(..))-import Control.Monad.Primitive (PrimMonad,PrimState) import Prelude hiding (sum)-import Statistics.Function (square)+import Statistics.Function (square) import Statistics.Sample.Internal (sum)-import System.Random.MWC (Gen, uniform)+import System.Random.Stateful (StatefulGen, uniformDouble01M) import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Generic as G @@ -50,13 +47,13 @@ class Distribution d where -- | Cumulative distribution function. The probability that a -- random variable /X/ is less or equal than /x/,- -- i.e. P(/X/≤/x/). Cumulative should be defined for+ -- i.e. P(/X/≤/x/). Cumulative should be defined for -- infinities as well: -- -- > cumulative d +∞ = 1 -- > cumulative d -∞ = 0 cumulative :: d -> Double -> Double-+ cumulative d x = 1 - complCumulative d x -- | One's complement of cumulative distribution: -- -- > complCumulative d x = 1 - cumulative d x@@ -67,17 +64,18 @@ -- encouraged to provide more precise implementation. complCumulative :: d -> Double -> Double complCumulative d x = 1 - cumulative d x+ {-# MINIMAL (cumulative | complCumulative) #-} + -- | Discrete probability distribution. class Distribution d => DiscreteDistr d where -- | Probability of n-th outcome. probability :: d -> Int -> Double probability d = exp . logProbability d- -- | Logarithm of probability of n-th outcome logProbability :: d -> Int -> Double logProbability d = log . probability d-+ {-# MINIMAL (probability | logProbability) #-} -- | Continuous probability distribution. --@@ -86,25 +84,23 @@ class Distribution d => ContDistr d where -- | Probability density function. Probability that random -- variable /X/ lies in the infinitesimal interval- -- [/x/,/x+/δ/x/) equal to /density(x)/⋅δ/x/+ -- [/x/,/x+/δ/x/) equal to /density(x)/⋅δ/x/ density :: d -> Double -> Double density d = exp . logDensity d-+ -- | Natural logarithm of density.+ logDensity :: d -> Double -> Double+ logDensity d = log . density d -- | Inverse of the cumulative distribution function. The value- -- /x/ for which P(/X/≤/x/) = /p/. If probability is outside+ -- /x/ for which P(/X/≤/x/) = /p/. If probability is outside -- of [0,1] range function should call 'error' quantile :: d -> Double -> Double-+ quantile d x = complQuantile d (1 - x) -- | 1-complement of @quantile@: -- -- > complQuantile x ≡ quantile (1 - x) complQuantile :: d -> Double -> Double complQuantile d x = quantile d (1 - x)-- -- | Natural logarithm of density.- logDensity :: d -> Double -> Double- logDensity d = log . density d-+ {-# MINIMAL (density | logDensity), (quantile | complQuantile) #-} -- | Type class for distributions with mean. 'maybeMean' should return -- 'Nothing' if it's undefined for current value of data@@ -126,9 +122,10 @@ -- Minimal complete definition is 'maybeVariance' or 'maybeStdDev' class MaybeMean d => MaybeVariance d where maybeVariance :: d -> Maybe Double- maybeVariance d = (*) <$> x <*> x where x = maybeStdDev d+ maybeVariance = fmap square . maybeStdDev maybeStdDev :: d -> Maybe Double- maybeStdDev = fmap sqrt . maybeVariance+ maybeStdDev = fmap sqrt . maybeVariance+ {-# MINIMAL (maybeVariance | maybeStdDev) #-} -- | Type class for distributions with variance. If distribution have -- finite variance for all valid parameter values it should be@@ -140,7 +137,9 @@ variance d = square (stdDev d) stdDev :: d -> Double stdDev = sqrt . variance+ {-# MINIMAL (variance | stdDev) #-} + -- | Type class for distributions with entropy, meaning Shannon entropy -- in the case of a discrete distribution, or differential entropy in the -- case of a continuous one. 'maybeEntropy' should return 'Nothing' if@@ -161,35 +160,30 @@ -- | Generate discrete random variates which have given -- distribution. class Distribution d => ContGen d where- genContVar :: PrimMonad m => d -> Gen (PrimState m) -> m Double+ genContVar :: (StatefulGen g m) => d -> g -> 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+ genDiscreteVar :: (StatefulGen g m) => d -> g -> m Int -- | Estimate distribution from sample. First parameter in sample is -- distribution type and second is element type. class FromSample d a where- -- | Estimate distribution from sample. Returns nothing is there's- -- not enough data to estimate or sample clearly doesn't come from- -- distribution in question. For example if there's negative- -- samples in exponential distribution.+ -- | Estimate distribution from sample. Returns 'Nothing' if there is+ -- not enough data, or if no usable fit results from the method+ -- used, e.g., the estimated distribution parameters would be+ -- invalid or inaccurate. fromSample :: G.Vector v a => v a -> Maybe d -- | Generate variates from continuous distribution using inverse -- transform rule.-genContinuous :: (ContDistr d, PrimMonad m) => d -> Gen (PrimState m) -> m Double+genContinuous :: (ContDistr d, StatefulGen g m) => d -> g -> m Double genContinuous d gen = do- x <- uniform gen+ x <- uniformDouble01M gen return $! quantile d x---- | Backwards compatibility with genContinuous.-genContinous :: (ContDistr d, PrimMonad m) => d -> Gen (PrimState m) -> m Double-genContinous = genContinuous-{-# DEPRECATED genContinous "Use genContinuous" #-} data P = P {-# UNPACK #-} !Double {-# UNPACK #-} !Double
Statistics/Distribution/Binomial.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PatternGuards #-} {-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-} -- | -- Module : Statistics.Distribution.Binomial@@ -31,7 +32,7 @@ import Data.Data (Data, Typeable) import GHC.Generics (Generic) import Numeric.SpecFunctions (choose,logChoose,incompleteBeta,log1p)-import Numeric.MathFunctions.Constants (m_epsilon)+import Numeric.MathFunctions.Constants (m_epsilon,m_tiny) import qualified Statistics.Distribution as D import qualified Statistics.Distribution.Poisson.Internal as I@@ -70,6 +71,7 @@ instance D.Distribution BinomialDistribution where cumulative = cumulative+ complCumulative = complCumulative instance D.DiscreteDistr BinomialDistribution where probability = probability@@ -104,9 +106,16 @@ | n == 0 = 1 -- choose could overflow Double for n >= 1030 so we switch to -- log-domain to calculate probability- | n < 1000 = choose n k * p^k * (1-p)^(n-k)- | otherwise = exp $ logChoose n k + log p * k' + log1p (-p) * nk'+ --+ -- We also want to avoid underflow when computing p^k &+ -- (1-p)^(n-k).+ | n < 1000+ , pK >= m_tiny+ , pNK >= m_tiny = choose n k * pK * pNK+ | otherwise = exp $ logChoose n k + log p * k' + log1p (-p) * nk' where+ pK = p^k+ pNK = (1-p)^(n-k) k' = fromIntegral k nk' = fromIntegral $ n - k @@ -119,7 +128,6 @@ k' = fromIntegral k nk' = fromIntegral $ n - k --- Summation from different sides required to reduce roundoff errors cumulative :: BinomialDistribution -> Double -> Double cumulative (BD n p) x | isNaN x = error "Statistics.Distribution.Binomial.cumulative: NaN input"@@ -130,6 +138,16 @@ where k = floor x +complCumulative :: BinomialDistribution -> Double -> Double+complCumulative (BD n p) x+ | isNaN x = error "Statistics.Distribution.Binomial.complCumulative: NaN input"+ | isInfinite x = if x > 0 then 0 else 1+ | k < 0 = 1+ | k >= n = 0+ | otherwise = incompleteBeta (fromIntegral (k+1)) (fromIntegral (n-k)) p+ where+ k = floor x+ mean :: BinomialDistribution -> Double mean (BD n p) = fromIntegral n * p @@ -157,7 +175,7 @@ -> Maybe BinomialDistribution binomialE n p | n < 0 = Nothing- | p >= 0 || p <= 1 = Just (BD n p)+ | p >= 0 && p <= 1 = Just (BD n p) | otherwise = Nothing errMsg :: Int -> Double -> String
Statistics/Distribution/CauchyLorentz.hs view
@@ -94,23 +94,43 @@ instance D.Distribution CauchyDistribution where- cumulative (CD m s) x = 0.5 + atan( (x - m) / s ) / pi+ cumulative (CD m s) x+ | y < -1 = atan (-1/y) / pi+ | otherwise = 0.5 + atan y / pi+ where+ y = (x - m) / s+ complCumulative (CD m s) x+ | y > 1 = atan (1/y) / pi+ | otherwise = 0.5 - atan y / pi+ where+ y = (x - m) / s instance D.ContDistr CauchyDistribution where density (CD m s) x = (1 / pi) / (s * (1 + y*y)) where y = (x - m) / s quantile (CD m s) p- | p > 0 && p < 1 = m + s * tan( pi * (p - 0.5) )- | p == 0 = -1 / 0- | p == 1 = 1 / 0- | otherwise =- error $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p+ | p == 0 = -1 / 0+ | p == 1 = 1 / 0+ | p == 0.5 = m+ | p < 0 = err+ | p < 0.5 = m - s / tan( pi * p )+ | p < 1 = m + s / tan( pi * (1 - p) )+ | otherwise = err+ where+ err = error+ $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p complQuantile (CD m s) p- | p > 0 && p < 1 = m + s * tan( pi * (0.5 - p) )- | p == 0 = 1 / 0- | p == 1 = -1 / 0- | otherwise =- error $ "Statistics.Distribution.CauchyLorentz.complQuantile: p must be in [0,1] range. Got: "++show p+ | p == 0 = 1 / 0+ | p == 1 = -1 / 0+ | p == 0.5 = m+ | p < 0 = err+ | p < 0.5 = m + s / tan( pi * p )+ | p < 1 = m - s / tan( pi * (1 - p) )+ | otherwise = err+ where+ err = error+ $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p+ instance D.ContGen CauchyDistribution where genContVar = D.genContinuous
Statistics/Distribution/DiscreteUniform.hs view
@@ -28,10 +28,11 @@ , rangeTo ) where -import Control.Applicative ((<$>), (<*>), empty)+import Control.Applicative (empty) import Data.Aeson (FromJSON(..), ToJSON, Value(..), (.:)) import Data.Binary (Binary(..)) import Data.Data (Data, Typeable)+import System.Random.Stateful (uniformRM) import GHC.Generics (Generic) import qualified Statistics.Distribution as D@@ -93,6 +94,12 @@ instance D.MaybeEntropy DiscreteUniform where maybeEntropy = Just . D.entropy++instance D.ContGen DiscreteUniform where+ genContVar d = fmap fromIntegral . D.genDiscreteVar d++instance D.DiscreteGen DiscreteUniform where+ genDiscreteVar (U a b) = uniformRM (a,b) -- | Construct discrete uniform distribution on support {1, ..., n}. -- Range /n/ must be >0.
Statistics/Distribution/Exponential.hs view
@@ -11,7 +11,7 @@ -- Portability : portable -- -- The exponential distribution. This is the continuous probability--- distribution of the times between events in a poisson process, in+-- distribution of the times between events in a Poisson process, in -- which events occur continuously and independently at a constant -- average rate. @@ -33,7 +33,6 @@ import Numeric.SpecFunctions (log1p,expm1) import Numeric.MathFunctions.Constants (m_neg_inf) import qualified System.Random.MWC.Distributions as MWC-import qualified Data.Vector.Generic as G import qualified Statistics.Distribution as D import qualified Statistics.Sample as S@@ -136,11 +135,9 @@ errMsg :: Double -> String errMsg l = "Statistics.Distribution.Exponential.exponential: scale parameter must be positive. Got " ++ show l --- | Create exponential distribution from sample. Returns @Nothing@ if--- sample is empty or contains negative elements. No other tests are--- made to check whether it truly is exponential.+-- | Create exponential distribution from sample. Estimates the rate+-- with the maximum likelihood estimator, which is biased. Returns+-- @Nothing@ if the sample mean does not exist or is not positive. instance D.FromSample ExponentialDistribution Double where- fromSample xs- | G.null xs = Nothing- | G.all (>= 0) xs = Nothing- | otherwise = Just $! ED (S.mean xs)+ fromSample xs = let m = S.mean xs+ in if m > 0 then Just (ED (1/m)) else Nothing
Statistics/Distribution/FDistribution.hs view
@@ -109,15 +109,36 @@ cumulative :: FDistribution -> Double -> Double cumulative (F n m _) x | x <= 0 = 0- | isInfinite x = 1 -- Only matches +∞- | otherwise = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y))+ -- Only matches +∞+ | isInfinite x = 1+ -- NOTE: Here we rely on implementation detail of incompleteBeta. It+ -- computes using series expansion for sufficiently small x+ -- and uses following identity otherwise:+ --+ -- I(x; a, b) = 1 - I(1-x; b, a)+ --+ -- Point is we can compute 1-x as m/(m+y) without loss of+ -- precision for large x. Sadly this switchover point is+ -- implementation detail.+ | n >= (n+m)*bx = incompleteBeta (0.5 * n) (0.5 * m) bx+ | otherwise = 1 - incompleteBeta (0.5 * m) (0.5 * n) bx1+ where+ y = n * x+ bx = y / (m + y)+ bx1 = m / (m + y) complCumulative :: FDistribution -> Double -> Double complCumulative (F n m _) x- | x <= 0 = 1- | isInfinite x = 0 -- Only matches +∞- | otherwise = let y = n*x- in incompleteBeta (0.5 * m) (0.5 * n) (m / (m + y))+ | x <= 0 = 1+ -- Only matches +∞+ | isInfinite x = 0+ -- See NOTE at cumulative+ | m >= (n+m)*bx = incompleteBeta (0.5 * m) (0.5 * n) bx+ | otherwise = 1 - incompleteBeta (0.5 * n) (0.5 * m) bx1+ where+ y = n*x+ bx = m / (m + y)+ bx1 = y / (m + y) logDensity :: FDistribution -> Double -> Double logDensity (F n m fac) x
Statistics/Distribution/Gamma.hs view
@@ -36,6 +36,7 @@ import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf) import Numeric.SpecFunctions (incompleteGamma, invIncompleteGamma, logGamma, digamma) import qualified System.Random.MWC.Distributions as MWC+import qualified Numeric.Sum as Sum import Statistics.Distribution.Poisson.Internal as Poisson import qualified Statistics.Distribution as D@@ -126,7 +127,11 @@ density = density logDensity (GD k theta) x | x <= 0 = m_neg_inf- | otherwise = log x * (k - 1) - (x / theta) - logGamma k - log theta * k+ | otherwise = Sum.sum Sum.kbn [ log x * (k - 1)+ , - (x / theta)+ , - logGamma k+ , - log theta * k+ ] quantile = quantile instance D.Variance GammaDistribution where
Statistics/Distribution/Geometric.hs view
@@ -39,7 +39,7 @@ import Data.Binary (Binary(..)) import Data.Data (Data, Typeable) import GHC.Generics (Generic)-import Numeric.MathFunctions.Constants (m_pos_inf, m_neg_inf)+import Numeric.MathFunctions.Constants (m_neg_inf) import Numeric.SpecFunctions (log1p,expm1) import qualified System.Random.MWC.Distributions as MWC @@ -81,10 +81,11 @@ instance D.DiscreteDistr GeometricDistribution where probability (GD s) n | n < 1 = 0- | otherwise = s * (1-s) ** (fromIntegral n - 1)+ | s >= 0.5 = s * (1 - s)^(n - 1)+ | otherwise = s * (exp $ log1p (-s) * (fromIntegral n - 1)) logProbability (GD s) n | n < 1 = m_neg_inf- | otherwise = log s + log (1-s) * (fromIntegral n - 1)+ | otherwise = log s + log1p (-s) * (fromIntegral n - 1) instance D.Mean GeometricDistribution where@@ -102,9 +103,8 @@ instance D.Entropy GeometricDistribution where entropy (GD s)- | s == 0 = m_pos_inf | s == 1 = 0- | otherwise = negate $ (s * log s + (1-s) * log (1-s)) / s+ | otherwise = -(s * log s + (1-s) * log1p (-s)) / s instance D.MaybeEntropy GeometricDistribution where maybeEntropy = Just . D.entropy@@ -120,14 +120,18 @@ | x < 1 = 0 | isInfinite x = 1 | isNaN x = error "Statistics.Distribution.Geometric.cumulative: NaN input"- | otherwise = negate $ expm1 $ fromIntegral (floor x :: Int) * log1p (-s)+ | s >= 0.5 = 1 - (1 - s)^k+ | otherwise = negate $ expm1 $ fromIntegral k * log1p (-s)+ where k = floor x :: Int complCumulative :: GeometricDistribution -> Double -> Double complCumulative (GD s) x | x < 1 = 1 | isInfinite x = 0- | isNaN x = error "Statistics.Distribution.Geometric.cumulative: NaN input"- | otherwise = exp $ fromIntegral (floor x :: Int) * log1p (-s)+ | isNaN x = error "Statistics.Distribution.Geometric.complCumulative: NaN input"+ | s >= 0.5 = (1 - s)^k+ | otherwise = exp $ fromIntegral k * log1p (-s)+ where k = floor x :: Int -- | Create geometric distribution.@@ -139,11 +143,11 @@ geometricE :: Double -- ^ Success rate -> Maybe GeometricDistribution geometricE x- | x >= 0 && x <= 1 = Just (GD x)+ | x > 0 && x <= 1 = Just (GD x) | otherwise = Nothing errMsg :: Double -> String-errMsg x = "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x+errMsg x = "Statistics.Distribution.Geometric.geometric: probability must be in (0,1] range. Got " ++ show x ----------------------------------------------------------------@@ -215,8 +219,8 @@ geometric0E :: Double -- ^ Success rate -> Maybe GeometricDistribution0 geometric0E x- | x >= 0 && x <= 1 = Just (GD0 x)+ | x > 0 && x <= 1 = Just (GD0 x) | otherwise = Nothing errMsg0 :: Double -> String-errMsg0 x = "Statistics.Distribution.Geometric.geometric0: probability must be in [0,1] range. Got " ++ show x+errMsg0 x = "Statistics.Distribution.Geometric.geometric0: probability must be in (0,1] range. Got " ++ show x
Statistics/Distribution/Hypergeometric.hs view
@@ -71,6 +71,7 @@ instance D.Distribution HypergeometricDistribution where cumulative = cumulative+ complCumulative = complCumulative instance D.DiscreteDistr HypergeometricDistribution where probability = probability@@ -133,10 +134,10 @@ errMsg :: Int -> Int -> Int -> String errMsg m l k- = "Statistics.Distribution.Hypergeometric.hypergeometric: "- ++ "m=" ++ show m- ++ "l=" ++ show l- ++ "k=" ++ show k+ = "Statistics.Distribution.Hypergeometric.hypergeometric:"+ ++ " m=" ++ show m+ ++ " l=" ++ show l+ ++ " k=" ++ show k ++ " should hold: l>0 & m in [0,l] & k in (0,l]" -- Naive implementation@@ -164,6 +165,18 @@ | n < minN = 0 | n >= maxN = 1 | otherwise = D.sumProbabilities d minN n+ where+ n = floor x+ minN = max 0 (mi+ki-li)+ maxN = min mi ki++complCumulative :: HypergeometricDistribution -> Double -> Double+complCumulative d@(HD mi li ki) x+ | isNaN x = error "Statistics.Distribution.Hypergeometric.complCumulative: NaN argument"+ | isInfinite x = if x > 0 then 0 else 1+ | n < minN = 1+ | n >= maxN = 0+ | otherwise = D.sumProbabilities d (n + 1) maxN where n = floor x minN = max 0 (mi+ki-li)
Statistics/Distribution/Laplace.hs view
@@ -151,9 +151,9 @@ errMsg _ s = "Statistics.Distribution.Laplace.laplace: scale parameter must be positive. Got " ++ show s --- | Create Laplace distribution from sample. No tests are made to--- check whether it truly is Laplace. Location of distribution--- estimated as median of sample.+-- | Create Laplace distribution from sample. The location is estimated+-- as the median of the sample, and the scale as the mean absolute+-- deviation of the median. instance D.FromSample LaplaceDistribution Double where fromSample xs | G.null xs = Nothing
+ Statistics/Distribution/Lognormal.hs view
@@ -0,0 +1,172 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}+-- |+-- Module : Statistics.Distribution.Lognormal+-- Copyright : (c) 2020 Ximin Luo+-- License : BSD3+--+-- Maintainer : infinity0@pwned.gg+-- Stability : experimental+-- Portability : portable+--+-- The log normal distribution. This is a continuous probability+-- distribution that describes data whose log is clustered around a+-- mean. For example, the multiplicative product of many independent+-- positive random variables.++module Statistics.Distribution.Lognormal+ (+ LognormalDistribution+ -- * Constructors+ , lognormalDistr+ , lognormalDistrErr+ , lognormalDistrMeanStddevErr+ , lognormalStandard+ ) where++import Data.Aeson (FromJSON, ToJSON)+import Data.Binary (Binary (..))+import Data.Data (Data, Typeable)+import GHC.Generics (Generic)+import Numeric.MathFunctions.Constants (m_huge, m_sqrt_2_pi)+import Numeric.SpecFunctions (expm1, log1p)+import qualified Data.Vector.Generic as G++import qualified Statistics.Distribution as D+import qualified Statistics.Distribution.Normal as N+import Statistics.Internal+++-- | The lognormal distribution.+newtype LognormalDistribution = LND N.NormalDistribution+ deriving (Eq, Typeable, Data, Generic)++instance Show LognormalDistribution where+ showsPrec i (LND d) = defaultShow2 "lognormalDistr" m s i+ where+ m = D.mean d+ s = D.stdDev d+instance Read LognormalDistribution where+ readPrec = defaultReadPrecM2 "lognormalDistr" $+ (either (const Nothing) Just .) . lognormalDistrErr++instance ToJSON LognormalDistribution+instance FromJSON LognormalDistribution++instance Binary LognormalDistribution where+ put (LND d) = put m >> put s+ where+ m = D.mean d+ s = D.stdDev d+ get = do+ m <- get+ sd <- get+ either fail return $ lognormalDistrErr m sd++instance D.Distribution LognormalDistribution where+ cumulative = cumulative+ complCumulative = complCumulative++instance D.ContDistr LognormalDistribution where+ logDensity = logDensity+ quantile = quantile+ complQuantile = complQuantile++instance D.MaybeMean LognormalDistribution where+ maybeMean = Just . D.mean++instance D.Mean LognormalDistribution where+ mean (LND d) = exp (m + v / 2)+ where+ m = D.mean d+ v = D.variance d++instance D.MaybeVariance LognormalDistribution where+ maybeStdDev = Just . D.stdDev+ maybeVariance = Just . D.variance++instance D.Variance LognormalDistribution where+ variance (LND d) = expm1 v * exp (2 * m + v)+ where+ m = D.mean d+ v = D.variance d++instance D.Entropy LognormalDistribution where+ entropy (LND d) = logBase 2 (s * exp (m + 0.5) * m_sqrt_2_pi)+ where+ m = D.mean d+ s = D.stdDev d++instance D.MaybeEntropy LognormalDistribution where+ maybeEntropy = Just . D.entropy++instance D.ContGen LognormalDistribution where+ genContVar d = D.genContinuous d++-- | Standard log normal distribution with mu 0 and sigma 1.+--+-- Mean is @sqrt e@ and variance is @(e - 1) * e@.+lognormalStandard :: LognormalDistribution+lognormalStandard = LND N.standard++-- | Create log normal distribution from parameters.+lognormalDistr+ :: Double -- ^ Mu+ -> Double -- ^ Sigma+ -> LognormalDistribution+lognormalDistr mu sig = either error id $ lognormalDistrErr mu sig++-- | Create log normal distribution from parameters.+lognormalDistrErr+ :: Double -- ^ Mu+ -> Double -- ^ Sigma+ -> Either String LognormalDistribution+lognormalDistrErr mu sig+ | sig >= sqrt (log m_huge - 2 * mu) = Left $ errMsg mu sig+ | otherwise = LND <$> N.normalDistrErr mu sig++errMsg :: Double -> Double -> String+errMsg mu sig =+ "Statistics.Distribution.Lognormal.lognormalDistr: sigma must be > 0 && < "+ ++ show lim ++ ". Got " ++ show sig+ where lim = sqrt (log m_huge - 2 * mu)++-- | Create log normal distribution from mean and standard deviation.+lognormalDistrMeanStddevErr+ :: Double -- ^ Mu+ -> Double -- ^ Sigma+ -> Either String LognormalDistribution+lognormalDistrMeanStddevErr m sd = LND <$> N.normalDistrErr mu sig+ where r = sd / m+ sig2 = log1p (r * r)+ sig = sqrt sig2+ mu = log m - sig2 / 2++-- | Variance is estimated using maximum likelihood method+-- (biased estimation) over the log of the data.+--+-- Returns @Nothing@ if sample contains less than one element or+-- variance is zero (all elements are equal)+instance D.FromSample LognormalDistribution Double where+ fromSample = fmap LND . D.fromSample . G.map log++logDensity :: LognormalDistribution -> Double -> Double+logDensity (LND d) x+ | x > 0 = let lx = log x in D.logDensity d lx - lx+ | otherwise = 0++cumulative :: LognormalDistribution -> Double -> Double+cumulative (LND d) x+ | x > 0 = D.cumulative d $ log x+ | otherwise = 0++complCumulative :: LognormalDistribution -> Double -> Double+complCumulative (LND d) x+ | x > 0 = D.complCumulative d $ log x+ | otherwise = 1++quantile :: LognormalDistribution -> Double -> Double+quantile (LND d) = exp . D.quantile d++complQuantile :: LognormalDistribution -> Double -> Double+complQuantile (LND d) = exp . D.complQuantile d
+ Statistics/Distribution/NegativeBinomial.hs view
@@ -0,0 +1,188 @@+{-# LANGUAGE OverloadedStrings, PatternGuards,+ DeriveDataTypeable, DeriveGeneric #-}+-- |+-- Module : Statistics.Distribution.NegativeBinomial+-- Copyright : (c) 2022 Lorenz Minder+-- License : BSD3+--+-- Maintainer : lminder@gmx.net+-- Stability : experimental+-- Portability : portable+--+-- The negative binomial distribution. This is the discrete probability+-- distribution of the number of failures in a sequence of independent+-- yes\/no experiments before a specified number of successes /r/. Each+-- Bernoulli trial has success probability /p/ in the range (0, 1]. The+-- parameter /r/ must be positive, but does not have to be integer.++module Statistics.Distribution.NegativeBinomial (+ NegativeBinomialDistribution+ -- * Constructors+ , negativeBinomial+ , negativeBinomialE+ -- * Accessors+ , nbdSuccesses+ , nbdProbability+) where++import Control.Applicative+import Data.Aeson (FromJSON(..), ToJSON, Value(..), (.:))+import Data.Binary (Binary(..))+import Data.Data (Data, Typeable)+import Data.Foldable (foldl')+import GHC.Generics (Generic)+import Numeric.SpecFunctions (incompleteBeta, log1p)+import Numeric.SpecFunctions.Extra (logChooseFast)+import Numeric.MathFunctions.Constants (m_epsilon, m_tiny)++import qualified Statistics.Distribution as D+import Statistics.Internal++-- Math helper functions++-- | Generalized binomial coefficients.+--+-- These computes binomial coefficients with the small generalization+-- that the /n/ need not be integer, but can be real.+gChoose :: Double -> Int -> Double+gChoose n k+ | k < 0 = 0+ | k' >= 50 = exp $ logChooseFast n k'+ | otherwise = foldl' (*) 1 factors+ where factors = [ (n - k' + j) / j | j <- [1..k'] ]+ k' = fromIntegral k+++-- Implementation of Negative Binomial++-- | The negative binomial distribution.+data NegativeBinomialDistribution = NBD {+ nbdSuccesses :: {-# UNPACK #-} !Double+ -- ^ Number of successes until stop+ , nbdProbability :: {-# UNPACK #-} !Double+ -- ^ Success probability.+ } deriving (Eq, Typeable, Data, Generic)++instance Show NegativeBinomialDistribution where+ showsPrec i (NBD r p) = defaultShow2 "negativeBinomial" r p i+instance Read NegativeBinomialDistribution where+ readPrec = defaultReadPrecM2 "negativeBinomial" negativeBinomialE++instance ToJSON NegativeBinomialDistribution+instance FromJSON NegativeBinomialDistribution where+ parseJSON (Object v) = do+ r <- v .: "nbdSuccesses"+ p <- v .: "nbdProbability"+ maybe (fail $ errMsg r p) return $ negativeBinomialE r p+ parseJSON _ = empty++instance Binary NegativeBinomialDistribution where+ put (NBD r p) = put r >> put p+ get = do+ r <- get+ p <- get+ maybe (fail $ errMsg r p) return $ negativeBinomialE r p++instance D.Distribution NegativeBinomialDistribution where+ cumulative = cumulative+ complCumulative = complCumulative++instance D.DiscreteDistr NegativeBinomialDistribution where+ probability = probability+ logProbability = logProbability++instance D.Mean NegativeBinomialDistribution where+ mean = mean++instance D.Variance NegativeBinomialDistribution where+ variance = variance++instance D.MaybeMean NegativeBinomialDistribution where+ maybeMean = Just . D.mean++instance D.MaybeVariance NegativeBinomialDistribution where+ maybeStdDev = Just . D.stdDev+ maybeVariance = Just . D.variance++instance D.Entropy NegativeBinomialDistribution where+ entropy = directEntropy++instance D.MaybeEntropy NegativeBinomialDistribution where+ maybeEntropy = Just . D.entropy++-- This could be slow for big n+probability :: NegativeBinomialDistribution -> Int -> Double+probability d@(NBD r p) k+ | k < 0 = 0+ -- Switch to log domain for large k + r to avoid overflows.+ --+ -- We also want to avoid underflow when computing (1-p)^k &+ -- p^r.+ | k' + r < 1000+ , pK >= m_tiny+ , pR >= m_tiny = gChoose (k' + r - 1) k * pK * pR+ | otherwise = exp $ logProbability d k+ where+ pK = exp $ log1p (-p) * k'+ pR = p**r+ k' = fromIntegral k++logProbability :: NegativeBinomialDistribution -> Int -> Double+logProbability (NBD r p) k+ | k < 0 = (-1)/0+ | otherwise = logChooseFast (k' + r - 1) k'+ + log1p (-p) * k'+ + log p * r+ where k' = fromIntegral k++cumulative :: NegativeBinomialDistribution -> Double -> Double+cumulative (NBD r p) x+ | isNaN x = error "Statistics.Distribution.NegativeBinomial.cumulative: NaN input"+ | isInfinite x = if x > 0 then 1 else 0+ | k < 0 = 0+ | otherwise = incompleteBeta r (fromIntegral (k+1)) p+ where+ k = floor x :: Integer++complCumulative :: NegativeBinomialDistribution -> Double -> Double+complCumulative (NBD r p) x+ | isNaN x = error "Statistics.Distribution.NegativeBinomial.complCumulative: NaN input"+ | isInfinite x = if x > 0 then 0 else 1+ | k < 0 = 1+ | otherwise = incompleteBeta (fromIntegral (k+1)) r (1 - p)+ where+ k = floor x :: Integer++mean :: NegativeBinomialDistribution -> Double+mean (NBD r p) = r * (1 - p)/p++variance :: NegativeBinomialDistribution -> Double+variance (NBD r p) = r * (1 - p)/(p * p)++directEntropy :: NegativeBinomialDistribution -> Double+directEntropy d =+ negate . sum $+ takeWhile (< -m_epsilon) $+ dropWhile (>= -m_epsilon) $+ [ let x = probability d k in x * log x | k <- [0..]]++-- | Construct negative binomial distribution. Number of successes /r/+-- must be positive and probability must be in (0,1] range+negativeBinomial :: Double -- ^ Number of successes.+ -> Double -- ^ Success probability.+ -> NegativeBinomialDistribution+negativeBinomial r p = maybe (error $ errMsg r p) id $ negativeBinomialE r p++-- | Construct negative binomial distribution. Number of successes /r/+-- must be positive and probability must be in (0,1] range+negativeBinomialE :: Double -- ^ Number of successes.+ -> Double -- ^ Success probability.+ -> Maybe NegativeBinomialDistribution+negativeBinomialE r p+ | r > 0 && 0 < p && p <= 1 = Just (NBD r p)+ | otherwise = Nothing++errMsg :: Double -> Double -> String+errMsg r p+ = "Statistics.Distribution.NegativeBinomial.negativeBinomial: r=" ++ show r+ ++ " p=" ++ show p ++ ", but need r>0 and p in (0,1]"
Statistics/Distribution/Normal.hs view
@@ -19,6 +19,7 @@ -- * Constructors , normalDistr , normalDistrE+ , normalDistrErr , standard ) where @@ -55,7 +56,7 @@ parseJSON (Object v) = do m <- v .: "mean" sd <- v .: "stdDev"- maybe (fail $ errMsg m sd) return $ normalDistrE m sd+ either fail return $ normalDistrErr m sd parseJSON _ = empty instance Binary NormalDistribution where@@ -63,7 +64,7 @@ get = do m <- get sd <- get- maybe (fail $ errMsg m sd) return $ normalDistrE m sd+ either fail return $ normalDistrErr m sd instance D.Distribution NormalDistribution where cumulative = cumulative@@ -111,7 +112,7 @@ normalDistr :: Double -- ^ Mean of distribution -> Double -- ^ Standard deviation of distribution -> NormalDistribution-normalDistr m sd = maybe (error $ errMsg m sd) id $ normalDistrE m sd+normalDistr m sd = either error id $ normalDistrErr m sd -- | Create normal distribution from parameters. --@@ -120,13 +121,20 @@ normalDistrE :: Double -- ^ Mean of distribution -> Double -- ^ Standard deviation of distribution -> Maybe NormalDistribution-normalDistrE m sd- | sd > 0 = Just ND { mean = m- , stdDev = sd- , ndPdfDenom = log $ m_sqrt_2_pi * sd- , ndCdfDenom = m_sqrt_2 * sd- }- | otherwise = Nothing+normalDistrE m sd = either (const Nothing) Just $ normalDistrErr m sd++-- | Create normal distribution from parameters.+--+normalDistrErr :: Double -- ^ Mean of distribution+ -> Double -- ^ Standard deviation of distribution+ -> Either String NormalDistribution+normalDistrErr m sd+ | sd > 0 = Right $ ND { mean = m+ , stdDev = sd+ , ndPdfDenom = log $ m_sqrt_2_pi * sd+ , ndCdfDenom = m_sqrt_2 * sd+ }+ | otherwise = Left $ errMsg m sd errMsg :: Double -> Double -> String errMsg _ sd = "Statistics.Distribution.Normal.normalDistr: standard deviation must be positive. Got " ++ show sd
Statistics/Distribution/Poisson.hs view
@@ -31,15 +31,18 @@ import Data.Binary (Binary(..)) import Data.Data (Data, Typeable) import GHC.Generics (Generic)++import qualified System.Random.MWC.Distributions as MWC+ import Numeric.SpecFunctions (incompleteGamma,logFactorial) import Numeric.MathFunctions.Constants (m_neg_inf) + import qualified Statistics.Distribution as D import qualified Statistics.Distribution.Poisson.Internal as I import Statistics.Internal - newtype PoissonDistribution = PD { poissonLambda :: Double } deriving (Eq, Typeable, Data, Generic)@@ -92,6 +95,14 @@ instance D.MaybeEntropy PoissonDistribution where maybeEntropy = Just . D.entropy++-- | @since 0.16.5.0+instance D.DiscreteGen PoissonDistribution where+ genDiscreteVar (PD lambda) = MWC.poisson lambda++-- | @since 0.16.5.0+instance D.ContGen PoissonDistribution where+ genContVar (PD lambda) gen = fromIntegral <$> MWC.poisson lambda gen -- | Create Poisson distribution. poisson :: Double -> PoissonDistribution
Statistics/Distribution/Poisson/Internal.hs view
@@ -33,7 +33,7 @@ (m_sqrt_2_pi * sqrt x) -- -- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy--- -- of the Poisson Law". This function is unused because 'directEntorpy'+-- -- of the Poisson Law". This function is unused because 'directEntropy' -- -- is just as accurate and is faster by about a factor of 4. -- alyThm1 :: Double -> Double -- alyThm1 lambda =@@ -61,7 +61,7 @@ 1.4189385332046727 + 0.5 * log lambda + zipCoefficients lambda coefficients --- | Returns the average of the upper and lower bounds accounding to+-- | Returns the average of the upper and lower bounds according to -- theorem 2. alyThm2 :: Double -> [Double] -> [Double] -> Double alyThm2 lambda upper lower =@@ -164,7 +164,7 @@ dropWhile (not . (< negate m_epsilon * lambda)) $ [ let x = probability lambda k in x * log x | k <- [0..]] --- | Compute the entropy of a poisson distribution using the best available+-- | Compute the entropy of a Poisson distribution using the best available -- method. poissonEntropy :: Double -> Double poissonEntropy lambda
Statistics/Distribution/StudentT.hs view
@@ -26,7 +26,7 @@ import Data.Data (Data, Typeable) import GHC.Generics (Generic) import Numeric.SpecFunctions (- logBeta, incompleteBeta, invIncompleteBeta, digamma)+ logBeta, incompleteBeta, invIncompleteBeta, digamma, log1p) import qualified Statistics.Distribution as D import Statistics.Distribution.Transform (LinearTransform (..))@@ -94,8 +94,9 @@ logDensityUnscaled :: StudentT -> Double -> Double-logDensityUnscaled (StudentT ndf) x =- log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)+logDensityUnscaled (StudentT ndf) x+ = log1p (x*x/ndf) * (-(0.5 * (1 + ndf)))+ - logBeta 0.5 (0.5 * ndf) quantile :: StudentT -> Double -> Double quantile (StudentT ndf) p
Statistics/Distribution/Transform.hs view
@@ -18,11 +18,9 @@ ) where import Data.Aeson (FromJSON, ToJSON)-import Control.Applicative ((<*>)) import Data.Binary (Binary) import Data.Binary (put, get) import Data.Data (Data, Typeable)-import Data.Functor ((<$>)) import GHC.Generics (Generic) import qualified Statistics.Distribution as D
Statistics/Distribution/Uniform.hs view
@@ -22,11 +22,11 @@ ) where import Control.Applicative-import Data.Aeson (FromJSON(..), ToJSON, Value(..), (.:))-import Data.Binary (Binary(..))-import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified System.Random.MWC as MWC+import Data.Aeson (FromJSON(..), ToJSON, Value(..), (.:))+import Data.Binary (Binary(..))+import Data.Data (Data, Typeable)+import System.Random.Stateful (uniformRM)+import GHC.Generics (Generic) import qualified Statistics.Distribution as D import Statistics.Internal@@ -117,4 +117,4 @@ maybeEntropy = Just . D.entropy instance D.ContGen UniformDistribution where- genContVar (UniformDistribution a b) gen = MWC.uniformR (a,b) gen+ genContVar (UniformDistribution a b) = uniformRM (a,b)
+ Statistics/Distribution/Weibull.hs view
@@ -0,0 +1,224 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}+-- |+-- Module : Statistics.Distribution.Lognormal+-- Copyright : (c) 2020 Ximin Luo+-- License : BSD3+--+-- Maintainer : infinity0@pwned.gg+-- Stability : experimental+-- Portability : portable+--+-- The Weibull distribution. This is a continuous probability+-- distribution that describes the occurrence of a single event whose+-- probability changes over time, controlled by the shape parameter.++module Statistics.Distribution.Weibull+ (+ WeibullDistribution+ -- * Constructors+ , weibullDistr+ , weibullDistrErr+ , weibullStandard+ , weibullDistrApproxMeanStddevErr+ ) where++import Control.Applicative+import Data.Aeson (FromJSON(..), ToJSON, Value(..), (.:))+import Data.Binary (Binary(..))+import Data.Data (Data, Typeable)+import GHC.Generics (Generic)+import Numeric.MathFunctions.Constants (m_eulerMascheroni)+import Numeric.SpecFunctions (expm1, log1p, logGamma)+import qualified Data.Vector.Generic as G++import qualified Statistics.Distribution as D+import qualified Statistics.Sample as S+import Statistics.Internal+++-- | The Weibull distribution.+data WeibullDistribution = WD {+ wdShape :: {-# UNPACK #-} !Double+ , wdLambda :: {-# UNPACK #-} !Double+ } deriving (Eq, Typeable, Data, Generic)++instance Show WeibullDistribution where+ showsPrec i (WD k l) = defaultShow2 "weibullDistr" k l i+instance Read WeibullDistribution where+ readPrec = defaultReadPrecM2 "weibullDistr" $+ (either (const Nothing) Just .) . weibullDistrErr++instance ToJSON WeibullDistribution+instance FromJSON WeibullDistribution where+ parseJSON (Object v) = do+ k <- v .: "wdShape"+ l <- v .: "wdLambda"+ either fail return $ weibullDistrErr k l+ parseJSON _ = empty++instance Binary WeibullDistribution where+ put (WD k l) = put k >> put l+ get = do+ k <- get+ l <- get+ either fail return $ weibullDistrErr k l++instance D.Distribution WeibullDistribution where+ cumulative = cumulative+ complCumulative = complCumulative++instance D.ContDistr WeibullDistribution where+ logDensity = logDensity+ quantile = quantile+ complQuantile = complQuantile++instance D.MaybeMean WeibullDistribution where+ maybeMean = Just . D.mean++instance D.Mean WeibullDistribution where+ mean (WD k l) = l * exp (logGamma (1 + 1 / k))++instance D.MaybeVariance WeibullDistribution where+ maybeStdDev = Just . D.stdDev+ maybeVariance = Just . D.variance++instance D.Variance WeibullDistribution where+ variance (WD k l) = l * l * (exp (logGamma (1 + 2 * invk)) - q * q)+ where+ invk = 1 / k+ q = exp (logGamma (1 + invk))++instance D.Entropy WeibullDistribution where+ entropy (WD k l) = m_eulerMascheroni * (1 - 1 / k) + log (l / k) + 1++instance D.MaybeEntropy WeibullDistribution where+ maybeEntropy = Just . D.entropy++instance D.ContGen WeibullDistribution where+ genContVar d = D.genContinuous d++-- | Standard Weibull distribution with scale factor (lambda) 1.+weibullStandard :: Double -> WeibullDistribution+weibullStandard k = weibullDistr k 1.0++-- | Create Weibull distribution from parameters.+--+-- If the shape (first) parameter is @1.0@, the distribution is equivalent to a+-- 'Statistics.Distribution.Exponential.ExponentialDistribution' with parameter+-- @1 / lambda@ the scale (second) parameter.+weibullDistr+ :: Double -- ^ Shape+ -> Double -- ^ Lambda (scale)+ -> WeibullDistribution+weibullDistr k l = either error id $ weibullDistrErr k l++-- | Create Weibull distribution from parameters.+--+-- If the shape (first) parameter is @1.0@, the distribution is equivalent to a+-- 'Statistics.Distribution.Exponential.ExponentialDistribution' with parameter+-- @1 / lambda@ the scale (second) parameter.+weibullDistrErr+ :: Double -- ^ Shape+ -> Double -- ^ Lambda (scale)+ -> Either String WeibullDistribution+weibullDistrErr k l | k <= 0 = Left $ errMsg k l+ | l <= 0 = Left $ errMsg k l+ | otherwise = Right $ WD k l++errMsg :: Double -> Double -> String+errMsg k l =+ "Statistics.Distribution.Weibull.weibullDistr: both shape and lambda must be positive. Got shape "+ ++ show k+ ++ " and lambda "+ ++ show l++-- | Create Weibull distribution from mean and standard deviation.+--+-- The algorithm is from "Methods for Estimating Wind Speed Frequency+-- Distributions", C. G. Justus, W. R. Hargreaves, A. Mikhail, D. Graber, 1977.+-- Given the identity:+--+-- \[+-- (\frac{\sigma}{\mu})^2 = \frac{\Gamma(1+2/k)}{\Gamma(1+1/k)^2} - 1+-- \]+--+-- \(k\) can be approximated by+--+-- \[+-- k \approx (\frac{\sigma}{\mu})^{-1.086}+-- \]+--+-- \(\lambda\) is then calculated straightforwardly via the identity+--+-- \[+-- \lambda = \frac{\mu}{\Gamma(1+1/k)}+-- \]+--+-- Numerically speaking, the approximation for \(k\) is accurate only within a+-- certain range. We arbitrarily pick the range \(0.033 \le \frac{\sigma}{\mu} \le 1.45\)+-- where it is good to ~6%, and will refuse to create a distribution outside of+-- this range. The paper does not cover these details but it is straightforward+-- to check them numerically.+weibullDistrApproxMeanStddevErr+ :: Double -- ^ Mean+ -> Double -- ^ Stddev+ -> Either String WeibullDistribution+weibullDistrApproxMeanStddevErr m s = if r > 1.45 || r < 0.033+ then Left msg+ else weibullDistrErr k l+ where r = s / m+ k = (s / m) ** (-1.086)+ l = m / exp (logGamma (1 + 1/k))+ msg = "Statistics.Distribution.Weibull.weibullDistr: stddev-mean ratio "+ ++ "outside approximation accuracy range [0.033, 1.45]. Got "+ ++ "stddev " ++ show s ++ " and mean " ++ show m++-- | Uses an approximation based on the mean and standard deviation in+-- 'weibullDistrEstMeanStddevErr', with standard deviation estimated+-- using maximum likelihood method (unbiased estimation).+--+-- Returns @Nothing@ if sample contains less than one element or+-- variance is zero (all elements are equal), or if the estimated mean+-- and standard-deviation lies outside the range for which the+-- approximation is accurate.+instance D.FromSample WeibullDistribution Double where+ fromSample xs+ | G.length xs <= 1 = Nothing+ | v == 0 = Nothing+ | otherwise = either (const Nothing) Just $+ weibullDistrApproxMeanStddevErr m (sqrt v)+ where+ (m,v) = S.meanVarianceUnb xs++logDensity :: WeibullDistribution -> Double -> Double+logDensity (WD k l) x+ | x < 0 = 0+ | otherwise = log k + (k - 1) * log x - k * log l - (x / l) ** k++cumulative :: WeibullDistribution -> Double -> Double+cumulative (WD k l) x | x < 0 = 0+ | otherwise = -expm1 (-(x / l) ** k)++complCumulative :: WeibullDistribution -> Double -> Double+complCumulative (WD k l) x | x < 0 = 1+ | otherwise = exp (-(x / l) ** k)++quantile :: WeibullDistribution -> Double -> Double+quantile (WD k l) p+ | p == 0 = 0+ | p == 1 = inf+ | p > 0 && p < 1 = l * (-log1p (-p)) ** (1 / k)+ | otherwise =+ error $ "Statistics.Distribution.Weibull.quantile: p must be in [0,1] range. Got: " ++ show p+ where inf = 1 / 0++complQuantile :: WeibullDistribution -> Double -> Double+complQuantile (WD k l) q+ | q == 0 = inf+ | q == 1 = 0+ | q > 0 && q < 1 = l * (-log q) ** (1 / k)+ | otherwise =+ error $ "Statistics.Distribution.Weibull.complQuantile: q must be in [0,1] range. Got: " ++ show q+ where inf = 1 / 0
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/Internal.hs view
@@ -25,8 +25,6 @@ import Control.Applicative import Control.Monad import Text.Read-import Data.Orphans ()- ----------------------------------------------------------------
Statistics/Quantile.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE DeriveFoldable #-} {-# LANGUAGE DeriveFunctor #-}@@ -55,7 +54,6 @@ import Data.Aeson (ToJSON,FromJSON) import Data.Data (Data,Typeable) import Data.Default.Class-import Data.Functor import qualified Data.Foldable as F import Data.Vector.Generic ((!)) import qualified Data.Vector as V@@ -196,7 +194,7 @@ | F.any (badQ nQ) qs = modErr "quantiles" "Wrong quantile number" | G.any isNaN xs = modErr "quantiles" "Sample contains NaNs" -- Doesn't matter what we put into empty container- | fnull qs = 0 <$ qs+ | null qs = 0 <$ qs | otherwise = fmap (estimateQuantile sortedXs) ks' where ks' = fmap (\q -> toPk param n q nQ) qs@@ -209,14 +207,6 @@ :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> U.Vector Double -> f Double #-} {-# SPECIALIZE quantiles :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> S.Vector Double -> f Double #-}---- COMPAT-fnull :: F.Foldable f => f a -> Bool-#if !MIN_VERSION_base(4,8,0)-fnull = F.foldr (\_ _ -> False) True-#else-fnull = null-#endif -- | O(/k·n/·log /n/). Same as quantiles but uses 'G.Vector' container -- instead of 'Foldable' one.
Statistics/Regression.hs view
@@ -13,7 +13,6 @@ , bootstrapRegress ) where -import Control.Applicative ((<$>)) import Control.Concurrent.Async (forConcurrently) import Control.DeepSeq (rnf) import Control.Monad (when)@@ -42,8 +41,15 @@ -- element than the list of predictors; the last element is the -- /y/-intercept value. ----- * /R²/, 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@@ -66,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@@ -93,7 +122,7 @@ where n = rows r l = U.length b --- | Compute /R²/, 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@@ -102,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/Resampling.hs view
@@ -1,5 +1,4 @@ {-# LANGUAGE BangPatterns #-}-{-# LANGUAGE CPP #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE DeriveFoldable #-} {-# LANGUAGE DeriveFunctor #-}@@ -40,7 +39,6 @@ ) where import Data.Aeson (FromJSON, ToJSON)-import Control.Applicative import Control.Concurrent.Async (forConcurrently_) import Control.Monad (forM_, forM, replicateM, liftM2) import Control.Monad.Primitive (PrimMonad(..))@@ -88,9 +86,7 @@ , resamples :: v a } deriving (Eq, Read, Show , Generic, Functor, T.Foldable, T.Traversable-#if __GLASGOW_HASKELL__ >= 708 , Typeable, Data-#endif ) instance (Binary a, Binary (v a)) => Binary (Bootstrap v a) where
Statistics/Resampling/Bootstrap.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} -- | -- Module : Statistics.Resampling.Bootstrap -- Copyright : (c) 2009, 2011 Bryan O'Sullivan@@ -31,9 +30,7 @@ import qualified Statistics.Resampling as R -#if !defined(__GHCJS__)-import Control.Monad.Par (parMap, runPar)-#endif+import Control.Parallel.Strategies (parMap, rdeepseq) data T = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double infixl 2 :<@@ -51,15 +48,7 @@ -- this. -> [Estimate ConfInt Double] bootstrapBCA confidenceLevel sample resampledData-#if defined(__GHCJS__)- -- monad-par causes seems to cause "thread blocked indefinitely on MVar"- -- on GHCJS still- --- -- I (phadej) would change the interface to return IO, and use mapConcurrently from async- = map e resampledData-#else- = runPar $ parMap e resampledData-#endif+ = parMap rdeepseq e resampledData where e (est, Bootstrap pt resample) | U.length sample == 1 || isInfinite bias =
Statistics/Sample.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE BangPatterns #-} -- | -- Module : Statistics.Sample -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan@@ -20,6 +21,7 @@ , range -- * Statistics of location+ , expectation , mean , welfordMean , meanWeighted@@ -51,20 +53,23 @@ , fastVarianceUnbiased , fastStdDev - -- * Joint distirbutions+ -- * 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 +81,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 +135,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 +151,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@@ -215,7 +228,7 @@ -- $variance ----- The variance—and hence the standard deviation—of a+-- The variance — and hence the standard deviation — of a -- sample of fewer than two elements are both defined to be zero. -- $robust@@ -354,42 +367,79 @@ -- | 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 -- samples have equal size. pair :: (G.Vector v a, G.Vector v b, G.Vector v (a,b)) => v a -> v b -> v (a,b)@@ -403,8 +453,9 @@ -- (^) operator from Prelude is just slow. (^) :: Double -> Int -> Double-x ^ 1 = x-x ^ n = x * (x ^ (n-1))+x0 ^ n0 = go (n0-1) x0 where+ go 0 !acc = acc+ go n acc = go (n-1) (acc*x0) {-# INLINE (^) #-} -- don't support polymorphism, as we can't get unboxed returns if we use it.
Statistics/Sample/Histogram.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE FlexibleContexts, BangPatterns #-}+{-# LANGUAGE FlexibleContexts, BangPatterns, ScopedTypeVariables #-} -- | -- Module : Statistics.Sample.Histogram@@ -19,6 +19,7 @@ , range ) where +import Control.Monad.ST import Numeric.MathFunctions.Constants (m_epsilon,m_tiny) import Statistics.Function (minMax) import qualified Data.Vector.Generic as G@@ -49,7 +50,7 @@ -- -- Interval (bin) sizes are uniform, based on the supplied upper -- and lower bounds.-histogram_ :: (Num b, RealFrac a, G.Vector v0 a, G.Vector v1 b) =>+histogram_ :: forall b a v0 v1. (Num b, RealFrac a, G.Vector v0 a, G.Vector v1 b) => Int -- ^ Number of bins. This value must be positive. A zero -- or negative value will cause an error.@@ -65,6 +66,7 @@ -> v1 b histogram_ numBins lo hi xs0 = G.create (GM.replicate numBins 0 >>= bin xs0) where+ bin :: forall s. v0 a -> G.Mutable v1 s b -> ST s (G.Mutable v1 s b) bin xs bins = go 0 where go i | i >= len = return bins@@ -75,7 +77,7 @@ go (i+1) write' bins' b !e = GM.write bins' b e len = G.length xs- d = ((hi - lo) * (1 + realToFrac m_epsilon)) / fromIntegral numBins+ d = ((hi - lo) / fromIntegral numBins) * (1 + realToFrac m_epsilon) {-# INLINE histogram_ #-} -- | /O(n)/ Compute decent defaults for the lower and upper bounds of
Statistics/Sample/Internal.hs view
@@ -14,9 +14,10 @@ ( robustSumVar , sum+ , sumF ) where -import Numeric.Sum (kbn, sumVector)+import qualified Numeric.Sum as Sum import Prelude hiding (sum) import Statistics.Function (square) import qualified Data.Vector.Generic as G@@ -26,5 +27,9 @@ {-# INLINE robustSumVar #-} sum :: (G.Vector v Double) => v Double -> Double-sum = sumVector kbn+sum = Sum.sumVector Sum.kbn {-# INLINE sum #-}++sumF :: Foldable f => f Double -> Double+sumF = Sum.sum Sum.kbn+{-# INLINE sumF #-}
+ Statistics/Test/Bartlett.hs view
@@ -0,0 +1,99 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE FlexibleContexts #-}+{-|+Module : Statistics.Test.Bartlett+Description : Bartlett's test for homogeneity of variances.+Copyright : (c) Praneya Kumar, Alexey Khudyakov, 2025+License : BSD-3-Clause++Bartlett's test is used to check that multiple groups of observations+come from distributions with equal variances. This test assumes that+samples come from normal distribution. If this is not the case it may+simple test for non-normality and Levene's ("Statistics.Test.Levene")+is preferred++>>> import qualified Data.Vector.Unboxed as VU+>>> import Statistics.Test.Bartlett+>>> :{+let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]+ b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]+ c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]+in bartlettTest [a,b,c]+:}+Right (Test {testSignificance = mkPValue 1.1254782518843598e-5, testStatistics = 22.789434813726768, testDistribution = chiSquared 2})++-}+module Statistics.Test.Bartlett (+ bartlettTest,+ module Statistics.Distribution.ChiSquared+) where++import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Primitive as VP+#if MIN_VERSION_vector(0,13,2)+import qualified Data.Vector.Strict as VV+#endif++import Statistics.Distribution (complCumulative)+import Statistics.Distribution.ChiSquared (chiSquared, ChiSquared(..))+import Statistics.Sample (varianceUnbiased)+import Statistics.Types (mkPValue)+import Statistics.Test.Types (Test(..))++-- | Perform Bartlett's test for equal variances. The input is a list+-- of vectors, where each vector represents a group of observations.+bartlettTest :: VG.Vector v Double => [v Double] -> Either String (Test ChiSquared)+bartlettTest groups+ | length groups < 2 = Left "At least two groups are required for Bartlett's test."+ | any ((< 2) . VG.length) groups = Left "Each group must have at least two observations."+ | any ((<= 0) . var) groupVariances = Left "All groups must have positive variance."+ | otherwise = Right Test+ { testSignificance = pValue+ , testStatistics = tStatistic+ , testDistribution = chiDist+ }+ where+ -- Number of groups+ k = length groups+ -- Sample sizes for each group+ ni = map (fromIntegral . VG.length) groups+ -- Total number of observations across all groups+ n_tot = sum $ fromIntegral . VG.length <$> groups+ -- Variance estimates+ groupVariances = toVar <$> groups+ sumWeightedVars = sum [ (n - 1) * v | Var{sampleN=n, var=v} <- groupVariances ]+ pooledVariance = sumWeightedVars / fromIntegral (n_tot - k)+ -- Numerator of Bartlett's statistic+ numerator =+ fromIntegral (n_tot - k) * log pooledVariance -+ sum [ (n - 1) * log v | Var{sampleN=n, var=v} <- groupVariances ]+ -- Denominator correction term+ sumReciprocals = sum [1 / (n - 1) | n <- ni]+ denomCorrection =+ 1 + (sumReciprocals - 1 / fromIntegral (n_tot - k)) / (3 * (fromIntegral k - 1))++ -- Test statistic and test distrubution+ tStatistic = max 0 $ numerator / denomCorrection+ chiDist = chiSquared (k - 1)+ pValue = mkPValue $ complCumulative chiDist tStatistic+{-# SPECIALIZE bartlettTest :: [V.Vector Double] -> Either String (Test ChiSquared) #-}+{-# SPECIALIZE bartlettTest :: [VU.Vector Double] -> Either String (Test ChiSquared) #-}+{-# SPECIALIZE bartlettTest :: [VS.Vector Double] -> Either String (Test ChiSquared) #-}+{-# SPECIALIZE bartlettTest :: [VP.Vector Double] -> Either String (Test ChiSquared) #-}+#if MIN_VERSION_vector(0,13,2)+{-# SPECIALIZE bartlettTest :: [VV.Vector Double] -> Either String (Test ChiSquared) #-}+#endif++-- Estimate of variance+data Var = Var+ { sampleN :: !Double -- ^ N of elements+ , var :: !Double -- ^ Sample variance+ }++toVar :: VG.Vector v Double => v Double -> Var+toVar xs = Var { sampleN = fromIntegral $ VG.length xs+ , var = varianceUnbiased xs+ }
Statistics/Test/ChiSquared.hs view
@@ -17,8 +17,8 @@ import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U--+import qualified Data.Vector.Fusion.Bundle as F+import qualified Numeric.Sum as Sum -- | Generic form of Pearson chi squared tests for binned data. Data -- sample is supplied in form of tuples (observed quantity,@@ -26,7 +26,7 @@ -- -- This test should be used only if all bins have expected values of -- at least 5.-chi2test :: (G.Vector v (Int,Double), G.Vector v Double)+chi2test :: (G.Vector v (Int,Double)) => Int -- ^ Number of additional degrees of -- freedom. One degree of freedom -- is due to the fact that the are@@ -39,12 +39,15 @@ | n > 0 = Just Test { testSignificance = mkPValue $ complCumulative d chi2 , testStatistics = chi2- , testDistribution = chiSquared ndf+ , testDistribution = chiSquared n } | otherwise = Nothing where n = G.length vec - ndf - 1- chi2 = sum $ G.map (\(o,e) -> square (fromIntegral o - e) / e) vec+ chi2 = Sum.kbn+ $ F.foldl' Sum.add Sum.zero+ $ F.map (\(o,e) -> square (fromIntegral o - e) / e)+ $ G.stream vec d = chiSquared n {-# INLINABLE chi2test #-} {-# SPECIALIZE@@ -56,7 +59,7 @@ -- | Chi squared test for data with normal errors. Data is supplied in -- form of pair (observation with error, and expectation). chi2testCont- :: (G.Vector v (Estimate NormalErr Double, Double), G.Vector v Double)+ :: (G.Vector v (Estimate NormalErr Double, Double)) => Int -- ^ Number of additional -- degrees of freedom. -> v (Estimate NormalErr Double, Double) -- ^ Observation and expectation.@@ -66,10 +69,13 @@ | n > 0 = Just Test { testSignificance = mkPValue $ complCumulative d chi2 , testStatistics = chi2- , testDistribution = chiSquared ndf+ , testDistribution = chiSquared n } | otherwise = Nothing where n = G.length vec - ndf - 1- chi2 = sum $ G.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s) vec+ chi2 = Sum.kbn+ $ F.foldl' Sum.add Sum.zero+ $ F.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s)+ $ G.stream vec d = chiSquared n
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/KolmogorovSmirnov.hs view
@@ -21,7 +21,7 @@ , kolmogorovSmirnovCdfD , kolmogorovSmirnovD , kolmogorovSmirnov2D- -- * Probablities+ -- * Probabilities , kolmogorovSmirnovProbability -- * References -- $references
Statistics/Test/KruskalWallis.hs view
@@ -17,7 +17,6 @@ ) where import Data.Ord (comparing)-import Data.Foldable (foldMap) import qualified Data.Vector.Unboxed as U import Statistics.Function (sort, sortBy, square) import Statistics.Distribution (complCumulative)
+ Statistics/Test/Levene.hs view
@@ -0,0 +1,153 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE FlexibleContexts #-}+{-|+Module : Statistics.Test.Levene+Description : Levene's test for homogeneity of variances.+Copyright : (c) Praneya Kumar, Alexey Khudyakov, 2025+License : BSD-3-Clause++Levene's test used to check whether samples have equal variance. Null+hypothesis is all samples are from distributions with same variance+(homoscedacity). Test is robust to non-normality, and versatile with+mean or median centering.++>>> import qualified Data.Vector.Unboxed as VU+>>> import Statistics.Test.Levene+>>> :{+let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]+ b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]+ c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]+in levenesTest Median [a, b, c]+:}+Right (Test {testSignificance = mkPValue 2.4315059672496814e-3, testStatistics = 7.584952754501659, testDistribution = fDistributionReal 2.0 27.0})+-}+module Statistics.Test.Levene (+ Center(..),+ levenesTest+) where++import Control.Monad+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Primitive as VP+#if MIN_VERSION_vector(0,13,2)+import qualified Data.Vector.Strict as VV+#endif+import Statistics.Distribution (complCumulative)+import Statistics.Distribution.FDistribution (fDistribution, FDistribution)+import Statistics.Types (mkPValue)+import Statistics.Test.Types (Test(..))+import Statistics.Function (gsort)+import Statistics.Sample (mean)++import qualified Statistics.Sample.Internal as IS+import Statistics.Quantile+++-- | Center calculation method+data Center+ = Mean -- ^ Use arithmetic mean+ | Median -- ^ Use median+ | Trimmed !Double -- ^ Trimmed mean with given proportion to cut from each end+ deriving (Eq, Show)++-- | Main Levene's test function with full error handling+levenesTest+ :: (VG.Vector v Double)+ => Center -- ^ Centering method+ -> [v Double] -- ^ Input samples+ -> Either String (Test FDistribution)+{-# INLINABLE levenesTest #-}+levenesTest center samples+ | length samples < 2 = Left "At least two samples required"+ -- NOTE: We don't have nice way of computing mean of a list!+ | otherwise = do+ let residuals = computeResiduals center <$> samples+ -- Average of all Z+ let n_tot = sum $ VG.length . vecZ <$> residuals -- Total number of samples+ let zbar = IS.sumF [ meanZ z * sampleN z+ | z <- residuals]+ / fromIntegral n_tot+ -- Numerator: Sum over (ni * (Z[i] - Z)^2)+ let numerator = IS.sumF [ sampleN z * sqr (meanZ z - zbar)+ | z <- residuals]+ -- Denominator: Sum over Σ((dev_ij - zbari)^2)+ let denominator = IS.sumF+ [ IS.sum $ VU.map (sqr . subtract (meanZ z)) (vecZ z)+ | z <- residuals+ ]+ -- Handle division by zero and invalid values+ when (denominator <= 0 || isNaN denominator || isInfinite denominator)+ $ Left "Invalid denominator in W-statistic calculation"+ let wStat = (fromIntegral (n_tot - k) / fromIntegral (k - 1)) * (numerator / denominator)+ fDist = fDistribution (k - 1) (n_tot - k)+ Right Test { testStatistics = wStat+ , testSignificance = mkPValue $ complCumulative fDist wStat+ , testDistribution = fDist+ }+ where+ k = length samples -- Number of groups+{-# SPECIALIZE levenesTest :: Center -> [V.Vector Double] -> Either String (Test FDistribution) #-}+{-# SPECIALIZE levenesTest :: Center -> [VU.Vector Double] -> Either String (Test FDistribution) #-}+{-# SPECIALIZE levenesTest :: Center -> [VS.Vector Double] -> Either String (Test FDistribution) #-}+{-# SPECIALIZE levenesTest :: Center -> [VP.Vector Double] -> Either String (Test FDistribution) #-}+#if MIN_VERSION_vector(0,13,2)+{-# SPECIALIZE levenesTest :: Center -> [VV.Vector Double] -> Either String (Test FDistribution) #-}+#endif++----------------------------------------------------------------+-- Implementation+----------------------------------------------------------------++-- | Trim data from both ends with error handling and performance optimization+trimboth :: (Ord a, Fractional a, VG.Vector v a)+ => v a+ -> Double+ -> v a+{-# INLINE trimboth #-}+trimboth vec p+ | p < 0 || p >= 0.5 = error "Statistics.Test.Levene: trimming: proportion must be between 0 and 0.5"+ | VG.null vec = vec+ | otherwise = VG.slice lowerCut (upperCut - lowerCut) sorted+ where+ n = VG.length vec+ sorted = gsort vec+ lowerCut = ceiling $ p * fromIntegral n+ upperCut = n - lowerCut++data Residuals = Residuals+ { sampleN :: !Double+ , meanZ :: !Double+ , vecZ :: !(VU.Vector Double)+ }++computeResiduals+ :: VG.Vector v Double+ => Center+ -> v Double+ -> Residuals+{-# INLINE computeResiduals #-}+computeResiduals method xs = case method of+ Mean ->+ let c = mean xs+ zs = VU.map (\x -> abs (x - c)) $ VU.convert xs+ in makeR zs+ Median ->+ let c = median medianUnbiased xs+ zs = VU.map (\x -> abs (x - c)) $ VU.convert xs+ in makeR zs+ Trimmed p ->+ let trimmed = trimboth xs p+ c = mean trimmed+ zs = VU.map (\x -> abs (x - c)) $ VU.convert trimmed+ in makeR zs+ where+ makeR zs = Residuals { sampleN = fromIntegral $ VU.length zs+ , meanZ = mean zs+ , vecZ = zs+ }++sqr :: Double -> Double+sqr x = x * x
Statistics/Test/MannWhitneyU.hs view
@@ -24,7 +24,6 @@ -- $references ) where -import Control.Applicative ((<$>)) import Data.List (findIndex) import Data.Ord (comparing) import Numeric.SpecFunctions (choose)
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)@@ -134,16 +134,16 @@ -- Calculate T-statistics for paired sample-tStatisticsPaired :: (G.Vector v (Double, Double), G.Vector v Double)+tStatisticsPaired :: (G.Vector v (Double, Double)) => v (Double, Double) -> (Double, Double) {-# INLINE tStatisticsPaired #-} tStatisticsPaired sample = (t, ndf) where -- t-statistics- t = let d = G.map (uncurry (-)) sample- sumd = G.sum d- in sumd / sqrt ((n * G.sum (G.map square d) - square sumd) / ndf)+ t = let d = U.map (uncurry (-)) $ G.convert sample+ sumd = U.sum d+ in sumd / sqrt ((n * U.sum (U.map square d) - square sumd) / ndf) -- degree of freedom ndf = n - 1 n = fromIntegral $ G.length sample
Statistics/Test/WilcoxonT.hs view
@@ -39,7 +39,6 @@ -- the sum of negative ranks (the ranks of the differences where the second parameter is higher). -- to the length of the shorter sample. -import Control.Applicative ((<$>)) import Data.Function (on) import Data.List (findIndex) import Data.Ord (comparing)
Statistics/Types.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE TypeFamilies #-}@@ -72,12 +71,6 @@ import Data.Vector.Unboxed (Unbox) import Data.Vector.Unboxed.Deriving (derivingUnbox) import GHC.Generics (Generic)--#if __GLASGOW_HASKELL__ == 704-import qualified Data.Vector.Generic-import qualified Data.Vector.Generic.Mutable-#endif- import Statistics.Internal import Statistics.Types.Internal import Statistics.Distribution@@ -101,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@.@@ -141,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")@@ -151,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)@@ -162,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 @@ -170,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@@ -306,6 +299,7 @@ -- > , confIntUDX = 6 -- > , confIntCL = cl95 -- > }+-- > } -- -- Prior to statistics 0.14 @Estimate@ data type used following definition: --@@ -324,9 +318,7 @@ , estError :: !(e a) -- ^ Confidence interval for estimate. } deriving (Eq, Read, Show, Generic-#if __GLASGOW_HASKELL__ >= 708 , Typeable, Data-#endif ) instance (Binary (e a), Binary a) => Binary (Estimate e a) where
+ 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/Main.hs view
@@ -0,0 +1,77 @@+module Main where++import Data.Complex+import Statistics.Sample+import Statistics.Transform+import Statistics.Correlation+import System.Random.MWC+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as MVU++import Bench+++-- Test sample+sample :: VU.Vector Double+sample = VU.create $ do g <- create+ MVU.replicateM 10000 (uniform g)++-- Weighted test sample+sampleW :: VU.Vector (Double,Double)+sampleW = VU.zip sample (VU.reverse sample)++-- Complex vector for FFT tests+sampleC :: VU.Vector (Complex Double)+sampleC = VU.zipWith (:+) sample (VU.reverse sample)+++-- Simple benchmark for functions from Statistics.Sample+main :: IO ()+main =+ defaultMain+ [ bgroup "sample"+ [ bench "range" $ nf (\x -> range x) sample+ -- Mean+ , bench "mean" $ nf (\x -> mean x) sample+ , bench "meanWeighted" $ nf (\x -> meanWeighted x) sampleW+ , bench "harmonicMean" $ nf (\x -> harmonicMean x) sample+ , bench "geometricMean" $ nf (\x -> geometricMean x) sample+ -- Variance+ , bench "variance" $ nf (\x -> variance x) sample+ , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample+ , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW+ -- Correlation+ , 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+ , bench "kurtosis" $ nf (\x -> kurtosis x) sample+ -- Central moments+ , bench "C.M. 2" $ nf (\x -> centralMoment 2 x) sample+ , bench "C.M. 3" $ nf (\x -> centralMoment 3 x) sample+ , bench "C.M. 4" $ nf (\x -> centralMoment 4 x) sample+ , bench "C.M. 5" $ nf (\x -> centralMoment 5 x) sample+ ]+ , bgroup "FFT"+ [ bgroup "fft"+ [ bench (show n) $ whnf fft (VU.take n sampleC) | n <- fftSizes ]+ , bgroup "ifft"+ [ bench (show n) $ whnf ifft (VU.take n sampleC) | n <- fftSizes ]+ , bgroup "dct"+ [ bench (show n) $ whnf dct (VU.take n sample) | n <- fftSizes ]+ , bgroup "dct_"+ [ bench (show n) $ whnf dct_ (VU.take n sampleC) | n <- fftSizes ]+ , bgroup "idct"+ [ bench (show n) $ whnf idct (VU.take n sample) | n <- fftSizes ]+ , bgroup "idct_"+ [ bench (show n) $ whnf idct_ (VU.take n sampleC) | n <- fftSizes ]+ ]+ ]+++fftSizes :: [Int]+fftSizes = [32,128,512,2048]
− benchmark/bench.hs
@@ -1,71 +0,0 @@-import Control.Monad.ST (runST)-import Criterion.Main-import Data.Complex-import Statistics.Sample-import Statistics.Transform-import Statistics.Correlation.Pearson-import System.Random.MWC-import qualified Data.Vector.Unboxed as U----- Test sample-sample :: U.Vector Double-sample = runST $ flip uniformVector 10000 =<< create---- Weighted test sample-sampleW :: U.Vector (Double,Double)-sampleW = U.zip sample (U.reverse sample)---- Comlex vector for FFT tests-sampleC :: U.Vector (Complex Double)-sampleC = U.zipWith (:+) sample (U.reverse sample)----- Simple benchmark for functions from Statistics.Sample-main :: IO ()-main =- defaultMain- [ bgroup "sample"- [ bench "range" $ nf (\x -> range x) sample- -- Mean- , bench "mean" $ nf (\x -> mean x) sample- , bench "meanWeighted" $ nf (\x -> meanWeighted x) sampleW- , bench "harmonicMean" $ nf (\x -> harmonicMean x) sample- , bench "geometricMean" $ nf (\x -> geometricMean x) sample- -- Variance- , bench "variance" $ nf (\x -> variance x) sample- , 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- -- Other- , bench "stdDev" $ nf (\x -> stdDev x) sample- , bench "skewness" $ nf (\x -> skewness x) sample- , bench "kurtosis" $ nf (\x -> kurtosis x) sample- -- Central moments- , bench "C.M. 2" $ nf (\x -> centralMoment 2 x) sample- , bench "C.M. 3" $ nf (\x -> centralMoment 3 x) sample- , bench "C.M. 4" $ nf (\x -> centralMoment 4 x) sample- , bench "C.M. 5" $ nf (\x -> centralMoment 5 x) sample- ]- , bgroup "FFT"- [ bgroup "fft"- [ bench (show n) $ whnf fft (U.take n sampleC) | n <- fftSizes ]- , bgroup "ifft"- [ bench (show n) $ whnf ifft (U.take n sampleC) | n <- fftSizes ]- , bgroup "dct"- [ bench (show n) $ whnf dct (U.take n sample) | n <- fftSizes ]- , bgroup "dct_"- [ bench (show n) $ whnf dct_ (U.take n sampleC) | n <- fftSizes ]- , bgroup "idct"- [ bench (show n) $ whnf idct (U.take n sample) | n <- fftSizes ]- , bgroup "idct_"- [ bench (show n) $ whnf idct_ (U.take n sampleC) | n <- fftSizes ]- ]- ]---fftSizes :: [Int]-fftSizes = [32,128,512,2048]
changelog.md view
@@ -1,11 +1,97 @@+## Changes in 0.16.5.0 [2026.01.09]++ * `ContGen` and `DiscreteGen` instances for `Poisson` distributions are added.+++## Changes in 0.16.4.0 [2025.10.23]++ * Bartlett's test (`Statistics.Test.Bartlett`) and Levene's test+ (`Statistics.Test.Levene`) for homogeneity of variances is added.++ * Improved performance in calculation of moments.++ * Improved precision in calculation of `logDensity` of Student T distribution.+++## 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`.++ * Compatibility with QuickCheck-2.14. Test suite doesn't fail every time.+++## Changes in 0.16.2.0++ * Improved precision for `complCumulative` for hypergeometric and binomial+ distributions. Precision improvements of geometric distribution++ * Negative binomial distribution added.+++## Changes in 0.16.1.2++ * Fixed bug in `fromSample` for exponential distribudion (#190)+++## Changes in 0.16.1.0++ * Dependency on monad-par is dropped. `parMap` from `parallel` is used instead.+++## Changes in 0.16.0.2++ * Bug in constructor of binomial distribution is fixed (#181). It accepted+ out-of range probability before.+++## Changes in 0.16.0.0++ * Random number generation switched to API introduced in random-1.2++ * Support of GHC<7.10 is dropped++ * Fix for chi-squared test (#167) which was completely wrong++ * Computation of CDF and quantiles of Cauchy distribution is now numerically+ stable.++ * Fix loss of precision in computing of CDF of gamma distribution++ * Log-normal and Weibull distributions added.++ * `DiscreteGen` instance added for `DiscreteUniform`++ ## Changes in 0.15.2.0 * Test suite is finally fixed (#42, #123). It took very-very-very long time but finally happened. - * Avoid loss of precision when computing CDF for exponential districution.+ * Avoid loss of precision when computing CDF for exponential distribution. - * Avoid loss of precision when computing CDF for geometric districution. Add+ * Avoid loss of precision when computing CDF for geometric distribution. Add complement of CDF. * Correctly handle case of n=0 in poissonCI@@ -59,7 +145,7 @@ ## Changes in 0.14.0.0 Breaking update. It seriously changes parts of API. It adds new data types for-dealing with with estimates, confidence intervals, confidence levels and+dealing with estimates, confidence intervals, confidence levels and p-value. Also API for statistical tests is changed. * Module `Statistis.Types` now contains new data types for estimates,@@ -231,19 +317,19 @@ * Bugs in DCT and IDCT are fixed. - * Accesors for uniform distribution are added.+ * Accessors for uniform distribution are added. - * ContGen instances for all continuous distribtuions are added.+ * ContGen instances for all continuous distributions are added. * Beta distribution is added. - * Constructor for improper gamma distribtuion is added.+ * Constructor for improper gamma distribution is added. * Binomial distribution allows zero trials. * Poisson distribution now accept zero parameter. - * Integer overflow in caculation of Wilcoxon-T test is fixed.+ * Integer overflow in calculation of Wilcoxon-T test is fixed. * Bug in 'ContGen' instance for normal distribution is fixed. @@ -284,7 +370,7 @@ * Root finding is added, in S.Math.RootFinding. * The complCumulative function is added to the Distribution- class in order to accurately assess probalities P(X>x) which are+ class in order to accurately assess probabilities P(X>x) which are used in one-tailed tests. * A stdDev function is added to the Variance class for@@ -299,7 +385,7 @@ * Bugs in quantile estimations for chi-square and gamma distribution are fixed. - * Integer overlow in mannWhitneyUCriticalValue is fixed. It+ * Integer overflow in mannWhitneyUCriticalValue is fixed. It produced incorrect critical values for moderately large samples. Something around 20 for 32-bit machines and 40 for 64-bit ones.@@ -321,18 +407,18 @@ * Mean and variance for gamma distribution are fixed. - * Much faster cumulative probablity functions for Poisson and+ * Much faster cumulative probability functions for Poisson and hypergeometric distributions. * Better density functions for gamma and Poisson distributions. * Student-T, Fisher-Snedecor F-distributions and Cauchy-Lorentz- distrbution are added.+ distribution are added. * The function S.Function.create is removed. Use generateM from the vector package instead. - * Function to perform approximate comparion of doubles is added to+ * Function to perform approximate comparison of doubles is added to S.Function.Comparison * Regularized incomplete beta function and its inverse are added to
statistics.cabal view
@@ -1,5 +1,8 @@+cabal-version: 3.0+build-type: Simple+ name: statistics-version: 0.15.2.0+version: 0.16.5.0 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful@@ -22,43 +25,52 @@ * Common statistical tests for significant differences between samples. -license: BSD2+license: BSD-2-Clause license-file: LICENSE-homepage: https://github.com/bos/statistics-bug-reports: https://github.com/bos/statistics/issues+homepage: https://github.com/haskell/statistics+bug-reports: https://github.com/haskell/statistics/issues author: Bryan O'Sullivan <bos@serpentine.com>, Alexey Khudaykov <alexey.skladnoy@gmail.com>-maintainer: Bryan O'Sullivan <bos@serpentine.com>, Alexey Khudaykov <alexey.skladnoy@gmail.com>+maintainer: Alexey Khudaykov <alexey.skladnoy@gmail.com> copyright: 2009-2014 Bryan O'Sullivan category: Math, Statistics-build-type: Simple-cabal-version: >= 1.8+ 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 +extra-doc-files:+ changelog.md+ tested-with:- GHC ==7.4.2- || ==7.6.3- || ==7.8.4- || ==7.10.3- || ==8.0.2- || ==8.2.2- || ==8.4.4- || ==8.6.5- || ==8.8.1- , GHCJS ==8.4+ GHC ==8.4.4+ || ==8.6.5+ || ==8.8.4+ || ==8.10.7+ || ==9.0.2+ || ==9.2.8+ || ==9.4.8+ || ==9.6.7+ || ==9.8.4+ || ==9.10.2+ || ==9.12.2 +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: Statistics.Autocorrelation Statistics.ConfidenceInt@@ -76,11 +88,14 @@ Statistics.Distribution.Geometric Statistics.Distribution.Hypergeometric Statistics.Distribution.Laplace+ Statistics.Distribution.Lognormal+ Statistics.Distribution.NegativeBinomial Statistics.Distribution.Normal Statistics.Distribution.Poisson Statistics.Distribution.StudentT Statistics.Distribution.Transform Statistics.Distribution.Uniform+ Statistics.Distribution.Weibull Statistics.Function Statistics.Quantile Statistics.Regression@@ -93,6 +108,8 @@ Statistics.Sample.KernelDensity.Simple Statistics.Sample.Normalize Statistics.Sample.Powers+ Statistics.Test.Bartlett+ Statistics.Test.Levene Statistics.Test.ChiSquared Statistics.Test.KolmogorovSmirnov Statistics.Test.KruskalWallis@@ -108,11 +125,11 @@ Statistics.Internal Statistics.Test.Internal Statistics.Types.Internal- build-depends: base >= 4.5 && < 5- , base-orphans >= 0.6 && <0.9+ build-depends: base >= 4.9 && < 5 --- , math-functions >= 0.3- , mwc-random >= 0.13.0.0+ , math-functions >= 0.3.4.1+ , mwc-random >= 0.15.3.0+ , random >= 1.2 -- , aeson >= 0.6.0.0 , async >= 2.2.2 && <2.3@@ -120,21 +137,21 @@ , binary >= 0.5.1.0 , primitive >= 0.3 , dense-linear-algebra >= 0.1 && <0.2+ , parallel >= 3.2.2.0 && <3.4 , vector >= 0.10 , vector-algorithms >= 0.4 , vector-th-unbox , vector-binary-instances >= 0.2.1 , data-default-class >= 0.1.2- if !impl(ghcjs)- build-depends:- monad-par >= 0.3.4+ -- Older GHC if impl(ghc < 7.6) build-depends: ghc-prim ghc-options: -O2 -Wall -fwarn-tabs -funbox-strict-fields -test-suite tests+test-suite statistics-tests+ default-language: Haskell2010 type: exitcode-stdio-1.0 hs-source-dirs: tests main-is: tests.hs@@ -142,6 +159,7 @@ Tests.ApproxEq Tests.Correlation Tests.Distribution+ Tests.ExactDistribution Tests.Function Tests.Helpers Tests.KDE@@ -156,6 +174,8 @@ Tests.Quantile ghc-options: -Wall -threaded -rtsopts -fsimpl-tick-factor=500+ if impl(ghc >= 9.8)+ ghc-options: -Wno-x-partial build-depends: base , statistics , dense-linear-algebra@@ -165,7 +185,6 @@ , aeson , ieee754 >= 0.7.3 , math-functions- , mwc-random , primitive , tasty , tasty-hunit@@ -174,6 +193,49 @@ , vector , vector-algorithms -source-repository head- type: git- location: https://github.com/bos/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.25++-- 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: Main.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: Main.hs+ Other-modules: Bench+ build-depends: tasty-papi >= 0.1.2
tests/Tests/Correlation.hs view
@@ -5,11 +5,11 @@ import Control.Arrow (Arrow(..)) import qualified Data.Vector as V+import Data.Maybe import Statistics.Correlation import Statistics.Correlation.Kendall-import Test.QuickCheck ((==>),Property,counterexample) import Test.Tasty-import Test.Tasty.QuickCheck+import Test.Tasty.QuickCheck hiding (sample) import Test.Tasty.HUnit import Tests.ApproxEq@@ -34,15 +34,19 @@ testPearson :: [(Double,Double)] -> Property testPearson sample- = (length sample > 1) ==> (exact ~= fast)+ = (length sample > 1 && isJust exact) ==> (case exact of+ Just e -> e ~= fast+ Nothing -> property False+ ) where (~=) = eql 1e-12 exact = exactPearson $ map (realToFrac *** realToFrac) sample fast = pearson $ V.fromList sample -exactPearson :: [(Rational,Rational)] -> Double+exactPearson :: [(Rational,Rational)] -> Maybe Double exactPearson sample- = realToFrac cov / sqrt (realToFrac (varX * varY))+ | varX == 0 || varY == 0 = Nothing+ | otherwise = Just $ realToFrac cov / sqrt (realToFrac (varX * varY)) where (xs,ys) = unzip sample n = fromIntegral $ length sample@@ -98,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@@ -113,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/Distribution.hs view
@@ -2,10 +2,10 @@ ViewPatterns #-} module Tests.Distribution (tests) where -import Control.Applicative ((<$), (<$>), (<*>)) import qualified Control.Exception as E import Data.List (find) import Data.Typeable (Typeable)+import Data.Word import Numeric.MathFunctions.Constants (m_tiny,m_huge,m_epsilon) import Numeric.MathFunctions.Comparison import Statistics.Distribution@@ -19,11 +19,14 @@ import Statistics.Distribution.Geometric import Statistics.Distribution.Hypergeometric import Statistics.Distribution.Laplace (LaplaceDistribution)+import Statistics.Distribution.Lognormal (LognormalDistribution)+import Statistics.Distribution.NegativeBinomial (NegativeBinomialDistribution) import Statistics.Distribution.Normal (NormalDistribution) import Statistics.Distribution.Poisson (PoissonDistribution) import Statistics.Distribution.StudentT import Statistics.Distribution.Transform (LinearTransform) import Statistics.Distribution.Uniform (UniformDistribution)+import Statistics.Distribution.Weibull (WeibullDistribution) import Statistics.Distribution.DiscreteUniform (DiscreteUniform) import Test.Tasty (TestTree, testGroup) import Test.Tasty.QuickCheck (testProperty)@@ -33,6 +36,7 @@ import Text.Printf (printf) import Tests.ApproxEq (ApproxEq(..))+import Tests.ExactDistribution (exactDistributionTests) import Tests.Helpers (T(..), Double01(..), testAssertion, typeName) import Tests.Helpers (monotonicallyIncreasesIEEE,isDenorm) import Tests.Orphanage ()@@ -46,8 +50,10 @@ , contDistrTests (T :: T ExponentialDistribution ) , contDistrTests (T :: T GammaDistribution ) , contDistrTests (T :: T LaplaceDistribution )+ , contDistrTests (T :: T LognormalDistribution ) , contDistrTests (T :: T NormalDistribution ) , contDistrTests (T :: T UniformDistribution )+ , contDistrTests (T :: T WeibullDistribution ) , contDistrTests (T :: T StudentT ) , contDistrTests (T :: T (LinearTransform NormalDistribution)) , contDistrTests (T :: T FDistribution )@@ -56,9 +62,11 @@ , discreteDistrTests (T :: T GeometricDistribution ) , discreteDistrTests (T :: T GeometricDistribution0 ) , discreteDistrTests (T :: T HypergeometricDistribution )+ , discreteDistrTests (T :: T NegativeBinomialDistribution ) , discreteDistrTests (T :: T PoissonDistribution ) , discreteDistrTests (T :: T DiscreteUniform ) + , exactDistributionTests , unitTests ] @@ -71,8 +79,7 @@ contDistrTests t = testGroup ("Tests for: " ++ typeName t) $ cdfTests t ++ [ testProperty "PDF sanity" $ pdfSanityCheck t- ] ++- [ (if quantileIsInvCDF_enabled t then id else ignoreTest)+ , (if quantileIsInvCDF_enabled t then id else ignoreTest) $ testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t , testProperty "log density check" $ logDensityCheck t@@ -86,7 +93,7 @@ [ testProperty "Prob. sanity" $ probSanityCheck t , testProperty "CDF is sum of prob." $ discreteCDFcorrect t , testProperty "Discrete CDF is OK" $ cdfDiscreteIsCorrect t- , testProperty "log probabilty check" $ logProbabilityCheck t+ , testProperty "log probability check" $ logProbabilityCheck t ] -- Tests for distributions which have CDF@@ -94,7 +101,8 @@ cdfTests t = [ testProperty "C.D.F. sanity" $ cdfSanityCheck t , testProperty "CDF limit at +inf" $ cdfLimitAtPosInfinity t- , testProperty "CDF limit at -inf" $ cdfLimitAtNegInfinity t+ , (if cdfLimitAtNegInfinity_enabled t then id else ignoreTest)+ $ testProperty "CDF limit at -inf" $ cdfLimitAtNegInfinity t , testProperty "CDF at +inf = 1" $ cdfAtPosInfinity t , testProperty "CDF at -inf = 1" $ cdfAtNegInfinity t , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing t@@ -143,11 +151,14 @@ -- CDF's complement is implemented correctly-cdfComplementIsCorrect :: (Distribution d, Param d) => T d -> d -> Double -> Bool+cdfComplementIsCorrect :: (Distribution d, Param d) => T d -> d -> Double -> Property cdfComplementIsCorrect _ d x- = 1 - (cumulative d x + complCumulative d x) <= tol+ = counterexample ("err. tolerance = " ++ show tol)+ $ counterexample ("difference = " ++ show delta)+ $ delta <= tol where- tol = prec_complementCDF d+ tol = prec_complementCDF d+ delta = 1 - (cumulative d x + complCumulative d x) -- CDF for discrete distribution uses <= for comparison cdfDiscreteIsCorrect :: (Param d, DiscreteDistr d) => T d -> d -> Property@@ -159,38 +170,43 @@ -- -- > CDF(i) - CDF(i-e) = P(i) --- -- Apporixmate equality is tricky here. Scale is set by maximum- -- value of CDF and probability. Case when all proabilities are- -- zero should be trated specially.+ -- Approximate equality is tricky here. Scale is set by maximum+ -- value of CDF and probability. Case when all probabilities are+ -- zero should be treated specially. badN = [ printf "N=%3i p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp) | i <- [0 .. 100] , let p = cumulative d $ fromIntegral i - 1e-6 p1 = cumulative d $ fromIntegral i dp = probability d i relerr = ((p1 - p) - dp) / max p1 dp- , not (p == 0 && p1 == 0 && dp == 0)- && relerr > tol+ , p > m_tiny || p == 0+ , p1 > m_tiny+ , dp > m_tiny+ , relerr > tol ] tol = prec_discreteCDF d -logDensityCheck :: (ContDistr d) => T d -> d -> Double -> Property+logDensityCheck :: (Param d, ContDistr d) => T d -> d -> Double -> Property logDensityCheck _ d x = not (isDenorm x) ==> ( counterexample (printf "density = %g" p) $ counterexample (printf "logDensity = %g" logP) $ counterexample (printf "log p = %g" (log p))- $ counterexample (printf "eps = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))+ $ counterexample (printf "ulps[log] = %i" ulpsLog)+ $ counterexample (printf "ulps[lin] = %i" ulpsLin) $ or [ p == 0 && logP == (-1/0) , p <= m_tiny && logP < log m_tiny -- To avoid problems with roundtripping error in case -- when density is computed as exponent of logDensity we -- accept either inequality- , (ulpDistance (log p) logP <= 32)- || (ulpDistance p (exp logP) <= 32)+ , (ulpsLog <= n) || (ulpsLin <= n) ]) where- p = density d x- logP = logDensity d x+ p = density d x+ logP = logDensity d x+ n = prec_logDensity d+ ulpsLog = ulpDistance (log p) logP+ ulpsLin = ulpDistance p (exp logP) -- PDF is positive pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool@@ -201,6 +217,8 @@ complQuantileCheck _ d (Double01 p) = counterexample (printf "x0 = %g" x0) $ counterexample (printf "x1 = %g" x1)+ $ counterexample (printf "abs err = %g" $ abs (x1 - x0))+ $ counterexample (printf "rel err = %g" $ relativeError x1 x0) -- We avoid extreme tails of distributions -- -- FIXME: all parameters are arbitrary at the moment@@ -208,7 +226,7 @@ , p < 0.99 , not $ isInfinite x0 , not $ isInfinite x1- ] ==> (abs (x1 - x0) < 1e-6)+ ] ==> (if x0 < 1e6 then abs (x1 - x0) < 1e-6 else relativeError x1 x0 < 1e-12) where x0 = quantile d (1 - p) x1 = complQuantile d p@@ -263,9 +281,9 @@ $ counterexample (printf "Sum = %g" p2) $ counterexample (printf "Delta = %g" (abs (p1 - p2))) $ abs (p1 - p2) < 3e-10- -- Avoid too large differeneces. Otherwise there is to much to sum+ -- Avoid too large differences. Otherwise there is to much to sum --- -- Absolute difference is used guard againist precision loss when+ -- Absolute difference is used guard against precision loss when -- close values of CDF are subtracted where n = min a b@@ -273,23 +291,26 @@ p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5) p2 = sum $ map (probability d) [n .. m] -logProbabilityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Property+logProbabilityCheck :: (Param d, DiscreteDistr d) => T d -> d -> Int -> Property logProbabilityCheck _ d x = counterexample (printf "probability = %g" p) $ counterexample (printf "logProbability = %g" logP) $ counterexample (printf "log p = %g" (log p))- $ counterexample (printf "eps = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))+ $ counterexample (printf "ulps[log] = %i" ulpsLog)+ $ counterexample (printf "ulps[lin] = %i" ulpsLin) $ or [ p == 0 && logP == (-1/0) , p < 1e-308 && logP < 609 -- To avoid problems with roundtripping error in case -- when density is computed as exponent of logDensity we -- accept either inequality- , (ulpDistance (log p) logP <= 32)- || (ulpDistance p (exp logP) <= 32)+ , (ulpsLog <= n) || (ulpsLin <= n) ] where p = probability d x logP = logProbability d x+ n = prec_logDensity d+ ulpsLog = ulpDistance (log p) logP+ ulpsLin = ulpDistance p (exp logP) -- | Parameters for distribution testing. Some distribution require@@ -298,6 +319,9 @@ -- | Whether quantileIsInvCDF is enabled quantileIsInvCDF_enabled :: T a -> Bool quantileIsInvCDF_enabled _ = True+ -- | Whether cdfLimitAtNegInfinity is enabled+ cdfLimitAtNegInfinity_enabled :: T a -> Bool+ cdfLimitAtNegInfinity_enabled _ = True -- | Precision for 'quantileIsInvCDF' test prec_quantile_CDF :: a -> (Double,Double) prec_quantile_CDF _ = (16,16)@@ -307,36 +331,57 @@ -- | Precision of CDF's complement prec_complementCDF :: a -> Double prec_complementCDF _ = 1e-14+ -- | Precision for logDensity check+ prec_logDensity :: a -> Word64+ prec_logDensity _ = 32 instance Param StudentT where -- FIXME: disabled unless incompleteBeta troubles are sorted out quantileIsInvCDF_enabled _ = False+ instance Param BetaDistribution where- -- FIXME: See https://github.com/bos/statistics/issues/161 for details+ -- FIXME: See https://github.com/haskell/statistics/issues/161 for details quantileIsInvCDF_enabled _ = False+ instance Param FDistribution where -- FIXME: disabled unless incompleteBeta troubles are sorted out quantileIsInvCDF_enabled _ = False+ -- We compute CDF and complement using same method so precision+ -- should be very good here.+ prec_complementCDF _ = 64 * m_epsilon instance Param ChiSquared where prec_quantile_CDF _ = (32,32) instance Param BinomialDistribution where- prec_discreteCDF _ = 1e-13-instance Param CauchyDistribution+ prec_discreteCDF _ = 1e-12+ prec_logDensity _ = 48+instance Param CauchyDistribution where+ -- Distribution is long-tailed enough that we may never get to zero+ cdfLimitAtNegInfinity_enabled _ = False+ instance Param DiscreteUniform instance Param ExponentialDistribution-instance Param GammaDistribution+instance Param GammaDistribution where+ -- We lose precision near `incompleteGamma 10` because of error+ -- introduced by exp . logGamma. This could only be fixed in+ -- math-function by implementing gamma+ prec_quantile_CDF _ = (24,24)+ prec_logDensity _ = 512 instance Param GeometricDistribution instance Param GeometricDistribution0 instance Param HypergeometricDistribution instance Param LaplaceDistribution+instance Param LognormalDistribution where+ prec_quantile_CDF _ = (64,64)+instance Param NegativeBinomialDistribution where+ prec_discreteCDF _ = 1e-12+ prec_logDensity _ = 48 instance Param NormalDistribution instance Param PoissonDistribution instance Param UniformDistribution+instance Param WeibullDistribution instance Param a => Param (LinearTransform a)-- ---------------------------------------------------------------- -- Unit tests
+ tests/Tests/ExactDistribution.hs view
@@ -0,0 +1,387 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE TypeFamilies #-}+-- |+-- Module : Tests.ExactDistribution+-- Copyright : (c) 2022 Lorenz Minder+-- License : BSD3+--+-- Maintainer : lminder@gmx.net+-- Stability : experimental+-- Portability : portable+--+-- Tests comparing distributions to exact versions.+--+-- This module provides exact versions of some distributions, and tests+-- to compare them to the production implementations in+-- Statistics.Distribution.*. It also contains the functionality to+-- test the production distributions against the exact versions. Errors+-- are flagged if data points are discovered where the probability mass+-- function, the cumulative probability function, or its complement+-- deviates too far (more than a prescribed tolerance) from the exact+-- calculation.+--+-- The distributions here are implemented with rational integer+-- arithmetic, using pretty much the textbook definitions formulas.+-- Numerical problems like overflow or rounding errors cannot occur with+-- this approach, making them are easy to write, read and verify. They+-- are, of course, substantially slower than the production+-- distributions in Statistics.Distribution.*. This makes them+-- unsuitable for most uses other than testing and debugging. (Also,+-- only a handful of distributions can be implemented exactly with+-- rational arithmetic.)+--+-- This module has the following sub-components:+-- +-- * Exact (rational) definitions of some distribution functions,+-- including both the probability mass as well as the CDF.+--+-- * QC.Arbitrary implementations to sample test cases (i.e.,+-- distribution parameters and evaluation points).+--+-- * "Linkage": a mechanism to construct a production distribution+-- corresponding to a test case for an exact distribution.+--+-- * A set of tests for the distributions derived using all of the above+-- components.+--+-- This module exports a number symbols which can be useful for+-- debugging and experimentation. For use in a test suite, only the+-- `exactDistributionTests` function is needed.++module Tests.ExactDistribution (+ -- * Exact math functions+ exactChoose++ -- * Exact distributions+ , ExactDiscreteDistr(..)++ , ExactBinomialDistr(..)+ , ExactDiscreteUniformDistr(..)+ , ExactGeometricDistr(..)+ , ExactHypergeomDistr(..)++ -- * Linking to production distributions+ , ProductionLinkage++ -- * Individual test routines+ , pmfMatch+ , cdfMatch+ , complCdfMatch++ -- * Test groups+ , Tag(..)+ , distTests+ , exactDistributionTests+) where++----------------------------------------------------------------++import Data.Foldable+import Data.Ratio++import Test.Tasty (TestTree, testGroup)+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+import Statistics.Distribution.DiscreteUniform+import Statistics.Distribution.Geometric+import Statistics.Distribution.Hypergeometric++----------------------------------------------------------------+--+-- Math functions.+--+-- Used for implementing the distributions below.+--+----------------------------------------------------------------++-- | Exactly compute binomial coefficient.+--+-- /n/ need not be an integer, can be fractional.+exactChoose :: Ratio Integer -> Integer -> Ratio Integer+exactChoose n k+ | k < 0 = 0+ | otherwise = foldl' (*) 1 factors+ where factors = [ (n - k' + j) / j | j <- [1..k'] ]+ k' = fromInteger k :: Ratio Integer++----------------------------------------------------------------+--+-- Exact distributions.+--+----------------------------------------------------------------++-- | Exact discrete distribution.+class ExactDiscreteDistr a where+ -- | Probability mass function.+ exactProb :: a -> Integer -> Ratio Integer+ exactProb d x = exactCumulative d x - exactCumulative d (x - 1)++ -- | Cumulative distribution function.+ exactCumulative :: a -> Integer -> Ratio Integer++-- | Exact Binomial distribution.+data ExactBinomialDistr = ExactBD Integer (Ratio Integer)+ deriving(Show)++instance ExactDiscreteDistr ExactBinomialDistr where+ -- Probability mass, computed with textbook formula.+ exactProb (ExactBD n p) k+ | k < 0 || k > n = 0+ | otherwise = exactChoose n' k * p^k * (1-p)^(n-k)+ where n' = fromIntegral n+ -- CDF + --+ -- Computed iteratively by summing up all the probabilities+ -- <= /k/. Rather than computing everything from scratch for each+ -- probability, we reuse previous results. The meanings of the+ -- variables in the "update" function are:+ -- + -- bc is the binomial coefficient (n choose j),+ -- pj is the term p^j,+ -- pnj is the term (1 - p)^(n - j)+ -- r is the (partial) sum of the probabilities + --+ exactCumulative (ExactBD n p) k+ | k < 0 = 0+ | k >= n = 1+ -- Special case for p = 1, since in the below fold we+ -- divide by (1 - p).+ | p == 1 = if k == n then 1 else 0+ | otherwise+ = result $ foldl' update (1, 1, (1 - p)^n, (1 - p)^n) [1..k]+ where update (!bc, !pj, !pnj, !r) !j =+ let bc' = bc * (n - j + 1) `div` j + pj' = pj * p+ pnj' = pnj / (1 - p)+ r' = r + (fromIntegral bc') * pj' * pnj'+ in (bc', pj', pnj', r')+ result (_, _, _, r) = r++-- | Exact Discrete Uniform distribution.+data ExactDiscreteUniformDistr = ExactDU Integer Integer+ deriving(Show)++instance ExactDiscreteDistr ExactDiscreteUniformDistr where+ exactProb (ExactDU lower upper) k+ | k < lower || k > upper = 0+ | otherwise = 1 % (upper - lower + 1)+ exactCumulative (ExactDU lower upper) k+ | k < lower = 0+ | k > upper = 1+ | otherwise =+ let d = (k - lower + 1)+ in d % (upper - lower + 1)++-- | Geometric distribution.+data ExactGeometricDistr = ExactGeom (Ratio Integer)+ deriving(Show)++instance ExactDiscreteDistr ExactGeometricDistr where+ exactProb (ExactGeom p) k+ | k < 1 = 0+ | otherwise = (1 - p)^(k - 1) * p++ exactCumulative (ExactGeom p) k = 1 - (1 - p)^k++-- | Hypergeometric distribution.+--+-- Parameters are /K/, /N/ and /n/, where:+-- - /N/ is the total sample space size.+-- - /K/ is number of "good" objects among /N/.+-- - /n/ is the number of draws without replacement.+data ExactHypergeomDistr = ExactHG Integer Integer Integer+ deriving(Show)++instance ExactDiscreteDistr ExactHypergeomDistr where+ exactProb (ExactHG nK nN n) k+ | k < 0 = 0+ | k > n || k > nN = 0+ | otherwise =+ exactChoose nK' k * exactChoose (nN' - nK') (n - k)+ / exactChoose nN' n+ where nN' = fromIntegral nN+ nK' = fromIntegral nK++ exactCumulative d k = sum [ exactProb d i | i <- [0..k] ]++----------------------------------------------------------------+--+-- TestCase construction.+--+-- Contains the TestCase data type which encapsulates an instance of an+-- exact distribution together with an evaluation point.+--+-- Then in contains the QC.Arbitrary implementations for TestCases of+-- the different exact distributions. As a general rule, we try the+-- sampling to be relatively efficient, i.e., we only want to sample+-- valid distribution parameters. The evaluation points are sampled+-- such that most points are within the support of the distribution.+--+----------------------------------------------------------------++-- Divisor to compute a rational number from an integer.+--+-- We want input parameters to be exactly representable as+-- Double values. This is so that the production distribution does not+-- mismatch the exact one simply because the input values don't exactly+-- match. (This can happen if the derivative of the distribution+-- function is large.) For this reason, the gd value needs to be a+-- power of 2, and <= 2^53, since the mantissa of a Double is 53 bits.+--+-- A value of 2^53 gives the most accurate and diverse tests, but the+-- cost is increased running times, as the computed numerators and+-- denominators will become quite large.+gd :: Integer+gd = 2^(16 :: Int)++-- TestCase+--+-- Combination of an exact distribution together with an evaluation point.+data TestCase a = TestCase a Integer deriving (Show)++instance QC.Arbitrary (TestCase ExactBinomialDistr) where+ arbitrary = do+ -- This somewhat odd sampling of /n/ is done so that lower+ -- values (<1000) are more often represented as the larger ones.+ n <- (*) <$> chooseInteger (1,1000) <*> chooseInteger(1,2)+ p <- (% gd) <$> chooseInteger (0, gd)+ k <- chooseInteger (-1, n + 1)+ return $ TestCase (ExactBD n p) k+ shrink _ = []++instance QC.Arbitrary (TestCase ExactDiscreteUniformDistr) where+ arbitrary = do+ a <- chooseInteger (-1000, 1000)+ sz <- chooseInteger (1, 1000)+ let b = a + sz+ k <- chooseInteger (a - 10, b + 10)+ return $ TestCase (ExactDU a b) k+ shrink _ = []++instance QC.Arbitrary (TestCase ExactGeometricDistr) where+ arbitrary = do+ p <- (% gd) <$> chooseInteger (1, gd)+ let lim = (floor $ 100 / p) :: Integer+ k <- chooseInteger (0, lim)+ return $ TestCase (ExactGeom p) k+ shrink _ = []++instance QC.Arbitrary (TestCase ExactHypergeomDistr) where+ arbitrary = do+ nN <- chooseInteger (1, 100) -- XXX lower bound should be 0+ nK <- chooseInteger (0, nN)+ n <- chooseInteger (1, nN) -- XXX lower bound should be 0+ k <- chooseInteger (0, min n nK)+ return $ TestCase (ExactHG nK nN n) k+ shrink _ = []++----------------------------------------------------------------+--+-- Linking to the production distributions+--+-- This section contains the ProductionLinkage typeclass and+-- implementation, that allows to obtain a functions for evaluating+-- the production distribution functions for a corresponding exact+-- distribution.+--+----------------------------------------------------------------++class (ExactDiscreteDistr a, DiscreteDistr (ProdDistrib a)+ ) => ProductionLinkage a where+ type ProdDistrib a+ toProd :: a -> ProdDistrib a++instance ProductionLinkage ExactBinomialDistr where+ type ProdDistrib ExactBinomialDistr = BinomialDistribution+ toProd (ExactBD n p) = binomial (fromIntegral n) (fromRational p)++instance ProductionLinkage ExactDiscreteUniformDistr where+ type ProdDistrib ExactDiscreteUniformDistr = DiscreteUniform+ toProd (ExactDU lower upper) = discreteUniformAB (fromIntegral lower) (fromIntegral upper)++instance ProductionLinkage ExactGeometricDistr where+ type ProdDistrib ExactGeometricDistr = GeometricDistribution+ toProd (ExactGeom p) = geometric $ fromRational p++instance ProductionLinkage ExactHypergeomDistr where+ 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 -> 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)+ = 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)+ = 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 :: 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 @a tol+ , testProperty "CDF match" $ cdfMatch @a tol+ , testProperty "1 - CDF match" $ complCdfMatch @a tol+ ]+++-- Test driver -------------------------------------------------++exactDistributionTests :: TestTree+exactDistributionTests = testGroup "Test distributions against exact"+ [ 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/Helpers.hs view
@@ -60,8 +60,8 @@ -- 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+-- In fact function is allowed to decrease less than one ulp in order+-- to guard against 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
− 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/Tests/Matrix.hs view
@@ -5,7 +5,6 @@ import Test.Tasty (TestTree, testGroup) import Test.Tasty.QuickCheck (testProperty) import Test.QuickCheck-import Tests.ApproxEq (ApproxEq(..)) import Tests.Matrix.Types import qualified Data.Vector.Unboxed as U @@ -27,9 +26,20 @@ t_transpose m = U.concat (map (column n) [0..rows m-1]) === toVector m where n = transpose m -t_qr :: Matrix -> Property-t_qr a = hasNaN p .||. eql 1e-10 a p- where p = uncurry multiply (qr a)+t_qr :: Property+t_qr = property $ do+ a <- do (r,c) <- arbitrary+ fromMat <$> arbMatWith r c (fromIntegral <$> choose (-10, 10::Int))+ let (q,r) = qr a+ a' = multiply q r+ pure $ counterexample ("A = \n"++show a)+ $ counterexample ("A' = \n"++show a')+ $ counterexample ("Q = \n"++show q)+ $ counterexample ("R = \n"++show r)+ $ dimension a == dimension a'+ && ( hasNaN a'+ || and (zipWith (\x y -> abs (x - y) < 1e-12) (toList a) (toList a'))+ ) tests :: TestTree tests = testGroup "Matrix"
tests/Tests/Matrix/Types.hs view
@@ -6,6 +6,8 @@ Mat(..) , fromMat , toMat+ , arbMat+ , arbMatWith ) where import Control.Monad (join)@@ -32,10 +34,21 @@ arbitrary = small $ join (arbMat <$> arbitrary <*> arbitrary) shrink (Mat r c xs) = Mat r c <$> shrinkFixedList (shrinkFixedList shrink) xs -arbMat :: (Arbitrary a) => Positive (Small Int) -> Positive (Small Int)- -> Gen (Mat a)-arbMat (Positive (Small r)) (Positive (Small c)) =- Mat r c <$> vectorOf r (vector c)+arbMat+ :: (Arbitrary a)+ => Positive (Small Int)+ -> Positive (Small Int)+ -> Gen (Mat a)+arbMat r c = arbMatWith r c arbitrary++arbMatWith+ :: (Arbitrary a)+ => Positive (Small Int)+ -> Positive (Small Int)+ -> Gen a+ -> Gen (Mat a)+arbMatWith (Positive (Small r)) (Positive (Small c)) genA =+ Mat r c <$> vectorOf r (vectorOf c genA) instance Arbitrary Matrix where arbitrary = fromMat <$> arbitrary
tests/Tests/Orphanage.hs view
@@ -16,11 +16,14 @@ import Statistics.Distribution.Geometric import Statistics.Distribution.Hypergeometric import Statistics.Distribution.Laplace (LaplaceDistribution, laplace)+import Statistics.Distribution.Lognormal (LognormalDistribution, lognormalDistr)+import Statistics.Distribution.NegativeBinomial (NegativeBinomialDistribution, negativeBinomial) import Statistics.Distribution.Normal (NormalDistribution, normalDistr) import Statistics.Distribution.Poisson (PoissonDistribution, poisson) import Statistics.Distribution.StudentT import Statistics.Distribution.Transform (LinearTransform, scaleAround) import Statistics.Distribution.Uniform (UniformDistribution, uniformDistr)+import Statistics.Distribution.Weibull (WeibullDistribution, weibullDistr) import Statistics.Distribution.DiscreteUniform (DiscreteUniform, discreteUniformAB) import Statistics.Types @@ -28,7 +31,7 @@ ------------------------------------------------------------------- Arbitrary instances for ditributions+-- Arbitrary instances for distributions ---------------------------------------------------------------- instance QC.Arbitrary BinomialDistribution where@@ -38,18 +41,23 @@ instance QC.Arbitrary LaplaceDistribution where arbitrary = laplace <$> QC.choose (-10,10) <*> QC.choose (0, 2) instance QC.Arbitrary GammaDistribution where- arbitrary = gammaDistr <$> QC.choose (0.1,10) <*> QC.choose (0.1,10)+ arbitrary = gammaDistr <$> QC.choose (0.1,100) <*> QC.choose (0.1,100) instance QC.Arbitrary BetaDistribution where arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10) instance QC.Arbitrary GeometricDistribution where- arbitrary = geometric <$> QC.choose (0,1)+ arbitrary = geometric <$> QC.choose (1e-10,1) instance QC.Arbitrary GeometricDistribution0 where- arbitrary = geometric0 <$> QC.choose (0,1)+ arbitrary = geometric0 <$> QC.choose (1e-10,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 LognormalDistribution where+ -- can't choose sigma too big, otherwise goes outside of double-float limit+ arbitrary = lognormalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-10, 20)+instance QC.Arbitrary NegativeBinomialDistribution where+ arbitrary = negativeBinomial <$> QC.choose (1,100) <*> QC.choose (1e-10,1) instance QC.Arbitrary NormalDistribution where arbitrary = normalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-3, 1e3) instance QC.Arbitrary PoissonDistribution where@@ -60,6 +68,8 @@ arbitrary = do a <- QC.arbitrary b <- QC.arbitrary `suchThat` (/= a) return $ uniformDistr a b+instance QC.Arbitrary WeibullDistribution where+ arbitrary = weibullDistr <$> QC.choose (1e-3,1e3) <*> QC.choose (1e-3, 1e3) instance QC.Arbitrary CauchyDistribution where arbitrary = cauchyDistribution <$> arbitrary
tests/Tests/Parametric.hs view
@@ -4,12 +4,18 @@ import Statistics.Test.StudentT import Statistics.Types import qualified Data.Vector.Unboxed as U-import Test.Tasty (testGroup)-import Tests.Helpers (testEquality)+import qualified Data.Vector as V+import Test.Tasty (testGroup, TestTree)+import Test.Tasty.HUnit (testCase, assertBool)+import Tests.Helpers (testEquality) import qualified Test.Tasty as Tst +import Statistics.Test.Levene+import Statistics.Test.Bartlett++ tests :: Tst.TestTree-tests = testGroup "Parametric tests" studentTTests+tests = testGroup "Parametric tests" [studentTTests, bartlettTests, leveneTests] -- 2 samples x 20 obs data --@@ -77,9 +83,9 @@ , testEquality name (isSignificant (mkPValue $ pValue pVal + 1e-5) test) Significant ]- -studentTTests :: [Tst.TestTree]-studentTTests = concat++studentTTests :: Tst.TestTree+studentTTests = testGroup "StudentT test" $ concat [ -- R: t.test(sample1, sample2, alt="two.sided", var.equal=T) testTTest "two-sample t-test SamplesDiffer Student" (mkPValue 0.03410) (fromJust $ studentTTest SamplesDiffer sample1 sample2)@@ -100,3 +106,119 @@ (mkPValue 0.01705) (fromJust $ pairedTTest BGreater sample12) ] where sample12 = U.zip sample1 sample2+++------------------------------------------------------------+-- Bartlett's Test+------------------------------------------------------------++bartlettTests :: TestTree+bartlettTests = testGroup "Bartlett's test"+ [ testCase "a,b,c" $ testBartlettTest [a,b,c] 1.8027132567760222 0.40601846976301237+ , testCase "a,b" $ testBartlettTest [a,b] 0.005221063776321886 0.9423974408021293+ , testCase "a,c" $ testBartlettTest [a,c] 1.1531619271845452 0.2828882244527482+ , testCase "a,a" $ testBartlettTest [a,a] 0.0 1.0+ ]+ where+ a = U.fromList [9.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]+ b = U.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 9.36, 9.18, 8.67, 9.05]+ c = U.fromList [8.95, 8.12, 8.95, 8.85, 8.03, 8.84, 8.07, 8.98, 8.86, 8.98]++testBartlettTest+ :: [U.Vector Double]+ -> Double+ -> Double+ -> IO ()+testBartlettTest samples w p = do+ r <- case bartlettTest samples of+ Left _ -> error "Bartlett's test failed"+ Right r -> pure r+ approxEqual "W" 1e-9 (testStatistics r) w+ approxEqual "p" 1e-9 (pValue $ testSignificance r) p++------------------------------------------------------------+-- Levene's Test (Trimmed Mean)+------------------------------------------------------------++leveneTests :: TestTree+leveneTests = testGroup "Levene test"+ -- Statistics' value and p-values are computed using + [ testCase "a,b,c Mean" $ testLeveneTest [a,b,c] Mean 7.905194483442054 0.001983795817472731+ , testCase "a,b Mean" $ testLeveneTest [a,b] Mean 8.83873787256358 0.008149720958328811+ , testCase "a,a Mean" $ testLeveneTest [a,a] Mean 0.0 1.0+ , testCase "a,b,c Median" $ testLeveneTest [a,b,c] Median 7.584952754501659 0.002431505967249681+ , testCase "a,b Median" $ testLeveneTest [a,b] Median 8.461374333228711 0.009364737715584399+ , testCase "aL,bL Mean" $ testLeveneTest [aL,bL] Mean 5.84424549939465 0.01653410652558999+ , testCase "aL,bL Trimmed" $ testLeveneTest [aL,bL] (Trimmed 0.05) 8.368311226366314 0.004294953946529551+ ]+ where+ a = V.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]+ b = V.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]+ c = V.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]+ -- Large samples for testing trimmed+ aL = V.fromList [+ -0.18919252, -1.62837673, 5.21332355, -0.00962043, -0.28417847,+ -0.88128233, 1.49698436, 6.1780359 , -1.22301348, 3.34598245,+ 5.33227264, -0.88732069, 0.14487346, 2.61060215, 4.22033907,+ 2.53139215, -0.72131061, 0.53063607, -0.60510374, -0.73230842,+ 1.54037043, -2.81103963, 3.40763063, 0.49005324, 2.13085513,+ 5.68650547, 4.16397279, -0.17325097, 1.12664972, 4.23297516,+ 4.15943436, -1.01452078, 2.40391646, 0.83019962, 0.29665879,+ -3.83031046, -1.98576933, 1.5356527 , 1.30773365, 0.292818 ,+ 2.45877828, 1.06482289, -0.63241873, 1.58465379, 1.96577614,+ 2.25791943, 4.13769848, -2.38595767, -0.65801423, -2.54007791,+ 3.17428087, 4.32096964, 0.92240335, -2.38101319, 1.35692587,+ 1.48279101, -0.04438309, 0.50296642, 2.08261495, 1.33181215,+ -1.95427198, 4.95406809, 1.51294898, -2.68536129, -0.2441218 ,+ 2.41142613, 4.71051493, 2.66618697, 1.12668301, -0.25732583,+ 1.25021838, -1.27523641, 5.01638744, 3.38864442, 0.17979744,+ -0.88481645, 3.89346357, -0.51512217, -1.60542888, 0.88378679,+ -2.12962732, -1.35989539, 5.09215112, -1.37442481, 0.83578405,+ 0.13829571, 1.25171481, 3.60552158, -3.24051591, -0.44301834,+ 0.78253445, 1.76098254, 1.79677434, -0.19010505, 3.07640466,+ 3.02853882, 1.24849063, 4.84505382, 6.82274999, 2.24063474]+ bL = V.fromList [+ 2.15584101, -2.74876744, -0.82231894, 1.97518087, 2.59280595,+ 1.28703417, 2.40450278, 1.9761031 , 2.35186598, 1.15611047,+ 2.26709318, 1.2832138 , -2.1486074 , 0.27563011, -0.51816861,+ 0.89658424, 3.27069545, 1.72846646, 3.84454277, 5.58301459,+ -0.40878188, 3.41602853, 1.1281526 , 0.9665913 , 0.76567084,+ 1.69522855, 1.69133014, 0.70529264, 2.65243202, -1.0088019 ,+ -0.62431026, 3.76667396, 3.66225181, 0.73217579, 0.04478736,+ 0.4169833 , 0.77065631, -1.31484093, 1.23858618, -0.08339456,+ 3.14154286, 1.84358218, -0.53511423, -3.4919477 , 0.24076997,+ 3.59381684, 1.99497806, 2.95499775, 1.67157731, 0.0214764 ,+ 3.32161612, -2.64762427, 0.06486472, 0.19653897, 1.34954235,+ 1.18568747, -0.54434597, -3.35544223, 1.41933109, 0.95100195,+ 2.7182116 , 1.1334068 , -0.95297806, -0.05421818, 1.42248799,+ -3.96201277, -3.21309254, -0.21209211, 0.9689551 , 0.13526401,+ -0.88656198, 0.41331783, -3.18766064, 4.34948246, 1.35656384,+ 0.41920101, -0.46578994, 1.55181583, 2.43937014, 2.49040644,+ 4.10505494, 1.68856296, 1.31503895, 0.41123368, 0.73242999,+ 0.2804349 , -1.83494592, -0.31073195, 2.61185513, 2.91645094,+ 1.26097638, 2.64197134, 3.88931972, 0.03783002, 2.55209729,+ 3.46869549, 0.96348003, 2.27658242, 2.7613171 , -0.1372434 ]++ +testLeveneTest+ :: [V.Vector Double]+ -> Center+ -> Double+ -> Double+ -> IO ()+testLeveneTest samples center w p = do+ r <- case levenesTest center samples of+ Left _ -> error "Levene's test failed"+ Right r -> pure r+ approxEqual "W" 1e-9 (testStatistics r) w+ approxEqual "p" 1e-9 (pValue $ testSignificance r) p+++----------------------------------------------------------------++approxEqual :: String -> Double -> Double -> Double -> IO ()+approxEqual name epsilon actual expected =+ assertBool (name ++ ": expected ≈ " ++ show expected ++ ", got " ++ show actual)+ (diff < epsilon)+ where+ diff = abs (actual - expected)
tests/Tests/Serialization.hs view
@@ -16,11 +16,14 @@ import Statistics.Distribution.Geometric import Statistics.Distribution.Hypergeometric import Statistics.Distribution.Laplace (LaplaceDistribution)+import Statistics.Distribution.Lognormal (LognormalDistribution)+import Statistics.Distribution.NegativeBinomial (NegativeBinomialDistribution) import Statistics.Distribution.Normal (NormalDistribution) import Statistics.Distribution.Poisson (PoissonDistribution) import Statistics.Distribution.StudentT import Statistics.Distribution.Transform (LinearTransform) import Statistics.Distribution.Uniform (UniformDistribution)+import Statistics.Distribution.Weibull (WeibullDistribution) import Statistics.Types import Test.Tasty (TestTree, testGroup)@@ -50,8 +53,11 @@ , serializationTests (T :: T ExponentialDistribution ) , serializationTests (T :: T GammaDistribution ) , serializationTests (T :: T LaplaceDistribution )+ , serializationTests (T :: T LognormalDistribution )+ , serializationTests (T :: T NegativeBinomialDistribution ) , serializationTests (T :: T NormalDistribution ) , serializationTests (T :: T UniformDistribution )+ , serializationTests (T :: T WeibullDistribution ) , serializationTests (T :: T StudentT ) , serializationTests (T :: T (LinearTransform NormalDistribution)) , serializationTests (T :: T FDistribution )
tests/Tests/Transform.hs view
@@ -8,7 +8,6 @@ import Data.Bits ((.&.), shiftL) import Data.Complex (Complex((:+)))-import Data.Functor ((<$>)) import Numeric.Sum (kbn, sumVector) import Statistics.Function (within) import Statistics.Transform (CD, dct, fft, idct, ifft)
+ tests/doctest.hs view
@@ -0,0 +1,5 @@+import Test.DocTest (doctest)++main :: IO ()+main = doctest ["-XHaskell2010", "Statistics"]+