summaryrefslogtreecommitdiff
path: root/examples/pca2.hs
diff options
context:
space:
mode:
Diffstat (limited to 'examples/pca2.hs')
-rw-r--r--examples/pca2.hs65
1 files changed, 0 insertions, 65 deletions
diff --git a/examples/pca2.hs b/examples/pca2.hs
deleted file mode 100644
index e7ea95f..0000000
--- a/examples/pca2.hs
+++ /dev/null
@@ -1,65 +0,0 @@
1-- Improved PCA, including illustrative graphics
2
3import Numeric.LinearAlgebra
4import Graphics.Plot
5import System.Directory(doesFileExist)
6import System.Process(system)
7import Control.Monad(when)
8
9type Vec = Vector Double
10type Mat = Matrix Double
11
12-- Vector with the mean value of the columns of a matrix
13mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
14
15-- covariance matrix of a list of observations stored as rows
16cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
17 where xc = x - asRow (mean x)
18
19
20type Stat = (Vec, [Double], Mat)
21-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)
22stat :: Mat -> Stat
23stat x = (m, toList s, trans v) where
24 m = mean x
25 (s,v) = eigSH' (cov x)
26
27-- creates the compression and decompression functions from the desired reconstruction
28-- quality and the statistics of a data set
29pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec)
30pca prec (m,s,v) = (encode,decode)
31 where
32 encode x = vp <> (x - m)
33 decode x = x <> vp + m
34 vp = takeRows n v
35 n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)
36 cumSum = tail . scanl (+) 0.0
37 prec' = if prec <=0.0 || prec >= 1.0
38 then error "the precision in pca must be 0<prec<1"
39 else prec
40
41shdigit :: Vec -> IO ()
42shdigit v = imshow (reshape 28 (-v))
43
44-- shows the effect of a given reconstruction quality on a test vector
45test :: Stat -> Double -> Vec -> IO ()
46test st prec x = do
47 let (pe,pd) = pca prec st
48 let y = pe x
49 print $ dim y
50 shdigit (pd y)
51
52main = do
53 ok <- doesFileExist ("mnist.txt")
54 when (not ok) $ do
55 putStrLn "\nTrying to download test datafile..."
56 system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
57 system("gunzip mnist.txt.gz")
58 return ()
59 m <- loadMatrix "mnist.txt"
60 let xs = takeColumns (cols m -1) m
61 let x = toRows xs !! 4 -- an arbitrary test vector
62 shdigit x
63 let st = stat xs
64 test st 0.90 x
65 test st 0.50 x