diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -49,3 +49,6 @@ added PDF.hs modified Information to take in PDF addition added Surrogate.hs + +0.2.2.1: + improved ICA decorrelation step diff --git a/hstatistics.cabal b/hstatistics.cabal --- a/hstatistics.cabal +++ b/hstatistics.cabal @@ -1,5 +1,5 @@ Name: hstatistics -Version: 0.2.1.1 +Version: 0.2.2.1 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010 @@ -8,7 +8,11 @@ Stability: provisional Homepage: http://code.haskell.org/hstatistics Synopsis: Statistics -Description: Purely functional interface for statistics based on hmatrix and hmatrix-gsl-stats +Description: + Purely functional interface for statistics based on hmatrix and hmatrix-gsl-stats + . + When hmatrix is installed with -fvector, the vector type is Data.Vector.Storable + from the vector package. Category: Math, Statistics tested-with: GHC ==6.12.1 diff --git a/lib/Numeric/Statistics.hs b/lib/Numeric/Statistics.hs --- a/lib/Numeric/Statistics.hs +++ b/lib/Numeric/Statistics.hs @@ -14,6 +14,7 @@ ----------------------------------------------------------------------------- module Numeric.Statistics ( + Sample,Samples, covarianceMatrix ) where @@ -29,8 +30,13 @@ ----------------------------------------------------------------------------- +type Sample a = Vector a +type Samples a = I.Array Int (Vector a) + +----------------------------------------------------------------------------- + -- | the covariance matrix -covarianceMatrix :: I.Array Int (Vector Double) -- ^ the dimensions of data (each vector being one dimension) +covarianceMatrix :: Samples 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 --- a/lib/Numeric/Statistics/ICA.hs +++ b/lib/Numeric/Statistics/ICA.hs @@ -138,15 +138,16 @@ cov = fromArray2D $ I.listArray ix $ map (\(m,n) -> covariance (mapVector g' (ys!!(m-1))) (ys!!(n-1))) $ I.range ix in w + (diag $ fromList ais) <> ((diag $ fromList bis) + cov) <> w -decorrelate :: NormType -> Double -> Matrix Double -> Matrix Double -decorrelate n t w = let w' = w / (scalar $ sqrt $ pnorm n (w <> trans w)) +decorrelate :: Matrix Double -> Matrix Double +decorrelate m = let (d',v') = eig m + d = fst $ fromComplex d' + v = fst $ fromComplex v' + in v <> (diag (d ** (-0.5))) <> trans v <> m +{-decorrelate n t w = let w' = w / (scalar $ sqrt $ pnorm n (w <> trans w)) in decorrelate' t w w' where decorrelate' t' m m' | converged t' m m' = m' | otherwise = decorrelate' t' m' ((scale 1.5 m') - (scale 0.5 (m' <> trans m' <> m'))) -{- 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 @@ -166,7 +167,7 @@ -> [Matrix Double] -- ^ input data in chunks -> Matrix Double -- ^ ica transform (weight matrix) ica' _ _ _ _ _ [] = error "no sample data" -ica' g g' n t w (x:xs) = let w' = normalise n $ decorrelate n t $ update g g' w x +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]) @@ -177,18 +178,18 @@ -> (Double -> Double) -- ^ derivative of transfer function -> NormType -- ^ type of normalisation: Infinity, PNorm1, PNorm2 -> Double -- ^ convergence tolerance for feature vectors - -> Int -- ^ output dimensions +-- -> 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 r (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 r g g' n t s a = let i = I.rangeSize $ I.bounds a + w = random_vector r (i,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,1) $ toRows y,w') ----------------------------------------------------------------------------- @@ -198,6 +199,6 @@ -> (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' PNorm1 0.0000001 (c-1) s a + in ica r sigmoid sigmoid' PNorm1 0.0000001 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 @@ -40,14 +40,16 @@ ----------------------------------------------------------------------------- -- | the entropy \sum p_i l\ln{p_i} of a sequence -entropy :: PDF a Double => a -- ^ the underlying distribution +entropy :: PDF a Double + => a -- ^ the underlying distribution -> Vector Double -- ^ the sequence -> Double -- ^ the entropy entropy p x = let ps = probability p x in negate $ dot ps (logE ps) -- | the mutual information \sum_x \sum_y p(x,y) \ln{\frac{p(x,y)}{p(x)p(y)}} -mutual_information :: (PDF a Double, PDF b (Double,Double)) => b -- ^ the underlying distribution +mutual_information :: (PDF a Double, PDF b (Double,Double)) + => b -- ^ the underlying distribution -> a -- ^ the first dimension distribution -> a -- ^ the second dimension distribution -> (Vector Double, Vector Double) -- ^ the sequence