statistics 0.13.2.3 → 0.13.3.0
raw patch · 17 files changed
+637/−55 lines, 17 files
Files
- Statistics/Correlation.hs +70/−0
- Statistics/Correlation/Kendall.hs +4/−4
- Statistics/Distribution.hs +2/−2
- Statistics/Distribution/Exponential.hs +1/−1
- Statistics/Distribution/Laplace.hs +125/−0
- Statistics/Matrix.hs +120/−3
- Statistics/Matrix/Mutable.hs +12/−0
- Statistics/Sample.hs +50/−0
- Statistics/Sample/KernelDensity.hs +21/−9
- Statistics/Test/Internal.hs +43/−3
- Statistics/Transform.hs +39/−13
- benchmark/bench.hs +5/−0
- examples/kde/KDE.hs +4/−3
- statistics.cabal +4/−2
- tests/Tests/Correlation.hs +131/−13
- tests/Tests/Distribution.hs +4/−0
- tests/Tests/Transform.hs +2/−2
+ Statistics/Correlation.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE BangPatterns #-}+-- |+-- Module : Statistics.Correlation.Pearson+--+module Statistics.Correlation+ ( -- * Pearson correlation+ pearson+ , pearsonMatByRow+ -- * Spearman correlation+ , spearman+ , spearmanMatByRow+ ) where++import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import Statistics.Matrix+import Statistics.Sample+import Statistics.Test.Internal (rankUnsorted)+++----------------------------------------------------------------+-- Pearson+----------------------------------------------------------------++-- | Pearson correlation for sample of pairs.+pearson :: (G.Vector v (Double, Double), G.Vector v Double)+ => v (Double, Double) -> Double+pearson = correlation+{-# INLINE pearson #-}++-- | Compute pairwise pearson correlation between rows of a matrix+pearsonMatByRow :: Matrix -> Matrix+pearsonMatByRow m+ = generateSym (rows m)+ (\i j -> pearson $ row m i `U.zip` row m j)+{-# INLINE pearsonMatByRow #-}++++----------------------------------------------------------------+-- Spearman+----------------------------------------------------------------++-- | 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)+ )+ => v (a, b)+ -> Double+spearman xy+ = pearson+ $ G.zip (rankUnsorted x) (rankUnsorted y)+ where+ (x, y) = G.unzip xy+{-# INLINE spearman #-}++-- | compute pairwise spearman correlation between rows of a matrix+spearmanMatByRow :: Matrix -> Matrix+spearmanMatByRow+ = pearsonMatByRow . fromRows . fmap rankUnsorted . toRows+{-# INLINE spearmanMatByRow #-}
Statistics/Correlation/Kendall.hs view
@@ -8,11 +8,11 @@ -- This module implementes Kendall's tau form b which allows ties in the data. -- This is the same formula used by other statistical packages, e.g., R, matlab. ----- $$\tau = \frac{n_c - n_d}{\sqrt{(n_0 - n_1)(n_0 - n_2)}}$$+-- > \tau = \frac{n_c - n_d}{\sqrt{(n_0 - n_1)(n_0 - n_2)}} ----- where $n_0 = n(n-1)/2$, $n_1 = number of pairs tied for the first quantify$,--- $n_2 = number of pairs tied for the second quantify$,--- $n_c = number of concordant pairs$, $n_d = number of discordant pairs$.+-- where n_0 = n(n-1)\/2, n_1 = number of pairs tied for the first quantify,+-- n_2 = number of pairs tied for the second quantify,+-- n_c = number of concordant pairs$, n_d = number of discordant pairs. module Statistics.Correlation.Kendall ( kendall
Statistics/Distribution.hs view
@@ -8,7 +8,7 @@ -- Stability : experimental -- Portability : portable ----- Types classes for probability distrubutions+-- Type classes for probability distributions module Statistics.Distribution (@@ -57,7 +57,7 @@ -- -- > complCumulative d x = 1 - cumulative d x --- -- It's useful when one is interested in P(/X/</x/) and+ -- It's useful when one is interested in P(/X/>/x/) and -- expression on the right side begin to lose precision. This -- function have default implementation but implementors are -- encouraged to provide more precise implementation.
Statistics/Distribution/Exponential.hs view
@@ -98,7 +98,7 @@ error $ "Statistics.Distribution.Exponential.quantile: p must be in [0,1] range. Got: "++show p -- | Create an exponential distribution.-exponential :: Double -- ^ λ (scale) parameter.+exponential :: Double -- ^ Rate parameter. -> ExponentialDistribution exponential l | l <= 0 =
+ Statistics/Distribution/Laplace.hs view
@@ -0,0 +1,125 @@+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}+-- |+-- Module : Statistics.Distribution.Laplace+-- Copyright : (c) 2015 Mihai Maruseac+-- License : BSD3+--+-- Maintainer : mihai.maruseac@maruseac.com+-- Stability : experimental+-- Portability : portable+--+-- The Laplace distribution. This is the continuous probability+-- defined as the difference of two iid exponential random variables+-- or a Brownian motion evaluated as exponentially distributed times.+-- It is used in differential privacy (Laplace Method), speech+-- recognition and least absolute deviations method (Laplace's first+-- law of errors, giving a robust regression method)+--++module Statistics.Distribution.Laplace+ (+ LaplaceDistribution+ -- * Constructors+ , laplace+ , laplaceFromSample+ -- * Accessors+ , ldLocation+ , ldScale+ ) where++import Data.Aeson (FromJSON, ToJSON)+import Data.Binary (Binary(..))+import Data.Data (Data, Typeable)+import GHC.Generics (Generic)+import qualified Data.Vector.Generic as G+import qualified Statistics.Distribution as D+import qualified Statistics.Quantile as Q+import qualified Statistics.Sample as S+import Statistics.Types (Sample)+import Control.Applicative ((<$>), (<*>))+++data LaplaceDistribution = LD {+ ldLocation :: {-# UNPACK #-} !Double+ -- ^ Location.+ , ldScale :: {-# UNPACK #-} !Double+ -- ^ Scale.+ } deriving (Eq, Read, Show, Typeable, Data, Generic)++instance FromJSON LaplaceDistribution+instance ToJSON LaplaceDistribution++instance Binary LaplaceDistribution where+ put (LD l s) = put l >> put s+ get = LD <$> get <*> get++instance D.Distribution LaplaceDistribution where+ cumulative = cumulative+ complCumulative = complCumulative++instance D.ContDistr LaplaceDistribution where+ density (LD l s) x = exp (- abs (x - l) / s) / (2 * s)+ logDensity (LD l s) x = - abs (x - l) / s - log 2 - log s+ quantile = quantile++instance D.Mean LaplaceDistribution where+ mean (LD l _) = l++instance D.Variance LaplaceDistribution where+ variance (LD _ s) = 2 * s * s++instance D.MaybeMean LaplaceDistribution where+ maybeMean = Just . D.mean++instance D.MaybeVariance LaplaceDistribution where+ maybeStdDev = Just . D.stdDev+ maybeVariance = Just . D.variance++instance D.Entropy LaplaceDistribution where+ entropy (LD _ s) = 1 + log (2 * s)++instance D.MaybeEntropy LaplaceDistribution where+ maybeEntropy = Just . D.entropy++instance D.ContGen LaplaceDistribution where+ genContVar = D.genContinous++cumulative :: LaplaceDistribution -> Double -> Double+cumulative (LD l s) x+ | x <= l = 0.5 * exp ( (x - l) / s)+ | otherwise = 1 - 0.5 * exp ( - (x - l) / s )++complCumulative :: LaplaceDistribution -> Double -> Double+complCumulative (LD l s) x+ | x <= l = 1 - 0.5 * exp ( (x - l) / s)+ | otherwise = 0.5 * exp ( - (x - l) / s )++quantile :: LaplaceDistribution -> Double -> Double+quantile (LD l s) p+ | p == 0 = -inf+ | p == 1 = inf+ | p == 0.5 = l+ | p > 0 && p < 0.5 = l + s * log (2 * p)+ | p > 0.5 && p < 1 = l - s * log (2 - 2 * p)+ | otherwise =+ error $ "Statistics.Distribution.Laplace.quantile: p must be in [0,1] range. Got: "++show p+ where+ inf = 1 / 0++-- | Create an Laplace distribution.+laplace :: Double -- ^ Location+ -> Double -- ^ Scale+ -> LaplaceDistribution+laplace l s+ | s <= 0 =+ error $ "Statistics.Distribution.Laplace.laplace: scale parameter must be positive. Got " ++ show s+ | otherwise = LD l 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.+laplaceFromSample :: Sample -> LaplaceDistribution+laplaceFromSample xs = LD s l+ where+ s = Q.continuousBy Q.medianUnbiased 1 2 xs+ l = S.mean $ G.map (\x -> abs $ x - s) xs
Statistics/Matrix.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE PatternGuards #-} -- | -- Module : Statistics.Matrix -- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan@@ -9,13 +10,25 @@ -- we implement the necessary minimum here. module Statistics.Matrix- (+ ( -- * Data types Matrix(..) , Vector- , fromList+ -- * Conversion from/to lists/vectors , fromVector+ , fromList+ , fromRowLists+ , fromRows+ , fromColumns , toVector , toList+ , toRows+ , toColumns+ , toRowLists+ -- * Other+ , generate+ , generateSym+ , ident+ , diag , dimension , center , multiply@@ -34,11 +47,22 @@ ) where import Prelude hiding (exponent, map, sum)+import Control.Applicative ((<$>))+import Control.Monad.ST+import qualified Data.Vector.Unboxed as U+import Data.Vector.Unboxed ((!))+import qualified Data.Vector.Unboxed.Mutable as UM+ import Statistics.Function (for, square) import Statistics.Matrix.Types+import Statistics.Matrix.Mutable (unsafeNew,unsafeWrite,unsafeFreeze) import Statistics.Sample.Internal (sum)-import qualified Data.Vector.Unboxed as U ++----------------------------------------------------------------+-- Conversion to/from vectors/lists+----------------------------------------------------------------+ -- | Convert from a row-major list. fromList :: Int -- ^ Number of rows. -> Int -- ^ Number of columns.@@ -46,6 +70,10 @@ -> Matrix fromList r c = fromVector r c . U.fromList +-- | create a matrix from a list of lists, as rows+fromRowLists :: [[Double]] -> Matrix+fromRowLists = fromRows . fmap U.fromList+ -- | Convert from a row-major vector. fromVector :: Int -- ^ Number of rows. -> Int -- ^ Number of columns.@@ -56,6 +84,22 @@ | otherwise = Matrix r c 0 v where len = U.length v +-- | create a matrix from a list of vectors, as rows+fromRows :: [Vector] -> Matrix+fromRows xs+ | [] <- xs = error "Statistics.Matrix.fromRows: empty list of rows!"+ | any (/=nCol) ns = error "Statistics.Matrix.fromRows: row sizes do not match"+ | nCol == 0 = error "Statistics.Matrix.fromRows: zero columns in matrix"+ | otherwise = fromVector nRow nCol (U.concat xs)+ where+ nCol:ns = U.length <$> xs+ nRow = length xs+++-- | create a matrix from a list of vectors, as columns+fromColumns :: [Vector] -> Matrix+fromColumns = transpose . fromRows+ -- | Convert to a row-major flat vector. toVector :: Matrix -> U.Vector Double toVector (Matrix _ _ _ v) = v@@ -64,6 +108,78 @@ toList :: Matrix -> [Double] toList = U.toList . toVector +-- | Convert to a list of lists, as rows+toRowLists :: Matrix -> [[Double]]+toRowLists (Matrix _ nCol _ v)+ = chunks $ U.toList v+ where+ chunks [] = []+ chunks xs = case splitAt nCol xs of+ (rowE,rest) -> rowE : chunks rest+++-- | Convert to a list of vectors, as rows+toRows :: Matrix -> [Vector]+toRows (Matrix _ nCol _ v) = chunks v+ where+ chunks xs+ | U.null xs = []+ | otherwise = case U.splitAt nCol xs of+ (rowE,rest) -> rowE : chunks rest++-- | Convert to a list of vectors, as columns+toColumns :: Matrix -> [Vector]+toColumns = toRows . transpose++++----------------------------------------------------------------+-- Other+----------------------------------------------------------------++-- | Generate matrix using function+generate :: Int -- ^ Number of rows+ -> Int -- ^ Number of columns+ -> (Int -> Int -> Double)+ -- ^ Function which takes /row/ and /column/ as argument.+ -> Matrix+generate nRow nCol f+ = Matrix nRow nCol 0 $ U.generate (nRow*nCol) $ \i ->+ let (r,c) = i `quotRem` nCol in f r c++-- | Generate symmetric square matrix using function+generateSym+ :: Int -- ^ Number of rows and columns+ -> (Int -> Int -> Double)+ -- ^ Function which takes /row/ and /column/ as argument. It must+ -- be symmetric in arguments: @f i j == f j i@+ -> Matrix+generateSym n f = runST $ do+ m <- unsafeNew n n+ for 0 n $ \r -> do+ unsafeWrite m r r (f r r)+ for (r+1) n $ \c -> do+ let x = f r c+ unsafeWrite m r c x+ unsafeWrite m c r x+ unsafeFreeze m+++-- | Create the square identity matrix with given dimensions.+ident :: Int -> Matrix+ident n = diag $ U.replicate n 1.0++-- | Create a square matrix with given diagonal, other entries default to 0+diag :: Vector -> Matrix+diag v+ = Matrix n n 0 $ U.create $ do+ arr <- UM.replicate (n*n) 0+ for 0 n $ \i ->+ UM.unsafeWrite arr (i*n + i) (v ! i)+ return arr+ where+ n = U.length v+ -- | Return the dimensions of this matrix, as a (row,column) pair. dimension :: Matrix -> (Int, Int) dimension (Matrix r c _ _) = (r, c)@@ -125,6 +241,7 @@ -> Double unsafeIndex = unsafeBounds U.unsafeIndex +-- | Apply function to every element of matrix map :: (Double -> Double) -> Matrix -> Matrix map f (Matrix r c e v) = Matrix r c e (U.map f v)
Statistics/Matrix/Mutable.hs view
@@ -12,6 +12,7 @@ , replicate , thaw , bounds+ , unsafeNew , unsafeFreeze , unsafeRead , unsafeWrite@@ -36,6 +37,17 @@ unsafeFreeze :: MMatrix s -> ST s Matrix unsafeFreeze (MMatrix r c e mv) = Matrix r c e <$> U.unsafeFreeze mv++-- | Allocate new matrix. Matrix content is not initialized hence unsafe.+unsafeNew :: Int -- ^ Number of row+ -> Int -- ^ Number of columns+ -> ST s (MMatrix s)+unsafeNew r c+ | r < 0 = error "Statistics.Matrix.Mutable.unsafeNew: negative number of rows"+ | c < 0 = error "Statistics.Matrix.Mutable.unsafeNew: negative number of columns"+ | otherwise = do+ vec <- M.new (r*c)+ return $ MMatrix r c 0 vec unsafeRead :: MMatrix s -> Int -> Int -> ST s Double unsafeRead mat r c = unsafeBounds mat r c M.unsafeRead
Statistics/Sample.hs view
@@ -50,6 +50,10 @@ , fastVarianceUnbiased , fastStdDev + -- * Joint distirbutions+ , covariance+ , correlation+ , pair -- * References -- $references ) where@@ -339,6 +343,52 @@ fastStdDev :: (G.Vector v Double) => v Double -> Double fastStdDev = sqrt . fastVariance {-# INLINE fastStdDev #-}++-- | Covariance of sample of pairs. For empty sample it's set to+-- zero+covariance :: (G.Vector v (Double,Double), G.Vector v 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)+ where+ n = G.length xy+ (xs,ys) = G.unzip xy+ muX = mean xs+ muY = mean ys+{-# 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)+ => 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)+{-# SPECIALIZE correlation :: U.Vector (Double,Double) -> Double #-}+{-# SPECIALIZE correlation :: V.Vector (Double,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)+pair va vb+ | G.length va == G.length vb = G.zip va vb+ | otherwise = error "Statistics.Sample.pair: vector must have same length"+{-# INLINE pair #-} ------------------------------------------------------------------------ -- Helper code. Monomorphic unpacked accumulators.
Statistics/Sample/KernelDensity.hs view
@@ -31,9 +31,10 @@ import Statistics.Math.RootFinding (fromRoot, ridders) import Statistics.Sample.Histogram (histogram_) import Statistics.Sample.Internal (sum)-import Statistics.Transform (dct, idct)-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Unboxed as U+import Statistics.Transform (CD, dct, idct)+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector as V -- | Gaussian kernel density estimator for one-dimensional data, using@@ -46,18 +47,23 @@ -- mesh interval, use 'kde_'.) -- -- * Density estimates at each mesh point.-kde :: Int+kde :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)+ => Int -- ^ The number of mesh points to use in the uniform discretization -- of the interval @(min,max)@. If this value is not a power of -- two, then it is rounded up to the next power of two.- -> U.Vector Double -> (U.Vector Double, U.Vector Double)+ -> v Double -> (v Double, v Double) kde n0 xs = kde_ n0 (lo - range / 10) (hi + range / 10) xs where (lo,hi) = minMax xs- range | U.length xs <= 1 = 1 -- Unreasonable guess+ range | G.length xs <= 1 = 1 -- Unreasonable guess | lo == hi = 1 -- All elements are equal | otherwise = hi - lo+{-# INLINABLE kde #-}+{-# SPECIAlIZE kde :: Int -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}+{-# SPECIAlIZE kde :: Int -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-} + -- | Gaussian kernel density estimator for one-dimensional data, using -- the method of Botev et al. --@@ -66,7 +72,8 @@ -- * The coordinates of each mesh point. -- -- * Density estimates at each mesh point.-kde_ :: Int+kde_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)+ => Int -- ^ The number of mesh points to use in the uniform discretization -- of the interval @(min,max)@. If this value is not a power of -- two, then it is rounded up to the next power of two.@@ -74,9 +81,10 @@ -- ^ Lower bound (@min@) of the mesh range. -> Double -- ^ Upper bound (@max@) of the mesh range.- -> U.Vector Double -> (U.Vector Double, U.Vector Double)+ -> v Double+ -> (v Double, v Double) kde_ n0 min max xs- | U.null xs = error "Statistics.KernelDensity.kde: empty sample"+ | G.null xs = error "Statistics.KernelDensity.kde: empty sample" | n0 <= 1 = error "Statistics.KernelDensity.kde: invalid number of points" | otherwise = (mesh, density) where@@ -103,6 +111,10 @@ const = (1 + 0.5 ** (s+0.5)) / 3 k0 = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi sqr x = x * x+{-# INLINABLE kde_ #-}+{-# SPECIAlIZE kde_ :: Int -> Double -> Double -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}+{-# SPECIAlIZE kde_ :: Int -> Double -> Double -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}+ -- $references --
Statistics/Test/Internal.hs view
@@ -1,11 +1,15 @@ {-# LANGUAGE FlexibleContexts #-} module Statistics.Test.Internal ( rank+ , rankUnsorted , splitByTags ) where -import qualified Data.Vector.Generic as G-+import Data.Ord+import Data.Vector.Generic ((!))+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as M+import Statistics.Function -- Private data type for unfolding@@ -16,7 +20,14 @@ , rankVec :: v a -- Remaining vector } --- | Calculate rank of sample. Sample should be already sorted.+-- | Calculate rank of every element of sample. In case of ties ranks+-- are averaged. Sample should be already sorted in ascending order.+--+-- >>> rank (==) (fromList [10,20,30::Int])+-- > fromList [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) => (a -> a -> Bool) -- ^ Equivalence relation -> v a -- ^ Vector to rank@@ -37,6 +48,35 @@ (h,rest) = G.span (eq $ G.head v) v go (Rank n val r v) = Just (val, Rank (n-1) val r v) {-# INLINE rank #-}++-- | Compute rank of every element of vector. Unlike rank it doesn't+-- require sample to be sorted.+rankUnsorted :: ( Ord a+ , G.Vector v a+ , G.Vector v Int+ , G.Vector v Double+ , G.Vector v (Int, a)+ )+ => v a+ -> v Double+rankUnsorted xs = G.create $ do+ -- Put ranks into their original positions+ -- NOTE: backpermute will do wrong thing+ vec <- M.new n+ for 0 n $ \i ->+ M.unsafeWrite vec (index ! i) (ranks ! i)+ return vec+ where+ n = G.length xs+ -- Calculate ranks for sorted array+ ranks = rank (==) sorted+ -- Sort vector and retain original indices of elements+ (index, sorted)+ = G.unzip+ $ sortBy (comparing snd)+ $ indexed xs+{-# INLINE rankUnsorted #-}+ -- | Split tagged vector splitByTags :: (G.Vector v a, G.Vector v (Bool,a)) => v (Bool,a) -> (v a, v a)
Statistics/Transform.hs view
@@ -34,23 +34,30 @@ import Data.Bits (shiftL, shiftR) import Data.Complex (Complex(..), conjugate, realPart) import Numeric.SpecFunctions (log2)-import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as M-import qualified Data.Vector.Unboxed as U-+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector as V type CD = Complex Double -- | Discrete cosine transform (DCT-II).-dct :: U.Vector Double -> U.Vector Double+dct :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v Double -> v Double dct = dctWorker . G.map (:+0)+{-# INLINABLE dct #-}+{-# SPECIAlIZE dct :: U.Vector Double -> U.Vector Double #-}+{-# SPECIAlIZE dct :: V.Vector Double -> V.Vector Double #-} -- | Discrete cosine transform (DCT-II). Only real part of vector is -- transformed, imaginary part is ignored.-dct_ :: U.Vector CD -> U.Vector Double+dct_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v CD -> v Double dct_ = dctWorker . G.map (\(i :+ _) -> i :+ 0)+{-# INLINABLE dct_ #-}+{-# SPECIAlIZE dct_ :: U.Vector CD -> U.Vector Double #-}+{-# SPECIAlIZE dct_ :: V.Vector CD -> V.Vector Double#-} -dctWorker :: U.Vector CD -> U.Vector Double+dctWorker :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v CD -> v Double+{-# INLINE dctWorker #-} dctWorker xs -- length 1 is special cased because shuffle algorithms fail for it. | G.length xs == 1 = G.map ((2*) . realPart) xs@@ -70,15 +77,22 @@ -- 'dct' only up to scale parameter: -- -- > (idct . dct) x = (* length x)-idct :: U.Vector Double -> U.Vector Double+idct :: (G.Vector v CD, G.Vector v Double) => v Double -> v Double idct = idctWorker . G.map (:+0)+{-# INLINABLE idct #-}+{-# SPECIAlIZE idct :: U.Vector Double -> U.Vector Double #-}+{-# SPECIAlIZE idct :: V.Vector Double -> V.Vector Double #-} -- | Inverse discrete cosine transform (DCT-III). Only real part of vector is -- transformed, imaginary part is ignored.-idct_ :: U.Vector CD -> U.Vector Double+idct_ :: (G.Vector v CD, G.Vector v Double) => v CD -> v Double idct_ = idctWorker . G.map (\(i :+ _) -> i :+ 0)+{-# INLINABLE idct_ #-}+{-# SPECIAlIZE idct_ :: U.Vector CD -> U.Vector Double #-}+{-# SPECIAlIZE idct_ :: V.Vector CD -> V.Vector Double #-} -idctWorker :: U.Vector CD -> U.Vector Double+idctWorker :: (G.Vector v CD, G.Vector v Double) => v CD -> v Double+{-# INLINE idctWorker #-} idctWorker xs | vectorOK xs = G.generate len interleave | otherwise = error "Statistics.Transform.dct: bad vector length"@@ -93,21 +107,29 @@ len = G.length xs + -- | Inverse fast Fourier transform.-ifft :: U.Vector CD -> U.Vector CD+ifft :: G.Vector v CD => v CD -> v CD ifft xs | vectorOK xs = G.map ((/fi (G.length xs)) . conjugate) . fft . G.map conjugate $ xs | otherwise = error "Statistics.Transform.ifft: bad vector length"+{-# INLINABLE ifft #-}+{-# SPECIAlIZE ifft :: U.Vector CD -> U.Vector CD #-}+{-# SPECIAlIZE ifft :: V.Vector CD -> V.Vector CD #-} -- | Radix-2 decimation-in-time fast Fourier transform.-fft :: U.Vector CD -> U.Vector CD+fft :: G.Vector v CD => v CD -> v CD fft v | vectorOK v = G.create $ do mv <- G.thaw v mfft mv return mv | otherwise = error "Statistics.Transform.fft: bad vector length"+{-# INLINABLE fft #-}+{-# SPECIAlIZE fft :: U.Vector CD -> U.Vector CD #-}+{-# SPECIAlIZE fft :: V.Vector CD -> V.Vector CD #-} -- Vector length must be power of two. It's not checked mfft :: (M.MVector v CD) => v s CD -> ST s ()+{-# INLINE mfft #-} mfft vec = bitReverse 0 0 where bitReverse i j | i == len-1 = stage 0 1@@ -138,13 +160,17 @@ len = M.length vec m = log2 len ++----------------------------------------------------------------+-- Helpers+----------------------------------------------------------------+ fi :: Int -> CD fi = fromIntegral halve :: Int -> Int halve = (`shiftR` 1) --vectorOK :: U.Unbox a => U.Vector a -> Bool+vectorOK :: G.Vector v a => v a -> Bool {-# INLINE vectorOK #-} vectorOK v = (1 `shiftL` log2 n) == n where n = G.length v
benchmark/bench.hs view
@@ -3,6 +3,7 @@ import Data.Complex import Statistics.Sample import Statistics.Transform+import Statistics.Correlation.Pearson import System.Random.MWC import qualified Data.Vector.Unboxed as U @@ -35,6 +36,10 @@ , 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
examples/kde/KDE.hs view
@@ -4,11 +4,12 @@ import Statistics.Sample.KernelDensity (kde) import Text.Hastache (MuType(..), defaultConfig, hastacheFile) import Text.Hastache.Context (mkStrContext)-import qualified Data.Attoparsec as B-import qualified Data.Attoparsec.Char8 as A+import qualified Data.Attoparsec.ByteString as B+import qualified Data.Attoparsec.ByteString.Char8 as A import qualified Data.ByteString as B import qualified Data.ByteString.Lazy as L import qualified Data.Vector.Unboxed as U+import qualified Data.Text.Lazy.IO as TL csv = do B.takeTill A.isEndOfLine@@ -20,4 +21,4 @@ let xs = map (\(a,b) -> [a,b]) . U.toList . uncurry U.zip . kde 64 $ waits context "data" = MuVariable . show $ xs s <- hastacheFile defaultConfig "kde.tpl" (mkStrContext context)- L.writeFile "kde.html" s+ TL.writeFile "kde.html" s
statistics.cabal view
@@ -1,5 +1,5 @@ name: statistics-version: 0.13.2.3+version: 0.13.3.0 synopsis: A library of statistical types, data, and functions description: This library provides a number of common functions and types useful@@ -49,6 +49,7 @@ exposed-modules: Statistics.Autocorrelation Statistics.Constants+ Statistics.Correlation Statistics.Correlation.Kendall Statistics.Distribution Statistics.Distribution.Beta@@ -60,6 +61,7 @@ Statistics.Distribution.Gamma Statistics.Distribution.Geometric Statistics.Distribution.Hypergeometric+ Statistics.Distribution.Laplace Statistics.Distribution.Normal Statistics.Distribution.Poisson Statistics.Distribution.StudentT@@ -134,7 +136,7 @@ Tests.Transform ghc-options:- -Wall -threaded -rtsopts+ -Wall -threaded -rtsopts -fsimpl-tick-factor=500 build-depends: HUnit,
tests/Tests/Correlation.hs view
@@ -3,35 +3,153 @@ module Tests.Correlation ( tests ) where +import Control.Arrow (Arrow(..))+import qualified Data.Vector as V+import Statistics.Matrix hiding (map)+import Statistics.Correlation+import Statistics.Correlation.Kendall+import Test.QuickCheck ((==>),Property,counterexample) import Test.Framework import Test.Framework.Providers.QuickCheck2 import Test.Framework.Providers.HUnit-import Test.HUnit (Assertion, (@=?))-import qualified Data.Vector as V-import Statistics.Correlation.Kendall+import Test.HUnit (Assertion, (@=?), assertBool) +import Tests.ApproxEq++----------------------------------------------------------------+-- Tests list+----------------------------------------------------------------+ tests :: Test tests = testGroup "Correlation"- [ testProperty "Kendall test -- general" testKendall- , testCase "Kendall test -- special cases" testKendallSpecial+ [ testProperty "Pearson correlation" testPearson+ , testProperty "Spearman correlation is scale invariant" testSpearmanScale+ , testProperty "Spearman correlation, nonlinear" testSpearmanNonlinear+ , testProperty "Kendall test -- general" testKendall+ , testCase "Kendall test -- special cases" testKendallSpecial ] ++----------------------------------------------------------------+-- Pearson's correlation+----------------------------------------------------------------++testPearson :: [(Double,Double)] -> Property+testPearson sample+ = (length sample > 1) ==> (exact ~= fast)+ where+ (~=) = eql 1e-12+ exact = exactPearson $ map (realToFrac *** realToFrac) sample+ fast = pearson $ V.fromList sample++exactPearson :: [(Rational,Rational)] -> Double+exactPearson sample+ = realToFrac cov / sqrt (realToFrac (varX * varY))+ where+ (xs,ys) = unzip sample+ n = fromIntegral $ length sample+ -- Mean+ muX = sum xs / n+ muY = sum ys / n+ -- Mean of squares+ muX2 = sum (map (\x->x*x) xs) / n+ muY2 = sum (map (\x->x*x) ys) / n+ -- Covariance+ cov = sum (zipWith (*) [x - muX | x<-xs] [y - muY | y<-ys]) / n+ varX = muX2 - muX*muX+ varY = muY2 - muY*muY++----------------------------------------------------------------+-- Spearman's correlation+----------------------------------------------------------------++-- Test that Spearman correlation is scale invariant+testSpearmanScale :: [(Double,Double)] -> Double -> Property+testSpearmanScale xs a+ = and [ length xs > 1 -- Enough to calculate underflow+ , a /= 0+ , not (isNaN c1)+ , not (isNaN c2)+ , not (isNaN c3)+ , not (isNaN c4)+ ]+ ==> ( counterexample (show xs2)+ $ counterexample (show xs3)+ $ counterexample (show xs4)+ $ counterexample (show (c1,c2,c3,c4))+ $ and [ c1 == c4+ , c1 == signum a * c2+ , c1 == signum a * c3+ ]+ )+ where+ xs2 = map ((*a) *** id ) xs+ xs3 = map (id *** (*a)) xs+ xs4 = map ((*a) *** (*a)) xs+ c1 = spearman $ V.fromList xs+ c2 = spearman $ V.fromList xs2+ c3 = spearman $ V.fromList xs3+ c4 = spearman $ V.fromList xs4++-- Test that Spearman correlation allows to transform sample with+testSpearmanNonlinear :: [(Double,Double)] -> Property+testSpearmanNonlinear sample0+ = and [ length sample0 > 1+ , not (isNaN c1)+ , not (isNaN c2)+ , not (isNaN c3)+ , not (isNaN c4)+ ]+ ==> ( counterexample (show sample0)+ $ counterexample (show sample1)+ $ counterexample (show sample2)+ $ counterexample (show sample3)+ $ counterexample (show sample4)+ $ counterexample (show (c1,c2,c3,c4))+ $ and [ c1 == c2+ , c1 == c3+ , c1 == c4+ ]+ )+ where+ -- 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 ]+ where+ a = minimum xs+ b = maximum xs+ sample1 = uncurry zip $ (stretch *** stretch) $ unzip sample0+ sample2 = map (exp *** id ) sample1+ sample3 = map (id *** exp) sample1+ sample4 = map (exp *** exp) sample1+ c1 = spearman $ V.fromList sample1+ c2 = spearman $ V.fromList sample2+ c3 = spearman $ V.fromList sample3+ c4 = spearman $ V.fromList sample4+++----------------------------------------------------------------+-- Kendall's correlation+----------------------------------------------------------------+ testKendall :: [(Double, Double)] -> Bool testKendall xy | isNaN r1 = isNaN r2 | otherwise = r1 == r2 where- r1 = kendallBruteForce xy + r1 = kendallBruteForce xy r2 = kendall $ V.fromList xy testKendallSpecial :: Assertion-testKendallSpecial = ys @=? map (kendall.V.fromList) xs- where - (xs, ys) = unzip testData- testData :: [([(Double, Double)], Double)]- testData = [ ( [(1,1), (2,2), (3,1), (1,5), (2,2)], -0.375 )- , ( [(1,3), (1,3), (1,3), (3,2), (3,5)], 0)+testKendallSpecial = vs @=? map (\(xs, ys) -> kendall $ V.fromList $ zip xs ys) d+ where+ (d, vs) = unzip testData+ testData :: [(([Double], [Double]), Double)]+ testData = [ (([1, 2, 3, 1, 2], [1, 2, 1, 5, 2]), -0.375)+ , (([1, 1, 1, 3, 3], [3, 3, 3, 2, 5]), 0) ]- + kendallBruteForce :: [(Double, Double)] -> Double kendallBruteForce xy = (n_c - n_d) / sqrt ((n_0 - n_1) * (n_0 - n_2))
tests/Tests/Distribution.hs view
@@ -17,6 +17,7 @@ import Statistics.Distribution.Gamma (GammaDistribution, gammaDistr) import Statistics.Distribution.Geometric import Statistics.Distribution.Hypergeometric+import Statistics.Distribution.Laplace (LaplaceDistribution, laplace) import Statistics.Distribution.Normal (NormalDistribution, normalDistr) import Statistics.Distribution.Poisson (PoissonDistribution, poisson) import Statistics.Distribution.StudentT@@ -42,6 +43,7 @@ , contDistrTests (T :: T ChiSquared ) , contDistrTests (T :: T ExponentialDistribution ) , contDistrTests (T :: T GammaDistribution )+ , contDistrTests (T :: T LaplaceDistribution ) , contDistrTests (T :: T NormalDistribution ) , contDistrTests (T :: T UniformDistribution ) , contDistrTests (T :: T StudentT )@@ -252,6 +254,8 @@ arbitrary = binomial <$> QC.choose (1,100) <*> QC.choose (0,1) instance QC.Arbitrary ExponentialDistribution where arbitrary = exponential <$> QC.choose (0,100)+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) instance QC.Arbitrary BetaDistribution where
tests/Tests/Transform.hs view
@@ -60,7 +60,7 @@ -- vector should be replicated in every real component of the result, -- and all the imaginary components should be zero. t_impulse :: Double -> Positive Int -> Bool-t_impulse k (Positive m) = G.all (c_near i) (fft v)+t_impulse k (Positive m) = U.all (c_near i) (fft v) where v = i `G.cons` G.replicate (n-1) 0 i = k :+ 0 n = 1 `shiftL` (m .&. 6)@@ -69,7 +69,7 @@ -- otherwise zero vector, the sum-of-squares of each component of the -- result should equal the square of the impulse. t_impulse_offset :: Double -> Positive Int -> Positive Int -> Bool-t_impulse_offset k (Positive x) (Positive m) = G.all ok (fft v)+t_impulse_offset k (Positive x) (Positive m) = U.all ok (fft v) where v = G.concat [G.replicate xn 0, G.singleton i, G.replicate (n-xn-1) 0] ok (re :+ im) = within ulps (re*re + im*im) (k*k) i = k :+ 0