diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Lennart Schmitt 2011 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Lennart Schmitt nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/linda.cabal b/linda.cabal new file mode 100644 --- /dev/null +++ b/linda.cabal @@ -0,0 +1,62 @@ +-- lda.cabal auto-generated by cabal init. For additional options, see +-- http://www.haskell.org/cabal/release/cabal-latest/doc/users-guide/authors.html#pkg-descr. +-- The name of the package. +Name: linda + +-- The package version. See the Haskell package versioning policy +-- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for +-- standards guiding when and how versions should be incremented. +Version: 0.1 + +-- A short (one-line) description of the package. +Synopsis: LINear Discriminant Analysis + +-- A longer description of the package. +Description: This package (mainly the module LDA) implements the linear discriminant analysis. It provides both data classification (according to Fisher) and data analysis (by discriminant criteria). + +-- The license under which the package is released. +License: BSD3 + +-- The file containing the license text. +License-file: LICENSE + +-- The package author(s). +Author: Lennart Schmitt + +-- An email address to which users can send suggestions, bug reports, +-- and patches. +Maintainer: lennart...schmitt@gmail.com + +-- A copyright notice. + Copyright: (c) Lennart Schmitt + +-- Stability of the pakcage (experimental, provisional, stable...) +Stability: Experimental + +Category: Math, Statistics + +Build-type: Simple + +-- Extra files to be distributed with the package, such as examples or +-- a README. +-- Extra-source-files: + +-- Constraint on the version of Cabal needed to build this package. +Cabal-version: >=1.2 + + +Library + -- Modules exported by the library. + Exposed-modules: Numeric.Function, Numeric.Statistics.LDA, Numeric.MatrixList, Numeric.Matrix, Numeric.Vector + + -- Packages needed in order to build this package. + Build-depends: base >= 2 && <= 4, hstats -any, hmatrix -any, HUnit >= 1 + + -- Modules not exported by this package. + Other-modules: Tests + + hs-source-dirs: src + + -- Extra tools (e.g. alex, hsc2hs, ...) needed to build the source. + -- Build-tools: + diff --git a/src/Numeric/Function.hs b/src/Numeric/Function.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Function.hs @@ -0,0 +1,47 @@ +-- While working on this module you are encouraged to remove it and fix +-- any warnings in the module. See +-- http://hackage.haskell.org/trac/ghc/wiki/WorkingConventions#Warnings +-- for details + +----------------------------------------------------------------------------- +-- | +-- Module : Function +-- Copyright : (c) Lennart Schmitt +-- License : BSD-style (see the file libraries/base/LICENSE) +-- +-- Maintainer : libraries@haskell.org +-- Stability : experimental +-- Portability : portable +-- +-- This module implements the simple functionality of multidimensional linear function calculation. +----------------------------------------------------------------------------- + +module Numeric.Function ( + -- * Data Types + LinFunction, + Values, + -- * Functions + calcLinFunction + )where + +-- | The function-type represents a function by its constants, e.g. +-- +-- > [b0,b1,...,bn] +-- +-- represents the function f = b0 + b1 * X1 + ... + bn * Xn +type LinFunction a = [a] + +-- | Similare to the function-type, but die value-type represents the values of the variables in a function, e.g. +-- +-- > [X1,...,Xn] +type Values a = [a] + +-- | Calculates the result of a given function with given values, e.g. +-- +-- > calcLinFunction [1,1,1] [1,2] == 1 + 1 * 1 + 1 * 2 == 4 +-- +-- > calcLinFunction [1,2,3] [1,1] == 1 + 2 * 1 + 3 * 1 == 6 +-- +-- > calcLinFunction [1,2,3] [1..] == 1 + 2 * 1 + 3 * 2 == 9 +calcLinFunction :: Num a => LinFunction a -> Values a -> a +calcLinFunction f = sum . (zipWith (*) f ) . ([1]++) diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Matrix.hs @@ -0,0 +1,117 @@ +-- While working on this module you are encouraged to remove it and fix +-- any warnings in the module. See +-- http://hackage.haskell.org/trac/ghc/wiki/WorkingConventions#Warnings +-- for details + +----------------------------------------------------------------------------- +-- | +-- Module : Matrix +-- Copyright : (c) Lennart Schmitt +-- License : BSD-style (see the file libraries/base/LICENSE) +-- +-- Maintainer : libraries@haskell.org +-- Stability : experimental +-- Portability : portable +-- +-- This module includes a few (standard) functions to work with matrixes. +----------------------------------------------------------------------------- +module Numeric.Matrix ( + -- * Data Types + Matrix, -- | Represents a matrix + RawMatrix, + Element, -- | The elements of a matrix + + -- * Convert a Matrix into another Type + flatten, -- | Flattens a matrix to a vector + toColumns, -- | Converts the representation of a matrix to a list of vectors by its columns. + toLists, -- | Converts the representation of a matrix to a list of lists. (NOT equivalent with toRows!!!) + toRows, -- | Converts the representation of a matrix to a list of vectors by its rows. + + -- * Convert some Data into Matrix-Type + asRow, -- | Converts the representation of a matrix from a list of vectors to a matrix (by rows) + fromLists, -- | Converts the representation of a matrix from a list of lists to a matrix (NOT equivalent with asRows!!!) + fromListToQuadraticMatrix, + + -- * Matrix calculations/transformations + (@@>), -- | Get the matrix-element in row x and col y (like "(!!)" for lists) + eigenvalue, + eigenvector, + inv, -- | Calculates the invariant matrix of the input matrix + identityMatrix, + mapMatrix, + reduceMatrix, + scalarMultiplication, + subtractMatrix, + trans, -- | Transposes a matrix + zipAllWith + ) where + +-- This module "extends" an existing matrix-module by using its datatypes and +-- it also uses linear algebra algorithms (so it is based on that modules too) +import Data.Packed.Matrix (Matrix, buildMatrix, fromLists, toLists, flatten, reshape, trans, rows,asRow,toRows, toColumns, (@@>), Element) +import Numeric.LinearAlgebra.Algorithms (eig, eigenvalues,inv) +-- There are also a few additional modules that are needed to work with the matrixes +import Numeric.Vector (Vector,RawVector,toList,fromList,count,maximumBy,transpose) +import Foreign.Marshal.Utils (fromBool) +import Data.Complex (realPart) + +-- | A matrix represented by a list of lists. +type RawMatrix a = [[a]] + +-- | Calculates the identity matrix (n x n) by given scale (n) +identityMatrix :: Int -> Matrix Double +identityMatrix i = buildMatrix i i (\ (x,y) -> fromBool (x == y)) + +-- | A simple map-Function which maps a given function on every element of the given matrix +mapMatrix :: (Double -> Double) -> Matrix Double -> Matrix Double +mapMatrix f m = fromLists.(map$map f).toLists $ m + +-- | Calculates the scalarproduct (with a scalar and matrix) +scalarMultiplication :: Double -> Matrix Double -> Matrix Double +scalarMultiplication x m = mapMatrix ((*) x) m + +-- | Calculates the difference (matrix) between two matrixes +subtractMatrix :: Matrix Double -> Matrix Double -> Matrix Double +subtractMatrix a b = fromListToQuadraticMatrix $ zipWith (-) (toList . flatten $a) (toList . flatten $b) + +-- | Builds a quadratic matrix out of a list +fromListToQuadraticMatrix :: [Double] -> Matrix Double +fromListToQuadraticMatrix xs = (reshape (round$ sqrt(count xs)) (fromList xs)) + +-- | Calculates the eigenvalue of a matrix, e.g. +-- +-- > eigenvalue (fromLists [[0.77143,-0.25714],[-0.42245,0.14082]]) +-- +-- returns +-- +-- > 0.9122456375784809 +eigenvalue :: Matrix Double -> Double +eigenvalue = ((maximumBy (compare . abs)) . (map realPart) . toList . eigenvalues) + +-- | Calculates one eigenvector of a given matrix, e.g. +-- +-- > eigenvector (fromLists [[-0.14081563757848092,-0.25714],[-0.42245,-0.7714256375784809]]) +-- +-- returns +-- +-- > [0.8770950095147589,-0.48031692067249215] +eigenvector :: Matrix Double -> Vector Double +eigenvector = fromList . (map realPart) . head . toLists . trans . snd . eig + +-- | Calculates the reduced matrix of a given matrix (by reducing the given matrix), e.g. +-- +-- > reduceMatrix (fromLists [[0.77143,-0.25714],[-0.42245,0.14082]]) +-- +-- returns +-- +-- > (2><2)[ -0.14081563757848092,-0.25714,-0.42245,-0.7714256375784809] +reduceMatrix :: Matrix Double -> Matrix Double +reduceMatrix a = subtractMatrix a (scalarMultiplication (eigenvalue a) (identityMatrix $rows a)) + +-- | Zipps a matrix col by col +-- +-- > zipAllWith sum [[1,2,3],[1,2,3],[1,2,3]] == [3,6,9] +zipAllWith :: (RawVector a -> b) -> RawMatrix a -> RawVector b +zipAllWith _ [] = [] +zipAllWith f xss = map f . transpose $ xss + diff --git a/src/Numeric/MatrixList.hs b/src/Numeric/MatrixList.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MatrixList.hs @@ -0,0 +1,108 @@ +-- While working on this module you are encouraged to remove it and fix +-- any warnings in the module. See +-- http://hackage.haskell.org/trac/ghc/wiki/WorkingConventions#Warnings +-- for details + +----------------------------------------------------------------------------- +-- | +-- Module : MatrixList +-- Copyright : (c) Lennart Schmitt +-- License : BSD-style (see the file libraries/base/LICENSE) +-- +-- Maintainer : libraries@haskell.org +-- Stability : experimental +-- Portability : portable +-- +-- This module implements a list of matrixes and a few functions to handle them. +----------------------------------------------------------------------------- + +module Numeric.MatrixList ( + -- * Data Types + MatrixList, + RawMatrixList, + -- * Functions + Numeric.MatrixList.toLists, + countMatrixElements, + countMatrixes, + countMatrixesCols, + averages, + crossProduct, + mapVectors, + foldVectors, + mapElements, + transposeAll, + countRows, + countAllRows + )where + +import Numeric.Vector +import Numeric.Matrix + +-- | A list of matrixes +type MatrixList a = [Matrix a] + +-- | A list of matrixes represented as a list of lists of lists +type RawMatrixList a = [[[a]]] + +-- | Transforms a Matrixlist to a RawMatrixList. +toLists :: (Numeric.Matrix.Element a) => MatrixList a -> RawMatrixList a +toLists = (map (map toList) ) . (map toColumns) + +-- | Counts the number of elements per matrix +-- @ unused +countMatrixElements :: RawMatrixList a -> RawMatrix Double +countMatrixElements = foldVectors count + +-- | Counts the number of matrixes in a list of matrixes +countMatrixes :: RawMatrixList a -> Double +countMatrixes = count + +-- | Counts the number of cols (based on the guess that all matrixes have the similare structure) +countMatrixesCols :: RawMatrixList a -> Double +countMatrixesCols = count . head . head + +-- | Calculate every cols averages +-- +-- > averages [[[1,2],[2,1]],[[2,3],[3,4]]] == [[1.5,1.5],[2.5,3.5]] +averages :: RawMatrixList Double -> RawMatrix Double +averages = foldVectors average + +-- | Calculates the cross-product of every matrix-cols +-- +-- > crossProduct [[[1,2],[2,1]],[[2,3],[3,4]]] == [[2.0,2.0],[6.0,12.0]] +crossProduct :: RawMatrixList Double -> RawMatrix Double +crossProduct = foldVectors product + +-- | maps a function over every vector of a list of matrixes +-- +-- > mapVectors (map (1+)) [[[1,2],[2,1]],[[2,3],[3,4]]] == [[[2.0,3.0],[3.0,2.0]],[[3.0,4.0],[4.0,5.0]]] +mapVectors :: ([a] -> [b]) -> RawMatrixList a -> RawMatrixList b +mapVectors f = map $ map $ f + +-- | folds every vector of a list of matrixes +-- +-- > foldVectors sum [[[1,2],[2,1]],[[2,3],[3,4]]] == [[[2.0,3.0],[3.0,2.0]],[[3.0,4.0],[4.0,5.0]]] +foldVectors :: ([a] -> b) -> RawMatrixList a -> RawMatrix b +foldVectors f = map $ map $ f + +-- | maps a function over every element of a list of matrixes +-- +-- > mapElements (1+) [[[1,2],[2,1]],[[2,3],[3,4]]] == [[[2.0,3.0],[3.0,2.0]],[[3.0,4.0],[4.0,5.0]]] +mapElements :: (a -> b) -> RawMatrixList a -> RawMatrixList b +mapElements f = mapVectors $ map f + +-- | Transposes every matrix in a lit of matrixes +transposeAll :: RawMatrixList a -> RawMatrixList a +transposeAll = map transpose + +-- | Counts the rows of every matrix in the list +-- +-- > countRows [[[1,2],[2,1]],[[2,3],[3,4],[1,1]]] == [2.0,3.0] +countRows :: RawMatrixList a -> RawVector Double +countRows = map count + +-- | Counts the sum of all matrixes-rows +-- +-- > countAllRows [[[1,2],[2,1]],[[2,3],[3,4],[1,1]]] == 5.0 +countAllRows :: RawMatrixList a -> Double +countAllRows = sum . countRows diff --git a/src/Numeric/Statistics/LDA.hs b/src/Numeric/Statistics/LDA.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Statistics/LDA.hs @@ -0,0 +1,245 @@ +-- While working on this module you are encouraged to remove it and fix +-- any warnings in the module. See +-- http://hackage.haskell.org/trac/ghc/wiki/WorkingConventions#Warnings +-- for details + +----------------------------------------------------------------------------- +-- | +-- Module : LDA +-- Copyright : (c) Lennart Schmitt +-- License : BSD-style (see the file libraries/base/LICENSE) +-- +-- Maintainer : libraries@haskell.org +-- Stability : experimental +-- Portability : portable +-- +-- This module implements some linear discriminant analysis functions. +-- Imagine you've made a poll and now you have values/attributes from every subscriber. +-- Further more you've grouped the subscribers into clusters. +-- The poll-datas are structured as follows: +-- +-- * poll-data of one subscriber = [value] --> Vector value +-- +-- * poll-data of one cluster/group of subscribers = [[values]] --> Matrix values +-- +-- * poll-data of all clusters/groups = [[[values]]] --> MatrixList values +-- +-- Now you want to check if you clustered right and/or how significant the values you asked for are... +-- +----------------------------------------------------------------------------- + +module Numeric.Statistics.LDA ( + fisher, + fisher', + fisherT, + fisherAll, + fisherClassificationFunction, + aprioriProbability, + discriminantCriteria, + isolatedDiscriminant +) where + +import Numeric.Matrix +import Numeric.MatrixList +import Numeric.Vector +import Numeric.Function +import Numeric.LinearAlgebra.LAPACK + + +-- | Calculates the difference between every element and the average of the matrixes row. +-- +-- > diffAverage [[[-1,1],[2,2]],[[1,3],[4,8]]] == [[[-1.0,1.0],[0.0,0.0]],[[-1.0,1.0],[-2.0,2.0]]] +diffAverage :: RawMatrixList Double -> RawMatrixList Double +diffAverage xs = ( zipWith $ zipWith (\ _xs y -> map (\x -> (-)x y) _xs ) ) xs (averages xs) + +-- | Calculates the square of the difference ("diffAverage"). +-- +-- > squareDiffAverage [[[-1,1],[2,2]],[[1,3],[4,8]]] == [[[1.0,1.0],[0.0,0.0]],[[1.0,1.0],[4.0,4.0]]] +squareDiffAverage :: RawMatrixList Double -> RawMatrixList Double +squareDiffAverage = ( mapElements (^2) ) . diffAverage + +-- | Calculates the average for every matrix/group. +-- +-- > totalAverages [[[-1,1],[2,2]],[[1,3],[4,8]]] == [1.5,3.5] +totalAverages :: RawMatrixList Double -> RawVector Double +totalAverages xs = map (flip(/) ( sum (countRows xs) )) ( map sum $ zipWith ( \x ys -> map ((*)x) ys) (countRows xs) (transpose $ averages $ transposeAll xs ) ) + +-- | Calculates the sum of the differences from the averages. +-- +-- > sumOfAverages [[[-1,1],[2,2]],[[1,3],[4,8]]] == [[2.0,0.0,0.0,0.0],[2.0,4.0,4.0,8.0]] +sumOfAverages :: RawMatrixList Double -> RawMatrix Double +sumOfAverages xss = ( foldVectors sum (map (\xs -> [zipWith (*) a b | x <- [0..((round (count xs))-1)] , y <- [0..((round (count xs))-1)], let a= xs !! x, let b = xs !! y]) ( diffAverage xss ))) + +-- | Calculates the spread within the cluster/group/matrix. +-- +-- > spreadWithinGroups [[[-1,1],[2,2]],[[1,3],[4,8]]] == (2><2) [ 9.0, 9.0, 9.0, 13.0 ] +spreadWithinGroups :: RawMatrixList Double -> Matrix Double +spreadWithinGroups = fromListToQuadraticMatrix . (zipAllWith sum) . sumOfAverages . transposeAll + +-- | Calculates the spread from total average. +-- +-- > spreadFromTotalAverages [[[-1,1],[2,2]],[[1,3],[4,8]]] == [[-1.0,-2.0],[1.0,2.0]] +spreadFromTotalAverages :: RawMatrixList Double -> RawMatrix Double +spreadFromTotalAverages x = map (zipWith (flip(-)) xs) xss + where + xs = (totalAverages x) + xss = (averages.transposeAll $ x) + +-- | Calculates the spread between the groups/clusters. +-- +-- > spreadBetweenGroups [[[-1,1],[2,2]],[[1,3],[4,8]]] == (2><2) [ 4.0, 8.0, 8.0, 16.0 ] +spreadBetweenGroups :: RawMatrixList Double -> Matrix Double +spreadBetweenGroups = fromListToQuadraticMatrix . (zipAllWith sum). d + where + d xss = [x|i <- [0..((count (spreadFromTotalAverages xss))-1)], let x = c ((spreadFromTotalAverages xss) !!i) ((b xss)!!i)] + c xs ys = ([a*b | x<- [0..((count xs)-1)],y<- [0..((count ys)-1)], let a= xs !! x, let b= ys !!y]) + b xs = map (zipWith (*) (countRows xs) ) $ spreadFromTotalAverages xs + + +-- | Calculates the isolated discriminants of every attribute. +-- +-- > isolatedDiscriminant [[[-1,1],[2,2]],[[1,3],[4,8]]] == [0.4444444444444444,1.2307692307692308] +isolatedDiscriminant :: RawMatrixList Double -> RawVector Double +isolatedDiscriminant xss = [b/w | i<- [0..c], let b = bM @@> (i,i), let w = wM @@> (i,i) ] + where + bM = spreadBetweenGroups xss + wM = spreadWithinGroups xss + c = count xss -1 + +{- ----------------------------------------------- +----------------------Utilities ------------------ +------------------------------------------------ -} + +-- | Calculates a analysis-matrix. +calcA :: RawMatrixList Double -> Matrix Double +calcA xss = cA (spreadWithinGroups xss) (spreadBetweenGroups xss) + where + cA w b = multiplyR (inv w) b + +-- | Calculation of the scaling factor of a matrix. +scalingFactor :: RawMatrixList Double -> Double +scalingFactor xsss = 1 / s + where + s = sqrt (v'Wv / (fallzahl-gruppen)) + v'Wv = head $ head $ (Numeric.Matrix.toLists) (v `multiplyR` w `multiplyR` (trans v)) + fallzahl = countAllRows xsss + gruppen = countMatrixes xsss + v = asRow.head . toRows . trans . calcA $ xsss + w = spreadWithinGroups xsss + +-- | Calculates the scaled discriminant coefficient. +scaledDiscriminantCoefficient :: Matrix Double -> Double -> Matrix Double +scaledDiscriminantCoefficient v s = scalarMultiplication s v + +-- | Calculates the constant element of the discriminant function. +constElement :: Matrix Double -> RawVector Double -> Double +constElement b x = -1 * (sum $ zipWith (*) (toList . flatten $ b) x) + +-- | Calculation of possible discriminant functions. +scaledDiscriminantFunctions :: RawMatrixList Double -> RawVector (LinFunction Double) +scaledDiscriminantFunctions xss = map (\x -> (cE x):x) normDiskKoefs + where + normDiskKoefs = (Numeric.Matrix.toLists).trans $ scaledDiscriminantCoefficient (calcA xss) (scalingFactor xss) + totAvg = totalAverages xss + cE xs = (flip) constElement totAvg (fromLists [xs]) + +-- | Calculates the scaled discriminant for an attribute by using a scaled discriminant function. +scaledDiscriminant :: LinFunction Double -> Values Double -> Double +scaledDiscriminant v x = calcLinFunction v x + +-- | Calculates all scaled discriminants. +scaledDiscriminants :: RawMatrixList Double -> RawMatrixList Double +scaledDiscriminants xss = map (\linFun -> foldVectors (scaledDiscriminant linFun) xss) normDiskFs + where + normDiskFs = scaledDiscriminantFunctions xss + +-- | Calculation of the centroid of a group/cluster by using the discriminant function an the groups values. +centroid :: LinFunction Double -> RawVector(Values Double) -> Double +centroid v m = (sum $ map (scaledDiscriminant v ) m) / (count m) + +-- | Calculate all centroids +centroids :: RawMatrixList Double -> RawMatrix Double +centroids xss = map (\linFun -> map (centroid linFun) xss) normDiskFs + where + normDiskFs = scaledDiscriminantFunctions xss + +-- | Calculates the averages of the scaled discriminants per discriminant function. +averageOfScaledDiscriminants :: RawMatrixList Double -> RawVector Double +averageOfScaledDiscriminants xss = map (\x -> (foldl1 (+) x) / (count x) ) $ map concat $ scaledDiscriminants xss + +-- | Calculates the total spread between all groups. +totalSpreadBetweenGroups :: RawMatrixList Double -> RawVector Double +totalSpreadBetweenGroups xss = [s y and |i <- [0..((count avgNormDisks)-1)], let y = ys !! i, let and = avgNormDisks !! i] + where + s y and = sum [i*(^2) w| g <- [0..((count is)-1)], let i = is !! g, let w = ((y !! g)-and ) ] + is = countRows xss + avgNormDisks = averageOfScaledDiscriminants xss + ys = centroids xss + +-- | Calculates the total spread within the groups. +totalSpreadWithinGroups :: RawMatrixList Double -> RawVector Double +totalSpreadWithinGroups xss = [s y nd | i <- [0..((count normDisks)-1)], let y = ys !! i, let nd = normDisks !! i] + where + s y and = sum $ map sum [map (\x-> (^2) (x-yy)) gnd | i <- [0..((count y)-1)], let yy = y !! i, let gnd = and !! i] + normDisks = scaledDiscriminants xss + ys = centroids xss + +-- | Calculates the discriminant criteria. +discriminantCriteria :: RawMatrixList Double -> RawVector Double +discriminantCriteria xss = [ssb/ssw | i <- [0..((count ssbs)-1)], let ssb = ssbs !! i, let ssw = ssws !! i] + where + ssbs = totalSpreadBetweenGroups xss + ssws = totalSpreadWithinGroups xss + +-- | Calculation of the a priori probability, more precisely the probability that an element belongs to a group. +aprioriProbability :: RawMatrixList Double -> RawVector Double +aprioriProbability xsss = [ i / i_ | g <- [0..g_], let i = is !! round g] + where + g_ = countMatrixes xsss-1 + is = countRows xsss + i_ = countAllRows xsss + +-- | Calculates the constant part of the classification function according to Fisher. +fisherClassificationFunctionConst :: RawMatrixList Double -> RawVector Double +fisherClassificationFunctionConst xsss = [ -0.5 * s + log p | g' <- [0..g_], let g = round g', let s = (ss g), let p = ps !! g] + where + ss g = sum [b * x | j' <- [0..j_], let j = round j', let b = bg !! g !! j, let x = x_ @@> (g,j)] + g_ = countMatrixes xsss - 1 + j_ = countMatrixesCols xsss -1 + ps = aprioriProbability xsss + bg = fisherClassificationFunctionVar xsss + x_ = fromLists.map (map (\x -> (foldl1 (+) x) / (count x) )) $ map transpose $ xsss + +-- | Calculates the variable parts of the classification function according to Fisher. +fisherClassificationFunctionVar :: RawMatrixList Double -> RawVector (LinFunction Double) +fisherClassificationFunctionVar xsss = [ b g | g <- [0..g_] ] + where + b g = [ig * b_ (round j) (round g) | j <- [0..j_] ] + b_ j g = sum [ (w * x) | rr <- [0..j_], let r = round rr, let w = w' @@> (r,j), let x = x_ @@> (g,r) ] + ig = i_ - g_ + is = countRows xsss + i_ = sum is -1 + g_ = countMatrixes xsss -1 + j_ = countMatrixesCols xsss -1 + w' = inv.trans.spreadWithinGroups $ xsss + x_ = fromLists.map (map (\x -> (foldl1 (+) x) / (count x) )) $ map transpose $ xsss + +-- | Calculates the classification function according to Fisher. +fisherClassificationFunction :: RawMatrixList Double -> RawVector (LinFunction Double) +fisherClassificationFunction xsss = zipWith (:) ( fisherClassificationFunctionConst xsss) (fisherClassificationFunctionVar xsss) + +-- | Calculation of the classification of a survey (or attributes) in a cluster. The function takes a vector/list of attributes/values and a context. The context consists of groups/clusters and its items values/attributes. The function returns the ID (starting with 0) of the cluster to which the given vector/list belongs to. This function uses the Fisher algorithm. +fisher :: RawMatrixList Double -> RawVector Double -> Int +fisher xsss attributes = fisher' (fisherClassificationFunction xsss) attributes + +-- | Calculates the ID of the cluster the given values belonging to. This function takes a list of clusters, representated by a tuple, and a list of values. The cluster-tuples consists of a ID of the cluster and the classification function (according to Fisher) of the cluster. This function uses the Fisher algorithm. +fisherT :: RawVector (Int, LinFunction Double) -> RawVector Double -> Int +fisherT clusterTupels obj = fst $ clusterTupels !! (fisher' (map snd clusterTupels) obj) + +-- | Calculates the ID (starting with 0) of the cluster the given list of attributes belongs to. The function takes a list of attributes and a list of clusters which are representated by there classification function. This function uses the Fisher algorithm. +fisher' :: RawVector (LinFunction Double) -> RawVector Double -> Int +fisher' clusterFunctions obj = maxPos $ map (flip calcLinFunction $ obj) clusterFunctions + +-- | Calculates the cluster of every survey of a poll. This function takes the data of a whole poll and classifies every survey of the poll. This function uses the Fisher algorithm. +fisherAll :: RawMatrixList Double -> RawMatrix Int +fisherAll xsss = foldVectors (\a -> fisher xsss a) xsss diff --git a/src/Numeric/Vector.hs b/src/Numeric/Vector.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Vector.hs @@ -0,0 +1,58 @@ +-- While working on this module you are encouraged to remove it and fix +-- any warnings in the module. See +-- http://hackage.haskell.org/trac/ghc/wiki/WorkingConventions#Warnings +-- for details + +----------------------------------------------------------------------------- +-- | +-- Module : Vector +-- Copyright : (c) Lennart Schmitt +-- License : BSD-style (see the file libraries/base/LICENSE) +-- +-- Maintainer : libraries@haskell.org +-- Stability : experimental +-- Portability : portable +-- +-- This module implements a few extensions for the vector-module. +----------------------------------------------------------------------------- + +module Numeric.Vector ( + -- * Data Types + RawVector, + Vector, + + -- * Functions + average, + count, + fromList,-- | Convertes the representation from a simple list to a vector + maximumBy, -- | Calculates a lists maximum depending on a given ordering-function + maxPos, + toList, -- | Convertes the representation from a vector to a simple list + transpose -- | Transposes a Vector + ) where + +import Data.Packed.Vector (Vector, toList, fromList) +import Data.List (transpose, maximumBy) + +-- | A Vector represented by a simple list. +type RawVector a = [a] + +-- | Calculates the lists elements average +-- +-- > average [1,3,2] == 2.0 +average :: Floating a => RawVector a -> a +average xs = (sum xs) / (count xs) + +-- | Calculates the position of a lists maximum +-- +-- > maxPos [1,10,8,3] == 1 +maxPos :: RawVector Double -> Int +maxPos xs = fst $ foldl1 f [(i,x)|(i,x)<- zip [0..] xs] + where + f a b = if snd a > snd b then a else b + +-- | Counts the elements of a given list +-- +-- > count [1,2,3,4,5] == 5 +count :: Num b => RawVector a -> b +count = sum . map (const 1) diff --git a/src/Tests.hs b/src/Tests.hs new file mode 100644 --- /dev/null +++ b/src/Tests.hs @@ -0,0 +1,82 @@ +module Tests where + +import Numeric.Vector +import Numeric.Matrix +import Numeric.MatrixList +import Numeric.Function +import Numeric.Statistics.LDA +import Test.HUnit ( + Test (..), + assertEqual, + runTestTT + ) + +{----------------------------------------------------------------- +-------------------- Init Testdata ------------------------------- +-----------------------------------------------------------------} +testMatrix :: [Matrix Double] +testMatrix = [testMatrix_A,testMatrix_B] + +testMatrix_A :: Matrix Double +testMatrix_A = fromLists testdata_A + +testMatrix_B :: Matrix Double +testMatrix_B = fromLists testdata_B + +testdata :: RawMatrixList Double +testdata = [testdata_A,testdata_B] + +-- [[1_Attribute1, ..., 1_AttributeN], ..., [N_Attribute1, ..., N_AttributeN]] +testdata_A :: RawMatrix Double +testdata_A = [[2,3],[3,4],[6,5],[4,4],[3,2],[4,7],[3,5],[2,4],[5,6],[3,6],[3,3],[4,5]] + +testdata_B :: RawMatrix Double +testdata_B = [[5,4],[4,3],[7,5],[3,3],[4,4],[5,2],[4,2],[5,5],[6,7],[5,3],[6,4],[6,6]] + +{----------------------------------------------------------------- +-------------------- Testing Matrix Module ----------------------- +-----------------------------------------------------------------} +test_eigenvalue = TestCase ( assertEqual "Testing eigenvalue-function" 0.9122456375784809 ( eigenvalue (fromLists [[0.77143,-0.25714],[-0.42245,0.14082]]) ) ) + +test_eigenvector = TestCase ( assertEqual "Testing eigenvector-function" [0.8770950095147589,-0.48031692067249215] (toList.eigenvector $ (fromLists [[-0.14081563757848092,-0.25714],[-0.42245,-0.7714256375784809]])) ) + +test_reduceMatrix = TestCase ( assertEqual "Testing reduceMatrix-function" ([[ -0.14081563757848092, -0.25714], [-0.42245, -0.7714256375784809 ]]) ( Numeric.Matrix.toLists.reduceMatrix $ (fromLists [[0.77143,-0.25714],[-0.42245,0.14082]]) ) ) + + +{----------------------------------------------------------------- +-------------------- Testing Util Module ------------------------- +-----------------------------------------------------------------} + +test_averages = TestCase ( assertEqual "Testing averages-function" [[2.5,3.5,5.5,4.0,2.5,5.5,4.0,3.0,5.5,4.5,3.0,4.5],[4.5,3.5,6.0,3.0,4.0,3.5,3.0,5.0,6.5,4.0,5.0,6.0]] + (averages testdata)) + +{----------------------------------------------------------------- +-------------------- Testing LDA Module -------------------------- +-----------------------------------------------------------------} + +test_fisher = TestCase ( assertEqual "Testing fisher-function" 0 (fisher testdata [2,3])) + +test_fisherAll = TestCase ( assertEqual "Testing fisherAll-function" [[0,0,1,0,0,0,0,0,0,0,0,0],[1,1,1,0,0,1,1,1,1,1,1,1]] (fisherAll testdata)) + +test_fisherClassificationFunction = TestCase ( assertEqual "Testing fisherClassificationFunction-function" [[-6.5972288132130075,1.7285714285714286,1.279591836734694],[-10.222739017294638,3.614285714285714,0.246938775510204]] (fisherClassificationFunction testdata)) + +test_aprioriProbability = TestCase ( assertEqual "Testing aprioriProbability-function" [0.5,0.5] (aprioriProbability testdata)) + +test_discriminantCriteria = TestCase ( assertEqual "Testing discriminantCriteria-function" [0.9122448979591838,0.9122448979591834] (discriminantCriteria testdata)) + +test_isolatedDiscriminant = TestCase ( assertEqual "Testing isolatedDiscriminant-function" [0.46551724137931033,3.0612244897959183e-2] (isolatedDiscriminant testdata)) + + +tests = TestList [ + TestLabel "Test eigenvalue" test_eigenvalue, + TestLabel "Test eigenvector" test_eigenvector, + TestLabel "Test reduceMatrix" test_reduceMatrix, + TestLabel "Test averages" test_averages, + TestLabel "Test fisher" test_fisher, + TestLabel "Test fisherAll" test_fisherAll, + TestLabel "Test fisherClassificationFunction" test_fisherClassificationFunction, + TestLabel "Test aprioriProbability" test_aprioriProbability, + TestLabel "Test discriminantCriteria" test_discriminantCriteria, + TestLabel "Test isolatedDiscriminant"test_isolatedDiscriminant] + +main = runTestTT tests