hstatistics 0.2.2.11 → 0.2.3

raw patch · 3 files changed

+225/−7 lines, 3 filesdep +vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: vector

API changes (from Hackage documentation)

Files

• CHANGES +3/−0
• hstatistics.cabal +5/−4
• lib/Numeric/Statistics.hs +217/−3

@@ -81,3 +81,6 @@  0.2.2.11: 		fixed bug in surrogate data sampling++0.2.3:+		added functions to Numeric.Statistics

@@ -1,8 +1,8 @@ Name:               hstatistics-Version:            0.2.2.11+Version:            0.2.3 License:            BSD3 License-file:       LICENSE-Copyright:          (c) A.V.H. McPhail 2010, 2011+Copyright:          (c) A.V.H. McPhail 2010, 2011, 2012 Author:             Vivian McPhail Maintainer:         haskell.vivian.mcphail <at> gmail <dot> com Stability:          provisional@@ -13,11 +13,11 @@      .      When hmatrix is installed with -fvector, the vector type is Data.Vector.Storable      from the vector package and compatible with the 'statistics' package -     <http://hackage.haskell.org/package/statistics +     <http://hackage.haskell.org/package/statistics>       .      Feature requests, suggestions, and bug fixes welcome. Category:           Math, Statistics-tested-with:        GHC ==7.0.1+tested-with:        GHC ==7.4.1  cabal-version:      >=1.8 @@ -30,6 +30,7 @@      Build-Depends:      base >= 4 && < 5,                         array, random,+                        vector,                         hmatrix >= 0.10.0.0,                         hmatrix-gsl-stats >= 0.1.2.9 

@@ -2,8 +2,8 @@ ----------------------------------------------------------------------------- -- | -- Module      :  Numeric.Statistics--- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010--- License     :  GPL-style+-- Copyright   :  (c) A. V. H. McPhail 2010, 2012+-- License     :  BSD -- -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com -- Stability   :  provisional@@ -15,22 +15,37 @@  module Numeric.Statistics (                            Sample,Samples-                          , covarianceMatrix+                          , covarianceMatrix, correlationCoefficientMatrix                           , meanList, meanArray, meanMatrix                           , varianceList, varianceArray, varianceMatrix+                          --+                          , centre, cloglog, corcoeff, cut+                          , ranks, kendall, logit+                          , mahalanobis+                          , mode, moment+                          , ols, percentile, range+                          , run_count+                          , spearman, studentize                           ) where   ----------------------------------------------------------------------------- +--import Debug.Trace+ --import Numeric.Vector --import Numeric.Matrix --import Numeric.Container import Numeric.LinearAlgebra  import qualified Data.Array.IArray as I +import qualified Data.List as DL+import qualified Data.Vector.Generic as GV +import Foreign.Storable+ import Numeric.GSL.Statistics+import Numeric.GSL.Sort(sort)  ----------------------------------------------------------------------------- @@ -45,7 +60,14 @@ covarianceMatrix d = let (s,f) = I.bounds d                       in fromArray2D $ I.array ((s,s),(f,f)) $ concat $ map (\(x,y) -> let c = covariance (d I.! x) (d I.! y) in if x == y then [((x,y),c)] else [((x,y),c),((y,x),c)]) $ filter (\(x,y) -> x <= y) $ I.range ((s,s),(f,f)) +-- | the correlation coefficient: (cov x y) / (std x) (std y)+correlationCoefficientMatrix :: Samples Double -> Matrix Double+correlationCoefficientMatrix d = let (s,f) = I.bounds d+                           in fromArray2D $ I.array ((s,s),(f,f)) $ concat $ map (\(x,y) -> let { x' = d I.! x ; y' = d I.! y ; c = (covariance x' y') / ((stddev x') * (stddev y')) } in if x == y then [((x,y),c)] else [((x,y),c),((y,x),c)]) $ filter (\(x,y) -> x <= y) $ I.range ((s,s),(f,f))+ -----------------------------------------------------------------------------+-----------------------------------------------------------------------------+-----------------------------------------------------------------------------  -- | the mean of a list of vectors meanList :: (Container Vector a, Num (Vector a)) => [Sample a] -> Sample a@@ -81,3 +103,195 @@ varianceMatrix a = varianceList $ toRows a  -----------------------------------------------------------------------------+-----------------------------------------------------------------------------+-----------------------------------------------------------------------------++-- | centre the data to 0: (x - (mean x))+centre :: Vector Double -> Vector Double+centre v = v - (realToFrac (mean v))++-----------------------------------------------------------------------------++-- | complementary log-log function+--cloglog :: Vector Double -> Vector Double+cloglog :: Floating a => a -> a+cloglog v = - log (- (log v))++-----------------------------------------------------------------------------++-- | corcoeff = covariance x / (std dev x * std dev y)+corcoeff :: Vector Double -> Vector Double -> Double+corcoeff x y = (covariance x y)/((stddev x)*(stddev y))++-----------------------------------------------------------------------------++-- | cut numerical data into intervals, data must fall inside the bounds+cut :: Vector Double +    -> Vector Double -- ^ intervals+    -> Vector Int    -- ^ data indexed by bin+cut v c  = let c' = sort c+           in mapVector (\x -> cut_helper 0 x c') v +    where+      cut_helper i x c +          | i >= dim c                       = error "Numeric.Statistics: cut: data point not within interval"+          | x >= (c @> i) && x <= (c @> (i+1)) = i+          | otherwise                       = cut_helper (i + 1) x c++-----------------------------------------------------------------------------++-- | return the rank of each element of the vector+--     multiple identical entries result in the average rank of those entries+--ranks :: Vector Double -> Vector Double+ranks :: (Fractional b, Storable b) => Vector Double -> Vector b+ranks v = let v' = sort v+          in mapVector (\x -> 1 + rank_helper x v') v+              where rank_helper x v' = let is = GV.elemIndices x v'+                                       in (realToFrac (GV.foldl (+) 0 is)) / (fromIntegral $ dim is)++-----------------------------------------------------------------------------++-- | kendall's rank correlation τ+kendall :: Vector Double -> Vector Double -> Matrix Double+kendall x y = let ln = dim x+                  rx = ranks x+                  ry = ranks y+                  r = fromColumns [rx,ry]+                  m = signum $ (kronecker r (asColumn $ constant 1.0 ln)) - (kronecker (asRow $ constant 1.0 ln) r)+                  c = rows m - 1+              in correlationCoefficientMatrix $ I.listArray (0,c) (toColumns m)++-----------------------------------------------------------------------------++-- | (logit p) = log(p/(1-p))+--logit :: Vector Double -> Vector Double+logit :: (Floating b, Storable b)+        => Vector b -> Vector b+logit v =  mapVector (\x -> - (log ((1 / x) - 1))) v++-----------------------------------------------------------------------------++-- | the Mahalanobis D-square distance between samples+--     columns are components and rows are observations (uses pseudoinverse)+mahalanobis :: Samples Double        -- ^ the data set+            -> Maybe (Sample Double) -- ^ (Just sample) to be measured or use mean when Nothing+            -> Double                -- ^ D^2 +mahalanobis x u = let (_,xr) = I.bounds x+                      xl     = I.elems x+                      s'     = pinv $ covarianceMatrix x+                      xu     = case u of+                                 Nothing -> fromList $ map mean xl+                                 Just m  -> m+                      xm     = fromRows $ map ((-) xu) $ toRows $ fromColumns xl+                      --um     = asColumn xu+                      --w      = ((trans xm) <> xm + (trans um) <> um)/(fromIntegral $ xr - 1)+                      --w'     = inv w+                  in ((xm <> s' <> (trans xm)) @@> (0,0)) ++-----------------------------------------------------------------------------++-- | a list of element frequencies+mode :: Vector Double -> [(Double,Integer)]+mode v = let w = sort v+         in DL.sortBy (\(_,n) (_,n') -> compare n' n) $ foldVector freqs [] w+            where freqs x []          = [(x,1)]+                  freqs x ((f,n):fns)+                      | f == x         = ((f,n+1):fns) +                      | otherwise     = ((x,1):(f,n):fns)++-----------------------------------------------------------------------------++-- | the p'th moment of a vector+moment :: Integral a +       => a             -- ^ moment+       -> Bool          -- ^ calculate central moment+       -> Bool          -- ^ calculate absolute moment+       -> Vector Double -- ^ data+       -> Double+moment p c a v +    | p <= 0     = error "Numeric.Statistics.moment: negative moment requested"+--    | p == 1     = mean v    +--    | p == 2     = variance v -- gives sample variance+    | otherwise = let u = if c then centre v else v+                      w = if a then abs u else u+                      x = mapVector (** (fromIntegral p)) w+                  in mean x++-----------------------------------------------------------------------------++-- | ordinary least squares estimation for the multivariate model+--   Y = X B + e        rows are observations, columns are elements+--   mean e = 0, cov e = kronecker s I+ols :: (Num (Vector t), Field t) +      => Matrix t         -- ^ X+    -> Matrix t           -- ^ Y+    -> (Matrix t, Matrix t, Matrix t) -- ^ (OLS estimator for B, OLS estimator for s, OLS residuals)+ols x y +    | rows x /= rows y = error "Numeric.Statistics: ols: incorrect matrix dimensions"+    | otherwise       = let (xr,xc) = (rows x,cols x)+                            (yr,yc) = (rows y,cols y)+                            z = (trans x) <> x+                            r = rank z+                            beta = if r == xc +                                      then (inv z) <> (trans x) <> y+                                      else (pinv x) <> y+                            rr = y - x <> beta+                            sigma = ((trans rr) <> rr) / (fromIntegral $ xr - r)+                        in (beta,rr,sigma)++-----------------------------------------------------------------------------++-- | compute quantiles in percent+percentile :: Double        -- ^ percentile (0 - 100)+           -> Vector Double -- ^ data+           -> Double        -- ^ result+percentile p d = quantile (0.01*p) d++-----------------------------------------------------------------------------++-- | the difference between the maximum and minimum of the input+range :: Container c e => c e -> e+range v = maxElement v - minElement v++-----------------------------------------------------------------------------++-- | count the number of runs greater than or equal to @n@ in the data+run_count :: (Num a, Num t, Ord b, Ord a, Storable b) +            => a             -- ^ longest run to count+          -> Vector b        -- ^ data+          -> [(a, t)]        -- ^ [(run length,count)]+run_count n v = let w = subVector 1 (dim v - 1) v+                    x = foldVector run_count' [(1,v @> 0)] w+                    y = map fst x+                    z = takeWhile (<= n) $  DL.sort y+                in foldr count [] z+    where run_count' m ((c,g):cs)+              | m < g             = ((c+1,m):cs)+              | otherwise         = ((1,m):(c,g):cs)+          count x []           = [(x,1)]+          count x ((yv,yc):ys)   +              | x == yv         = ((yv,yc+1):ys)+              | otherwise      = ((x,1):(yv,yc):ys)++-----------------------------------------------------------------------------++-- | Spearman's rank correlation coefficient+spearman :: Vector Double -> Vector Double -> Double+spearman x y = corcoeff (ranks x) (ranks y)++-----------------------------------------------------------------------------++-- | centre and normalise a vector+studentize :: Vector Double -> Vector Double+studentize x = (centre x)/(fromList $ [stddev x])++-----------------------------------------------------------------------------++--table++-----------------------------------------------------------------------------+++++-----------------------------------------------------------------------------+