diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -94,3 +94,6 @@
 0.2.4.2:
 		Conrad Parker pointed out that hmatrix defaults to 
 		  Data.Vector.Storable
+
+0.2.5:
+		Added pcaN as requested by Marcel Ruegenberg
diff --git a/hstatistics.cabal b/hstatistics.cabal
--- a/hstatistics.cabal
+++ b/hstatistics.cabal
@@ -1,5 +1,5 @@
 Name:               hstatistics
-Version:            0.2.4.2
+Version:            0.2.5
 License:            BSD3
 License-file:       LICENSE
 Copyright:          (c) A.V.H. McPhail 2010, 2011, 2012
diff --git a/lib/Numeric/Statistics/PCA.hs b/lib/Numeric/Statistics/PCA.hs
--- a/lib/Numeric/Statistics/PCA.hs
+++ b/lib/Numeric/Statistics/PCA.hs
@@ -13,13 +13,15 @@
 -----------------------------------------------------------------------------
 
 module Numeric.Statistics.PCA (
-                               pca, pcaTransform, pcaReduce
+                               pca, pcaN, pcaTransform, pcaReduce
                           ) where
 
 
 -----------------------------------------------------------------------------
 
 import qualified Data.Array.IArray as I 
+import Data.List(sortBy)
+import Data.Ord(comparing)
 
 import Numeric.LinearAlgebra
 
@@ -29,7 +31,8 @@
 
 -----------------------------------------------------------------------------
 
--- | find the n principal components of multidimensional data
+-- | find the principal components of multidimensional data greater than
+--    the threshhold
 pca :: I.Array Int (Vector Double)    -- the data
     -> Double                         -- eigenvalue threshold
     -> Matrix Double
@@ -41,6 +44,20 @@
               v' = zip val vec
               v = filter (\(x,_) -> x > q) v'  -- keep only eigens > than parameter
           in fromColumns $ snd $ unzip v
+
+-- | find N greatest principal components of multidimensional data
+--    according to size of the eigenvalue
+pcaN :: I.Array Int (Vector Double)    -- the data
+     -> Int                            -- number of components to return
+     -> Matrix Double
+pcaN d n = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
+               cv = covarianceMatrix d'
+               (val',vec') = eigSH cv           -- the covariance matrix is real symmetric
+               val = toList val'
+               vec = toColumns vec'
+               v' = zip val vec
+               v = take n $ reverse $ sortBy (comparing fst) v'
+           in fromColumns $ snd $ unzip v
 
 -- | perform a PCA transform of the original data (remove mean)
 -- |     Final = M^T Data^T