diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -31,3 +31,7 @@ 0.2.0.5: added Histogram + +0.2.0.6: + added PCA + added ICA diff --git a/INSTALL b/INSTALL --- a/INSTALL +++ b/INSTALL @@ -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 diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,2 +1,2 @@ -Copyright Alberto Ruiz 2006-2007 +Copyright A.V.H. McPhail 2010 GPL license diff --git a/hstatistics.cabal b/hstatistics.cabal --- a/hstatistics.cabal +++ b/hstatistics.cabal @@ -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: diff --git a/lib/Numeric/Statistics.hs b/lib/Numeric/Statistics.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Statistics.hs @@ -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 gmail 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)) + +----------------------------------------------------------------------------- diff --git a/lib/Numeric/Statistics/ICA.hs b/lib/Numeric/Statistics/ICA.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Statistics/ICA.hs @@ -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 gmail 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 + +----------------------------------------------------------------------------- diff --git a/lib/Numeric/Statistics/Information.hs b/lib/Numeric/Statistics/Information.hs --- a/lib/Numeric/Statistics/Information.hs +++ b/lib/Numeric/Statistics/Information.hs @@ -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)) ----------------------------------------------------------------------------- diff --git a/lib/Numeric/Statistics/PCA.hs b/lib/Numeric/Statistics/PCA.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Statistics/PCA.hs @@ -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 gmail 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) + +-----------------------------------------------------------------------------