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
diffAverage :: RawMatrixList Double -> RawMatrixList Double
diffAverage xs = ( zipWith $ zipWith (\ _xs y -> map (\x -> ()x y) _xs ) ) xs (averages xs)
squareDiffAverage :: RawMatrixList Double -> RawMatrixList Double
squareDiffAverage = ( mapElements (^2) ) . diffAverage
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 ) )
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 )))
spreadWithinGroups :: RawMatrixList Double -> Matrix Double
spreadWithinGroups = fromListToQuadraticMatrix . (zipAllWith sum) . sumOfAverages . transposeAll
spreadFromTotalAverages :: RawMatrixList Double -> RawMatrix Double
spreadFromTotalAverages x = map (zipWith (flip()) xs) xss
where
xs = (totalAverages x)
xss = (averages.transposeAll $ x)
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
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
calcA :: RawMatrixList Double -> Matrix Double
calcA xss = cA (spreadWithinGroups xss) (spreadBetweenGroups xss)
where
cA w b = multiplyR (inv w) b
scalingFactor :: RawMatrixList Double -> Double
scalingFactor xsss = 1 / s
where
s = sqrt (v'Wv / (fallzahlgruppen))
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
scaledDiscriminantCoefficient :: Matrix Double -> Double -> Matrix Double
scaledDiscriminantCoefficient v s = scalarMultiplication s v
constElement :: Matrix Double -> RawVector Double -> Double
constElement b x = 1 * (sum $ zipWith (*) (toList . flatten $ b) x)
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])
scaledDiscriminant :: LinFunction Double -> Values Double -> Double
scaledDiscriminant v x = calcLinFunction v x
scaledDiscriminants :: RawMatrixList Double -> RawMatrixList Double
scaledDiscriminants xss = map (\linFun -> foldVectors (scaledDiscriminant linFun) xss) normDiskFs
where
normDiskFs = scaledDiscriminantFunctions xss
centroid :: LinFunction Double -> RawVector(Values Double) -> Double
centroid v m = (sum $ map (scaledDiscriminant v ) m) / (count m)
centroids :: RawMatrixList Double -> RawMatrix Double
centroids xss = map (\linFun -> map (centroid linFun) xss) normDiskFs
where
normDiskFs = scaledDiscriminantFunctions xss
averageOfScaledDiscriminants :: RawMatrixList Double -> RawVector Double
averageOfScaledDiscriminants xss = map (\x -> (foldl1 (+) x) / (count x) ) $ map concat $ scaledDiscriminants xss
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
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) (xyy)) gnd | i <- [0..((count y)1)], let yy = y !! i, let gnd = and !! i]
normDisks = scaledDiscriminants xss
ys = centroids xss
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
aprioriProbability :: RawMatrixList Double -> RawVector Double
aprioriProbability xsss = [ i / i_ | g <- [0..g_], let i = is !! round g]
where
g_ = countMatrixes xsss1
is = countRows xsss
i_ = countAllRows xsss
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
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
fisherClassificationFunction :: RawMatrixList Double -> RawVector (LinFunction Double)
fisherClassificationFunction xsss = zipWith (:) ( fisherClassificationFunctionConst xsss) (fisherClassificationFunctionVar xsss)
fisher :: RawMatrixList Double -> RawVector Double -> Int
fisher xsss attributes = fisher' (fisherClassificationFunction xsss) attributes
fisherT :: RawVector (Int, LinFunction Double) -> RawVector Double -> Int
fisherT clusterTupels obj = fst $ clusterTupels !! (fisher' (map snd clusterTupels) obj)
fisher' :: RawVector (LinFunction Double) -> RawVector Double -> Int
fisher' clusterFunctions obj = maxPos $ map (flip calcLinFunction $ obj) clusterFunctions
fisherAll :: RawMatrixList Double -> RawMatrix Int
fisherAll xsss = foldVectors (\a -> fisher xsss a) xsss