hstatistics 0.2.0.8 → 0.2.1.1
raw patch · 8 files changed
+158/−32 lines, 8 files
Files
- CHANGES +4/−1
- hstatistics.cabal +3/−1
- lib/Numeric/Statistics/Histogram.hs +1/−1
- lib/Numeric/Statistics/ICA.hs +14/−13
- lib/Numeric/Statistics/Information.hs +14/−15
- lib/Numeric/Statistics/PCA.hs +1/−1
- lib/Numeric/Statistics/PDF.hs +61/−0
- lib/Numeric/Statistics/Surrogate.hs +60/−0
CHANGES view
@@ -45,4 +45,7 @@ error in decorrelation fixed -0.2.0.9:+0.2.1.1:+ added PDF.hs+ modified Information to take in PDF addition+ added Surrogate.hs
hstatistics.cabal view
@@ -1,5 +1,5 @@ Name: hstatistics-Version: 0.2.0.8+Version: 0.2.1.1 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010@@ -31,8 +31,10 @@ Exposed-modules: Numeric.Statistics Numeric.Statistics.PCA Numeric.Statistics.ICA+ Numeric.Statistics.PDF Numeric.Statistics.Information Numeric.Statistics.Histogram+ Numeric.Statistics.Surrogate other-modules: C-sources:
lib/Numeric/Statistics/Histogram.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.Statistics.Histogram--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com
lib/Numeric/Statistics/ICA.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.Statistics.ICA--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -13,11 +13,11 @@ -- -- 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 ----- Aapo Hyvärinen and Erkki Oja--- Independent Component Analysis: Algorithms and Applications--- Neural Networks, 13(4-5):411-430, 2000+-- <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> -- ----------------------------------------------------------------------------- @@ -138,12 +138,12 @@ 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 :: Matrix Double -> Matrix Double-decorrelate w = let w' = w / (scalar $ sqrt $ pnorm PNorm2 (w <> trans w))- in decorrelate' w w'- where decorrelate' m m' - | converged 0.000001 m m' = m'- | otherwise = decorrelate' m' ((scale 1.5 m') - (scale 0.5 (m' <> trans m' <> m')))+decorrelate :: NormType -> Double -> Matrix Double -> Matrix Double+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@@ -166,11 +166,12 @@ -> [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 $ update g g' w x+ica' g g' n t w (x:xs) = let w' = normalise n $ decorrelate n t $ update g g' w x in if converged t w w' then w' else ica' g g' n t w' (xs ++ [x]) +-- | perform an ICA transform ica :: Int -- ^ random seed -> (Double -> Double) -- ^ transfer function (tanh,u exp(u^2/2), etc...) -> (Double -> Double) -- ^ derivative of transfer function@@ -197,6 +198,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' PNorm2 0.0000001 (c-1) s a+ in ica r sigmoid sigmoid' PNorm1 0.0000001 (c-1) s a -----------------------------------------------------------------------------
lib/Numeric/Statistics/Information.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.Statistics.Information--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com@@ -23,8 +23,7 @@ import Data.Packed.Vector -import qualified Numeric.GSL.Histogram as H-import qualified Numeric.GSL.Histogram2D as H2+import Numeric.Statistics.PDF import Numeric.LinearAlgebra.Algorithms import Numeric.LinearAlgebra.Interface()@@ -41,21 +40,21 @@ ----------------------------------------------------------------------------- -- | the entropy \sum p_i l\ln{p_i} of a sequence-entropy :: H.Histogram -- ^ the underlying distribution- -> Vector Double -- ^ the sequence (expected to fall within bounds of Histogram)+entropy :: PDF a Double => a -- ^ the underlying distribution+ -> Vector Double -- ^ the sequence -> Double -- ^ the entropy-entropy p x = let ps = H.prob p x+entropy p x = let ps = probability p x 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- -> H.Histogram -- ^ the first dimension distribution- -> H.Histogram -- ^ the second dimension distribution- -> (Vector Double, Vector Double) -- ^ the sequence (expected to fall within bounds of Histogram)+-- | 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+ -> a -- ^ the first dimension distribution+ -> a -- ^ the second dimension distribution+ -> (Vector Double, Vector Double) -- ^ the sequence -> Double -- ^ the mutual information-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 negate $ dot ps (logE ps - logE (xs*ys)) +mutual_information p px py (x,y) = let ps = probability p $ zipVector (,) x y+ xs = probability px x+ ys = probability py y+ in negate $ dot ps (logE ps - logE (xs*ys)) -----------------------------------------------------------------------------
lib/Numeric/Statistics/PCA.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.Statistics.PCA--- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com
+ lib/Numeric/Statistics/PDF.hs view
@@ -0,0 +1,61 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Statistics.PDF+-- Copyright : (c) A. V. H. McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : portable+--+-- Probability Density Function interface+--+-----------------------------------------------------------------------------++module Numeric.Statistics.PDF (+ -- * PDF+ PDFFunction,+ PDF(..),+ -- * Creation+ pdfFromFunction+ ) where+++-----------------------------------------------------------------------------++import qualified Data.Packed.Vector as V++import qualified Numeric.GSL.Histogram as H+import qualified Numeric.GSL.Histogram2D as H2++import Foreign.Storable+--import Numeric.Statistics++-----------------------------------------------------------------------------++-- ^ a probability density function+data PDFFunction a = P_Func (a -> Double) -- p(x)++-----------------------------------------------------------------------------++-- ^ a PDF interface+class PDF b a where+ -- ^ calculate a probability+ probability :: b -> V.Vector a -> V.Vector Double++instance Storable b => PDF (PDFFunction b) b where+ probability (P_Func f) = V.mapVector f++instance PDF H.Histogram Double where+ probability = H.prob++instance PDF H2.Histogram2D (Double,Double) where+ probability = H2.probPaired++-----------------------------------------------------------------------------++-- | create a PDF from an arbtrary function f :-> [0,1]+pdfFromFunction :: (a -> Double) -> PDFFunction a+pdfFromFunction = P_Func+
+ lib/Numeric/Statistics/Surrogate.hs view
@@ -0,0 +1,60 @@+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Statistics.Surrogate+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : portable+--+-- Methods for tests using surrogate data+--+-----------------------------------------------------------------------------++module Numeric.Statistics.Surrogate (+ surrogate+ ) where+++-----------------------------------------------------------------------------++import Data.Packed.Vector+--import Data.Packed.Matrix++import qualified Data.Array.IArray as I ++import Numeric.GSL.Permutation++import System.Random++-----------------------------------------------------------------------------++-- | perform an analysis using surrogate data+surrogate :: Int -- ^ random seed+ -> Int -- ^ number of repetitions+ -> (I.Array Int (Vector Double) -> a) -- ^ the evaluation function+ -> I.Array Int (Vector Double) -- ^ the data+ -> I.Array Int a -- ^ the results, with the evaluated real data in position 1 and the rest of the array containing the evaluated surrogate data+ +surrogate r n f d = I.listArray (1,n+1) $ (f d) : (surrogate' (mkStdGen r) n f d)++surrogate' :: StdGen -> Int-> (I.Array Int (Vector Double) -> a) -> I.Array Int (Vector Double) -> [a]+surrogate' _ 0 _ _ = []+surrogate' g (n+1) f d = let (g',g'') = split g+ d' = permute_data g' d+ in (f d) : (surrogate' g'' n f d)++randomList :: StdGen -> Int -> [Int]+randomList _ 0 = []+randomList g (n+1) = let (r,g') = random g+ in r : (randomList g' n)++permute_data :: StdGen -> I.Array Int (Vector Double) -> I.Array Int (Vector Double)+permute_data g d = let s = I.rangeSize $ I.bounds d+ rs = randomList g s+ ds = zip rs $ I.elems d+ in I.listArray (I.bounds d) $ map (\(r,v) -> permute (random_permute r (dim v)) v) ds+ +-----------------------------------------------------------------------------