hstatistics 0.2.0.5 → 0.2.0.6
raw patch · 8 files changed
+338/−34 lines, 8 filesdep +arraydep +randomdep ~hmatrix-gsl-stats
Dependencies added: array, random
Dependency ranges changed: hmatrix-gsl-stats
Files
- CHANGES +4/−0
- INSTALL +9/−27
- LICENSE +1/−1
- hstatistics.cabal +8/−4
- lib/Numeric/Statistics.hs +38/−0
- lib/Numeric/Statistics/ICA.hs +203/−0
- lib/Numeric/Statistics/Information.hs +2/−2
- lib/Numeric/Statistics/PCA.hs +73/−0
CHANGES view
@@ -31,3 +31,7 @@ 0.2.0.5: added Histogram++0.2.0.6:+ added PCA+ added ICA
INSTALL view
@@ -1,35 +1,17 @@ ------------------------------------------------ A simple signal processing library for Haskell+ A statistics library for Haskell ----------------------------------------------- INSTALLATION -Recommended method (ok in Ubuntu/Debian systems):- $ cabal install hsignal--INSTALLATION ON WINDOWS ------------------------------------------1) Install a recent ghc (e.g. ghc-6.10.3)--2) Install cabal-install. A binary for windows can be obtained from:-- http://www.haskell.org/cabal/release/cabal-install-0.6.2/cabal.exe-- Put it somewhere in the path, for instance in c:\ghc\ghc-6.10.3\bin--3) Download and uncompress hmatrix-x.y.z.tar.gz from Hackage:-- http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix--4) Open a terminal, cd to the hmatrix folder, and run-- > cabal install--5) Download and uncompress hsignal-x.y.z.tar.gz from Hackage:-- http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hsignal+cabal install hstatistics -6) Open a terminal, cd to the hsignal folder, and run+OR - > cabal install+tar xzf hstatistics-x.y.z.tar.gz+cd hstatistics+runhaskell Setup.lhs configure+runhaskell Setup.lhs build+runhaskell Setup.lhs hadock+runhaskell Setup.lhs install
LICENSE view
@@ -1,2 +1,2 @@-Copyright Alberto Ruiz 2006-2007+Copyright A.V.H. McPhail 2010 GPL license
hstatistics.cabal view
@@ -1,5 +1,5 @@ Name: hstatistics-Version: 0.2.0.5+Version: 0.2.0.6 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010@@ -9,7 +9,7 @@ Homepage: http://code.haskell.org/hstatistics Synopsis: Statistics Description: Purely functional interface for statistics based on hmatrix and hmatrix-gsl-stats-Category: Math+Category: Math, Statistics tested-with: GHC ==6.12.1 cabal-version: >=1.2@@ -22,12 +22,16 @@ library Build-Depends: base >= 3 && < 5,- hmatrix >= 0.10.0, hmatrix-gsl-stats >= 0.1.1.4+ array, random,+ hmatrix >= 0.10.0, hmatrix-gsl-stats >= 0.1.1.5 Extensions: hs-source-dirs: lib- Exposed-modules: Numeric.Statistics.Information+ Exposed-modules: Numeric.Statistics+ Numeric.Statistics.PCA+ Numeric.Statistics.ICA+ Numeric.Statistics.Information Numeric.Statistics.Histogram other-modules: C-sources:
+ lib/Numeric/Statistics.hs view
@@ -0,0 +1,38 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Statistics+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : portable+--+-- Useful statistical functions+--+-----------------------------------------------------------------------------++module Numeric.Statistics (+ covarianceMatrix+ ) where+++-----------------------------------------------------------------------------++import Data.Packed.Vector+import Data.Packed.Matrix++import qualified Data.Array.IArray as I ++import Numeric.GSL.Statistics++-----------------------------------------------------------------------------++-- | the covariance matrix+covarianceMatrix :: I.Array Int (Vector Double) -- ^ the dimensions of data (each vector being one dimension)+ -> Matrix Double -- ^ the symmetric covariance matrix+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))++-----------------------------------------------------------------------------
+ lib/Numeric/Statistics/ICA.hs view
@@ -0,0 +1,203 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Statistics.ICA+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : portable+--+-- Independent Components Analysis+--+-- implements the FastICA algorithm found in:+--+-- http://www.google.com/url?sa=t&source=web&cd=2&ved=0CBgQFjAB&url=http%3A%2F%2Fciteseerx.ist.psu.edu%2Fviewdoc%2Fdownload%3Fdoi%3D10.1.1.79.7003%26rep%3Drep1%26type%3Dpdf&ei=RQozTJb6L4_fcbCV6cMD&usg=AFQjCNGClLIB9MAvbrEj45SyUx9cYubLyA&sig2=hg5Wnfy3dLPkoIc1hqSfjg+--+-- Aapo Hyvärinen and Erkki Oja+-- Independent Component Analysis: Algorithms and Applications+-- Neural Networks, 13(4-5):411-430, 2000+--+-----------------------------------------------------------------------------++module Numeric.Statistics.ICA (+ sigmoid, sigmoid',+ demean, whiten,+ ica, icaDefaults+ ) where+++-----------------------------------------------------------------------------++import qualified Data.Array.IArray as I ++import Data.Packed.Vector+import Data.Packed.Matrix+--import Data.Packed.Random++import Numeric.LinearAlgebra.Interface+import Numeric.LinearAlgebra.Algorithms++import Numeric.GSL.Statistics++import Numeric.Statistics++import Control.Monad(replicateM)++import System.Random++-----------------------------------------------------------------------------++-- | sigmoid transfer function+sigmoid :: Double -> Double+sigmoid u = u * exp((-u**2)/2)++-- | derivative of sigmoid transfer function+sigmoid' :: Double -> Double+sigmoid' u = -u**2 * exp((-u**2)/2)++-----------------------------------------------------------------------------++-- preprocessing:+-- demean+-- whiten+-- eigenvalue decomposition of covariance matrix E{xx^T} = EDE^T+-- E orthogonal matrix of eigenvectors+-- D diagonal matrix of eigenvalues, D = diag(d_1,...,d_n)+-- x_white = ED^{-1/2}E^Tx+-- D^{-1/2} = diag{d_1^{-1/2},...}+--++-----------------------------------------------------------------------------++-- | remove the mean from data+demean :: I.Array Int (Vector Double) -- ^ the data+ -> (I.Array Int (Vector Double),Vector Double) -- ^ (demeaned data,mean)+demean d = let u = I.elems $ fmap mean d+ d' = I.listArray (I.bounds d) (zipWith (-) (I.elems d) (map scalar u))+ u' = fromList u+ in (d',u')++-- | whiten data+whiten :: I.Array Int (Vector Double) -- ^ the data+ -> Double -- ^ eigenvalue threshold+ -> (I.Array Int (Vector Double),Matrix Double) -- ^ (whitened data,transform)+whiten d q = let cv = covarianceMatrix d+ (val',vec') = eigSH cv -- the covariance matrix is real symmetric+ val = toList val'+ vec = toColumns vec'+ v' = zip val vec+ v = filter ((> q) . fst) v' -- keep only eigens > than parameter+ (dd',e') = unzip v+ dd = diag $ (** (-0.5)) $ fromList dd' -- square root of eigenvalues diagonalised+ e = fromColumns e'+ x = fromRows $ I.elems d+ t = e <> dd <> trans e -- the actual mathematics+ x' = t <> x -- the actual mathematics+ d' = I.listArray (I.bounds d) (toRows x') + in (d',t)+ +-----------------------------------------------------------------------------++-- assuming that a weight vector is a row++-- algorithm:+-- 1 initial random weight vectors w_i+-- 2 w_i^+ = E{xg(w^Tx)} - E{g'(w^Tx)}w (newton phase)+-- 3 W = W = (WW^T)^{-1/2)W (decorrelation) W = ( ..., w_i, ...)^T+-- WW^T = FDF^T (eigenvalue decomposition)+-- 4 w_i = w^+/norm(w^+) (normalisation) (almost any norm but not Frobenius)+-- 5 if not converged (dot w w^+ ~ 1 implies convergence) go to step 2+--+-- in matrix form, 2 becomes:+-- W^+ = W + (diag a_i)[(diag b_i) + E{g(y)y^T}]W+--+-- where+-- y = Wx+-- b_i = -E{y_ig(y_i)}+-- a_i = -1/(b_i-E{g'(y_i)})+--+-- g(u) = tanh(au) 0<=a<=2, often a = 1+-- g(u) = u exp(-u^2/2)++-----------------------------------------------------------------------------++unconcat 0 _ _ = []+unconcat (r+1) c xs = [take c xs] ++ unconcat r c (drop c xs)++random_vector :: Int -> (Int,Int) -> Matrix Double+random_vector s (r,c) = fromLists $ unconcat r c $ randomRs (-1,1) (mkStdGen s)++-- g g' w x -> w'+update :: (Double -> Double) -> (Double -> Double) -> Matrix Double -> Matrix Double -> Matrix Double+update g g' w x = let y = w <> x+ ys = toRows y+ bis = map (\y' -> - mean (y' * (mapVector sigmoid y'))) ys+ ais = zipWith (\b y' -> -1 / (b - mean (mapVector sigmoid y'))) bis ys+ r = rows y+ ix = ((1,1),(r,r))+ cov = fromArray2D $ I.listArray ix $ map (\(m,n) -> covariance (mapVector sigmoid' (ys!!(m-1))) (ys!!(n-1))) $ I.range ix+ in w + (diag $ fromList ais) <> ((diag $ fromList bis) + cov) <> w ++decorrelate :: Matrix Double -> Matrix Double+decorrelate w = let w' = w / (scalar $ sqrt $ pnorm PNorm2 (w <> trans w))+ in decorrelate w w'+ where decorrelate w w' + | converged 0.000001 w w' = w'+ | otherwise = decorrelate w' ((scale 1.5 w') - (scale 0.5 (w <> trans w <> w)))+{- don't know how to do svd of non-square matrices+decorrelate m = let (u,d,v) = svd m+ in u <> (diag (d ** (-0.5))) <> trans v <> m+-}++normalise :: NormType -> Matrix Double -> Matrix Double+normalise t m = fromRows $ map (\v -> v / (scalar $ pnorm t v)) (toRows m)++converged :: Double -> Matrix Double -> Matrix Double -> Bool+converged t m m' = let d' = map ((-) 1) $ zipWith dot (toRows m) (toRows m')+ in maximum d' <= t++-----------------------------------------------------------------------------++ica' :: (Double -> Double) -- ^ transfer function (tanh,u exp(u^2/2), etc...)+ -> (Double -> Double) -- ^ derivative of transfer function+ -> NormType -- ^ type of normalisation: Infinity, PNorm1, PNorm2+ -> Double -- ^ convergence tolerance for feature vectors+ -> Matrix Double -- ^ weight matrix+ -> [Matrix Double] -- ^ input data in chunks+ -> Matrix Double -- ^ ica transform (weight matrix)+ica' g g' n t w (x:xs) = let w' = normalise n $ decorrelate $ update g g' w x+ in if converged t w w' + then w'+ else ica' g g' n t w' (xs ++ [x])++ica :: Int -- ^ random seed+ -> (Double -> Double) -- ^ transfer function (tanh,u exp(u^2/2), etc...)+ -> (Double -> Double) -- ^ derivative of transfer function+ -> NormType -- ^ type of normalisation: Infinity, PNorm1, PNorm2+ -> Double -- ^ convergence tolerance for feature vectors+ -> Int -- ^ output dimensions+ -> Int -- ^ sampling size (must be smaller than length of data)+ -> I.Array Int (Vector Double) -- ^ data+ -> (I.Array Int (Vector Double),Matrix Double) -- ^ transformed data, ica transform+ica r g g' n t o s a = let i = I.rangeSize $ I.bounds a+ w = random_vector s (o,i)+ x' = fromRows $ I.elems a+ -- next line is BAD if distribution not stationary+ x = concat $ toBlocksEvery i s x'+ w' = ica' g g' n t w x+ y = w' <> x'+ in (I.listArray (1,o) $ toRows y,w') ++-----------------------------------------------------------------------------++-- | ICA with default values: no dimension reduction, euclidean norms, 16 sample groups, sigmoid+icaDefaults :: Int -- ^ random seed+ -> I.Array Int (Vector Double) -- ^ data+ -> (I.Array Int (Vector Double),Matrix Double) -- ^ transformed data, ica transform+icaDefaults r a = let c = I.rangeSize $ I.bounds a+ s = (dim $ (a I.! 1)) `div` 16+ in ica r sigmoid sigmoid' PNorm2 0.0000001 (c-1) s a++-----------------------------------------------------------------------------
lib/Numeric/Statistics/Information.hs view
@@ -45,7 +45,7 @@ -> Vector Double -- ^ the sequence (expected to fall within bounds of Histogram) -> Double -- ^ the entropy entropy p x = let ps = H.prob p x- in dot ps (logE ps)+ in negate $ dot ps (logE ps) -- | the mutual information \sum_x \sum_y \ln{\frac{p(x,y)}{p(x)p(y)}} mutual_information :: H2.Histogram2D -- ^ the underlying distribution@@ -56,6 +56,6 @@ mutual_information p px py z@(x,y) = let ps = H2.prob p z xs = H.prob px x ys = H.prob py y- in dot ps (logE ps - logE (xs*ys)) + in negate $ dot ps (logE ps - logE (xs*ys)) -----------------------------------------------------------------------------
+ lib/Numeric/Statistics/PCA.hs view
@@ -0,0 +1,73 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Statistics.PCA+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : portable+--+-- Principal Components Analysis+--+-----------------------------------------------------------------------------++module Numeric.Statistics.PCA (+ pca, pcaTransform, pcaReduce+ ) where+++-----------------------------------------------------------------------------++import qualified Data.Array.IArray as I ++import Data.Packed.Vector+import Data.Packed.Matrix++import Numeric.LinearAlgebra.Interface+import Numeric.LinearAlgebra.Algorithms++import Numeric.GSL.Statistics++import Numeric.Statistics++-----------------------------------------------------------------------------++-- | find the n principal components of multidimensional data+pca :: I.Array Int (Vector Double) -- the data+ -> Double -- eigenvalue threshold+ -> Matrix Double+pca d q = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension+ cv = covarianceMatrix d'+ (val',vec') = eigSH cv -- the covariance matrix is real symmetric+ val = toList val'+ vec = toColumns vec'+ v' = zip val vec+ v = filter (\(x,_) -> x > q) v' -- keep only eigens > than parameter+ in fromColumns $ snd $ unzip v++-- | perform a PCA transform of the original data (remove mean)+-- | Final = M^T Data^T+pcaTransform :: I.Array Int (Vector Double) -- ^ the data+ -> Matrix Double -- ^ the principal components+ -> I.Array Int (Vector Double) -- ^ the transformed data+pcaTransform d m = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension+ in I.listArray (1,cols m) $ toRows $ (trans m) <> (fromRows $ I.elems d')++-- | perform a dimension-reducing PCA modification+pcaReduce :: I.Array Int (Vector Double) -- ^ the data+ -> Double -- ^ eigenvalue threshold+ -> I.Array Int (Vector Double) -- ^ the reduced data, with n principal components+pcaReduce d q = let u = fmap (scalar . mean) d+ d' = zipWith (-) (I.elems d) (I.elems u)+ cv = covarianceMatrix $ I.listArray (I.bounds d) d'+ (val',vec') = eigSH cv -- the covariance matrix is real symmetric+ val = toList val'+ vec = toColumns vec'+ v' = zip val vec+ v = filter (\(x,_) -> x > q) v' -- keep only eigens > than parameter+ m = fromColumns $ snd $ unzip v+ in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (trans m) <> fromRows d') (I.elems u) ++-----------------------------------------------------------------------------