Merge branch 'oldmaster'

Commit which were missing after fixing problems with
HG<->git mirroring

Conflicts:
	.hgtags
	statistics.cabal

.hgtags

@@ -44,4 +44,5 @@ c75dd236003cf7eadab302d31f7523767407ea92 0.11.0.3
 2eed7bb5d817365a264ffc9b51f60e54d12cf328 0.13.1.1
 74572dc8487579c11f69400b91863eb1a150e507 0.13.2.0
 e9452e9e69cedccbcbd13e407c38cca6286ac66f 0.13.2.2
+b4307c59867021ffc58f5afa62121f99b8d352ff 0.13.2.1
 4302388f1acb398ec24f959175542be7addee36f 0.13.2.3

Statistics/Correlation.hs

@@ -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/Distribution/Exponential.hs

@@ -98,7 +98,7 @@ quantile (ED l) p
     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

@@ -0,0 +1,111 @@
+{-# 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
+    -- * Accessors
+    , ldLocation
+    , ldScale
+    ) where
+
+import Data.Aeson (FromJSON, ToJSON)
+import Data.Binary (Binary(..))
+import Data.Data (Data, Typeable)
+import GHC.Generics (Generic)
+import qualified Statistics.Distribution         as D
+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

Statistics/Matrix.hs

@@ -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,10 +47,21 @@ module Statistics.Matrix
     ) 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.
@@ -46,6 +70,10 @@ fromList :: Int                 -- ^ Number of rows.
          -> 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 @@ fromVector r c v
   | 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 @@ toVector (Matrix _ _ _ v) = v
 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 @@ unsafeIndex :: Matrix
             -> 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)