statistics 0.16.3.0 → 0.16.4.0
raw patch · 14 files changed
+520/−120 lines, 14 files
Files
- README.markdown +3/−10
- Statistics/Distribution/NegativeBinomial.hs +4/−4
- Statistics/Distribution/StudentT.hs +4/−3
- Statistics/Sample.hs +4/−2
- Statistics/Sample/Internal.hs +7/−2
- Statistics/Test/Bartlett.hs +99/−0
- Statistics/Test/ChiSquared.hs +12/−6
- Statistics/Test/Levene.hs +153/−0
- benchmark/Main.hs +77/−0
- benchmark/bench.hs +0/−75
- changelog.md +13/−3
- statistics.cabal +15/−8
- tests/Tests/Distribution.hs +1/−1
- tests/Tests/Parametric.hs +128/−6
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
Statistics/Distribution/NegativeBinomial.hs view
@@ -47,7 +47,7 @@ gChoose :: Double -> Int -> Double gChoose n k | k < 0 = 0- | k' >= 50 = exp $ logChooseFast n k' + | k' >= 50 = exp $ logChooseFast n k' | otherwise = foldl' (*) 1 factors where factors = [ (n - k' + j) / j | j <- [1..k'] ] k' = fromIntegral k@@ -151,7 +151,7 @@ | k < 0 = 1 | otherwise = incompleteBeta (fromIntegral (k+1)) r (1 - p) where- k = (floor x)::Integer+ k = floor x :: Integer mean :: NegativeBinomialDistribution -> Double mean (NBD r p) = r * (1 - p)/p@@ -166,14 +166,14 @@ dropWhile (>= -m_epsilon) $ [ let x = probability d k in x * log x | k <- [0..]] --- | Construct negative binomial distribution. Number of failures /r/+-- | 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 failures /r/+-- | 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.
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/Sample.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE BangPatterns #-} -- | -- Module : Statistics.Sample -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan@@ -452,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/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@@ -44,7 +44,10 @@ | 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.@@ -71,5 +74,8 @@ | 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/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
+ 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,75 +0,0 @@-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]
changelog.md view
@@ -1,3 +1,13 @@+## 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`,@@ -61,11 +71,11 @@ * 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` @@ -130,7 +140,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,
statistics.cabal view
@@ -2,7 +2,7 @@ build-type: Simple name: statistics-version: 0.16.3.0+version: 0.16.4.0 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful@@ -36,7 +36,6 @@ extra-source-files: README.markdown- changelog.md examples/kde/KDE.hs examples/kde/data/faithful.csv examples/kde/kde.html@@ -44,6 +43,9 @@ tests/utils/Makefile tests/utils/fftw.c +extra-doc-files:+ changelog.md+ tested-with: GHC ==8.4.4 || ==8.6.5@@ -52,9 +54,10 @@ || ==9.0.2 || ==9.2.8 || ==9.4.8- || ==9.6.6+ || ==9.6.7 || ==9.8.4- || ==9.10.1+ || ==9.10.2+ || ==9.12.2 source-repository head type: git@@ -105,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@@ -132,7 +137,7 @@ , binary >= 0.5.1.0 , primitive >= 0.3 , dense-linear-algebra >= 0.1 && <0.2- , parallel >= 3.2.2.0 && <3.3+ , parallel >= 3.2.2.0 && <3.4 , vector >= 0.10 , vector-algorithms >= 0.4 , vector-th-unbox@@ -169,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@@ -201,7 +208,7 @@ build-depends: base -any , statistics -any- , doctest >=0.15 && <0.24+ , 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@@ -219,7 +226,7 @@ import: bench-stanza type: exitcode-stdio-1.0 hs-source-dirs: benchmark bench-time- main-is: bench.hs+ main-is: Main.hs Other-modules: Bench build-depends: tasty-bench >= 0.3 @@ -229,6 +236,6 @@ if impl(ghcjs) || !flag(BenchPAPI) buildable: False hs-source-dirs: benchmark bench-papi- main-is: bench.hs+ main-is: Main.hs Other-modules: Bench build-depends: tasty-papi >= 0.1.2
tests/Tests/Distribution.hs view
@@ -340,7 +340,7 @@ 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
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)