diff --git a/Learning.cabal b/Learning.cabal --- a/Learning.cabal +++ b/Learning.cabal @@ -2,13 +2,13 @@ -- -- see: https://github.com/sol/hpack -- --- hash: 41211cf12c83c4bc7d6b1152f8d6adc0e35b8d8a8aee00c41bf06728458510ac +-- hash: 6c4606a45d47a42344ba884c903f3f574904cf74f1bf8ba2b085dc33882900dd name: Learning -version: 0.0.1 -synopsis: Most frequently used machine learning tools +version: 0.0.2 +synopsis: The most frequently used machine learning tools description: Please see the README on Github at -category: ML +category: Machine Learning homepage: https://github.com/masterdezign/Learning#readme bug-reports: https://github.com/masterdezign/Learning/issues author: Bogdan Penkovsky @@ -32,6 +32,7 @@ src build-depends: base >=4.7 && <5 + , containers , hmatrix >=0.18.0.0 , vector exposed-modules: @@ -40,20 +41,38 @@ Paths_Learning default-language: Haskell2010 -executable Learning-exe - main-is: Main.hs +executable learning-pca + main-is: MainPCA.lhs hs-source-dirs: app ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends: Learning , base >=4.7 && <5 + , containers , hmatrix >=0.18.0.0 , vector other-modules: + MainPCA2 Paths_Learning default-language: Haskell2010 +executable learning-pca-advanced + main-is: MainPCA2.lhs + hs-source-dirs: + app + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + Learning + , base >=4.7 && <5 + , containers + , hmatrix >=0.18.0.0 + , vector + other-modules: + MainPCA + Paths_Learning + default-language: Haskell2010 + test-suite Learning-test type: exitcode-stdio-1.0 main-is: Spec.hs @@ -63,6 +82,7 @@ build-depends: Learning , base >=4.7 && <5 + , containers , hmatrix >=0.18.0.0 , vector other-modules: diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,4 +1,44 @@ # Learning -A micro library containing the most common machine learning tools -written in Haskell. +A Haskell micro library containing the most common machine learning tools. + +The name of the package can be interpreted in two ways: + +1. Either as "Learning" in "Machine Learning". +2. Or "Learning" meaning that [examples](https://github.com/masterdezign/Learning/tree/master/app) +are written in [literate style](https://en.wikipedia.org/wiki/Literate_programming) +and can be used to discover machine learning techniques. + + +## Features + +* Supervised learning + * Ridge regression + * Linear classifier +* Evaluation metrics +* Principal components analysis + + +## Getting Started + +Use [Stack](http://haskellstack.org). + + $ git clone https://github.com/masterdezign/Learning.git && cd Learning + $ stack build --install-ghc + +### Demo 1: principal components analysis (PCA) + +Launch the [PCA demo](https://github.com/masterdezign/Learning/blob/master/app/MainPCA.lhs) + + $ stack exec learning-pca + +### Demo 2: advanced principal components analysis (PCA) + +Launch the advanced [PCA demo](https://github.com/masterdezign/Learning/blob/master/app/MainPCA2.lhs) + + $ stack exec learning-pca-advanced + +### What's next? + +Check the [documentation](https://hackage.haskell.org/package/Learning/docs/Learning.html) +or [open an issue](https://github.com/masterdezign/Learning/issues). diff --git a/app/Main.hs b/app/Main.hs deleted file mode 100644 --- a/app/Main.hs +++ /dev/null @@ -1,4 +0,0 @@ -module Main where - -main :: IO () -main = putStrLn "No demo yet" diff --git a/app/MainPCA.lhs b/app/MainPCA.lhs new file mode 100644 --- /dev/null +++ b/app/MainPCA.lhs @@ -0,0 +1,62 @@ +Principal Components Analysis (PCA) demo +---------------------------------------- + + +The tutorial is based on http://setosa.io/ev/principal-component-analysis/ + +Suppose, we study nutrition habits of the citizens of four countries. +Here, we provide the food consumption data among those countries. + +> import Learning +> import qualified Numeric.LinearAlgebra as LA + +> england = [375, 57, 245, 1472, 105, 54, 193, 147, 1102, +> 720, 253, 685, 488, 198, 360, 1374, 156] + +> northernIreland = [135, 47, 267, 1494, 66, 41, 209, 93, 674, +> 1033, 143, 586, 355, 187, 334, 1506, 139] + +> scotland = [458, 53, 242, 1462, 103, 62, 184, 122, 957, +> 566, 171, 750, 418, 220, 337, 1572, 147] + +> wales = [475, 73, 227, 1582, 103, 64, 235, 160, 1137, +> 874, 265, 803, 570, 203, 365, 1256, 175] + +We want to know how differ the countries based on those data. +For that purpose, we would like to reduce the redundant information +or, in other words, perform PCA. + +We create a single list of feature vectors (each country) used later +for the analysis. + +> countries = map LA.fromList [england, northernIreland, scotland, wales] + +We perform PCA, i.e. calculate the compression (dimensionality reduction) +function `compress`. The `pca` function is given +`principalComponents` parameter. Here it's 1, that means that +`countries` vectors of 17 features will be reduced into scalars (1D vectors). + +> compress = let principalComponents = 1 +> pca1 = pca principalComponents countries +> in _compress pca1 + +Output the resulting scalar values for each country + +> main = mapM_ (print. compress) countries + +Here is a summary: + + England -702.9850482521952 + Northern Ireland -80.6002572540017 + Scotland -649.8612350689822 + Wales -798.521043705295 + +Wales + \/ England Northern Ireland + \/ \/ +..o....o..o.................................o + /\ + Scotland + +Now, we can clearly see that there exists a difference between +Northern Ireland and the rest of the countries. diff --git a/app/MainPCA2.lhs b/app/MainPCA2.lhs new file mode 100644 --- /dev/null +++ b/app/MainPCA2.lhs @@ -0,0 +1,77 @@ +Advanced Principal Components Analysis (PCA) demo +------------------------------------------------- + + +The tutorial is a continuation of ./MainPCA.lhs + +Previously, we were able to quickly determine that the food ration +in Northern Ireland is somewhat different comparing to the other three +countries. We have projected data from 17 dimensions into +one dimension. However, we could also make a projection into +two or more dimensions. So how do we determine which dimensionality +reduction does preserve the most of the information? How do we make +sure that only redundant information was removed? + +For that purpose, we will calculate retained variance [1] +depending on the number of the principal components. + +[1] http://www.dsc.ufcg.edu.br/~hmg/disciplinas/posgraduacao/rn-copin-2014.3/material/SignalProcPCA.pdf + +Let's start with the imports and data definition. + +> import Learning ( pca' ) +> import qualified Numeric.LinearAlgebra as LA +> import Data.List ( scanl' ) +> import Text.Printf ( printf ) + +> england = [375, 57, 245, 1472, 105, 54, 193, 147, 1102, +> 720, 253, 685, 488, 198, 360, 1374, 156] + +> northernIreland = [135, 47, 267, 1494, 66, 41, 209, 93, 674, +> 1033, 143, 586, 355, 187, 334, 1506, 139] + +> scotland = [458, 53, 242, 1462, 103, 62, 184, 122, 957, +> 566, 171, 750, 418, 220, 337, 1572, 147] + +> wales = [475, 73, 227, 1582, 103, 64, 235, 160, 1137, +> 874, 265, 803, 570, 203, 365, 1256, 175] + +> countries = map LA.fromList [england, northernIreland, scotland, wales] + +In order to compute the eigenvectors u and eigenvalues eig of +a covariance matrix, we use function pca'. In fact, pca' was +already called under the hood in the previous tutorial. + +> (u, eig) = pca' countries + +Sums of the first N eigenvalues: + +> cumul = drop 1 $ scanl' (+) 0 $ LA.toList eig + +Sum of all eigenvalues: + +> total = last cumul + +> main = mapM_ (\(i, s) -> +> let retained = s / total * 100 :: Double +> msg = "%d principal component(s): Retained variance %.1f%%" +> in putStrLn $ printf msg i retained) +> $ zip [1::Int ..] cumul + + 1 principal component(s): Retained variance 67.4% + 2 principal component(s): Retained variance 96.5% + 3 principal component(s): Retained variance 100.0% + 4 principal component(s): Retained variance 100.0% + + ... + + 17 principal component(s): Retained variance 100.0% + +From this data we can conclude that the first two principal components +contain 96.5% of information. Therefore, we will loose 3.5% of information +after projecting into two orthogonal axes in the transformed coordinate system +obtained after PCA. + +Hint: to compute compression (dimensionality reduction) and +decompression functions for specified variance to retain, use +(_compress. pcaVariance) and (_decompress. pcaVariance) functions. diff --git a/src/Learning.hs b/src/Learning.hs --- a/src/Learning.hs +++ b/src/Learning.hs @@ -14,6 +14,7 @@ , PCA (..) , pca , pca' + , pcaVariance -- * Supervised learning , Teacher @@ -28,15 +29,26 @@ , winnerTakesAll -- * Evaluation + -- ** Classification + , accuracy , errorRate , errors - , accuracy + , showConfusion + , confusion + , Normalize (..) + , confusion' + + -- ** Regression , nrmse ) where import Numeric.LinearAlgebra import qualified Data.Vector.Storable as V +import qualified Data.Map as M +import Data.List ( nub, sort ) +import Text.Printf ( printf ) + -- | A dataset representation for supervised learning data Dataset a b = Dataset { _samples :: [a] @@ -65,9 +77,8 @@ -> (Matrix Double, Vector Double) pca' xs = (u', s) where - xs' = fromBlocks $ map ((: []). tr. reshape 1) xs -- Covariance matrix - sigma = snd $ meanCov xs' + sigma = snd $ meanCov $ fromRows xs -- Eigenvectors matrix u' and eigenvalues vector s (u', s, _) = svd $ unSym sigma @@ -83,18 +94,35 @@ -- | Principal components analysis resulting in `PCA` tools pca :: Int -- ^ Number of principal components to preserve - -> [Vector Double] -- ^ Analyzed data samples + -> [Vector Double] -- ^ Observations -> PCA pca maxDim xs = let (u', _) = pca' xs - u = takeColumns maxDim u' - in PCA - { _u = u - , _compress = (tr u <>). reshape 1 - , _decompress = flatten. (u <>) - } + in _pca maxDim u' --- | Classifier function that maps some network state with measurements as matrix columns --- and features as rows, into a categorical output. +-- | Perform PCA using the minimal number of principal +-- components required to retain given variance +pcaVariance :: Double -- ^ Retained variance, % + -> [Vector Double] -- ^ Observations + -> PCA +pcaVariance var xs = let (u', eig) = pca' xs + cumul = V.drop 1 $ V.scanl' (+) 0 eig + var' = var / 100 -- Scale 100% -> 1.0 + total = V.last cumul + isRetained = map (\v -> let retained = v / total + in retained >= var') $ V.toList cumul + dim = fst $ head $ filter snd $ zip [1..] isRetained + in _pca dim u' + +_pca :: Int -> Matrix Double -> PCA +_pca maxDim u' = let u = takeColumns maxDim u' + in PCA + { _u = u + , _compress = (tr u <>). reshape 1 + , _decompress = flatten. (u <>) + } + +-- | Classifier function that maps some measurements as matrix columns +-- and corresponding features as rows, into a categorical output. newtype Classifier a = Classifier { classify :: Matrix Double -> a } -- | Regressor function that maps some feature matrix @@ -157,7 +185,7 @@ -- Similar to `learnRegressor`, but instead of a `Regressor` function -- a (already transposed) `Readout` matrix may be returned. learn' - :: Matrix Double -- ^ Network state (nonlinear response) + :: Matrix Double -- ^ Measurements (feature matrix) -> Matrix Double -- ^ Horizontally concatenated `Teacher` matrices -> Maybe Readout learn' a b = case ridgeRegression 1e-4 a b of @@ -230,20 +258,110 @@ where errNo = length $ errors $ zip tgtLbls cLbls {-# SPECIALIZE errorRate :: [Int] → [Int] → Double #-} --- | Accuracy of classification, @100% - errorRate@ +-- | Accuracy of classification, @100% - `errorRate`@ -- -- >>> accuracy [1,2,3,4] [1,2,3,7] -- 75.0 -accuracy :: (Eq a, Fractional acc) => [a] -> [a] -> acc +accuracy :: (Eq lab, Fractional acc) => [lab] -> [lab] -> acc accuracy tgt clf = let erate = errorRate tgt clf in 100 - erate {-# SPECIALIZE accuracy :: [Int] → [Int] → Double #-} +-- | Confusion matrix for arbitrary number of classes (not normalized) +confusion' :: (Ord lab, Eq lab) + => [lab] + -- ^ Target labels + -> [lab] + -- ^ Predicted labels + -> M.Map (lab, lab) Int + -- ^ Map keys: (target, predicted), values: confusion count +confusion' tgtlab lab = mp + where + -- Count all possible pairs of labels + mp = foldr g M.empty $ zip tgtlab lab + g k mp = M.alter f k mp + + f Nothing = Just 1 + f (Just x) = Just (x + 1) + +-- | Normalization strategies for `confusion` matrix +data Normalize = ByRow | ByColumn deriving (Show, Eq) + +-- | Normalized confusion matrix for arbitrary number of classes +confusion :: (Ord lab, Eq lab) + => Normalize + -- ^ Normalize `ByRow` or `ByColumn` + -> [lab] + -- ^ Target labels + -> [lab] + -- ^ Predicted labels + -> M.Map (lab, lab) Double + -- ^ Map keys: (target, predicted), values: normalized confusion +confusion by tgtlab lab = res + where + allLabels = sort $ nub tgtlab + mp = confusion' tgtlab lab + lookup2 k' mp' = case M.lookup k' mp' of + Just x -> x + _ -> 0 + + res = foldr (\i mp' -> let key j = if by == ByRow + then (i, j) + else (j, i) + -- Find sum + grp = map (\j -> let k = key j in (k, lookup2 k mp)) allLabels + total = fromIntegral $ sum $ map snd grp + -- Normalize + in foldr (\(k, v) mp'' -> M.insert k (fromIntegral v / total * 100) mp'') mp' grp + ) M.empty allLabels + +-- | Confusion matrix normalized by row: ASCII representation. +-- +-- Note: it is assumed that target (true) labels list contains +-- all possible labels. +-- +-- @ +-- | Predicted +-- ---+------------ +-- | \_ \_ \_ \_ \_ +-- True | \_ \_ \_ \_ \_ +-- | \_ \_ \_ \_ \_ +-- label | \_ \_ \_ \_ \_ +-- | \_ \_ \_ \_ \_ +-- @ +-- +-- >>> putStr $ showConfusion [1, 2, 3, 1] [1, 2, 3, 2] +-- 1 2 3 +-- 1 50.0 50.0 0.0 +-- 2 0.0 100.0 0.0 +-- 3 0.0 0.0 100.0 +showConfusion :: (Ord lab, Eq lab, Show lab) + => [lab] -- ^ Target labels + -> [lab] -- ^ Predicted labels + -> String +showConfusion tgtlab lab = unlines $ predictedLabels: "": table + where + allLabels = sort $ nub tgtlab + mp = confusion ByRow tgtlab lab + table = map (fmtRow mp) allLabels + + predictedLabels = let spc1 = replicate 2 ' ' + spc2 = replicate 4 ' ' + in spc1 ++ (unwords $ map ((spc2 ++). show) allLabels) + + -- Tabulate row + fmtRow mp i = unwords (show i: "": line) + where + fmt x = let s = printf "%.1f" x + l = length s + in replicate (5 - l) ' ' ++ s + line = map (\j -> fmt $ mp M.! (i, j)) allLabels + -- | Pairs of misclassified and correct values -- -- >>> errors $ zip ['x','y','z'] ['x','b','a'] -- [('y','b'),('z','a')] -errors :: Eq a => [(a, a)] -> [(a, a)] +errors :: Eq lab => [(lab, lab)] -> [(lab, lab)] errors = filter (uncurry (/=)) {-# SPECIALIZE errors :: [(Int, Int)] -> [(Int, Int)] #-} @@ -254,8 +372,8 @@ cov :: (V.Storable a, Fractional a) => Vector a -> Vector a -> a cov xs ys = V.sum (V.zipWith (*) xs' ys') / fromIntegral (V.length xs') where - xs' = V.map (`subtract` (mean xs)) xs - ys' = V.map (`subtract` (mean ys)) ys + xs' = V.map (`subtract` mean xs) xs + ys' = V.map (`subtract` mean ys) ys {-# SPECIALISE cov :: Vector Double -> Vector Double -> Double #-} var :: (V.Storable a, Fractional a) => Vector a -> a