diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -103,3 +103,6 @@
 
 0.2.5.2:
 		exposed pcaReduceN as requested by Tom Nielsen
+
+0.2.5.3:
+		updated for hmatrix 0.17
diff --git a/hstatistics.cabal b/hstatistics.cabal
--- a/hstatistics.cabal
+++ b/hstatistics.cabal
@@ -1,8 +1,8 @@
 Name:               hstatistics
-Version:            0.2.5.2
+Version:            0.2.5.3
 License:            BSD3
 License-file:       LICENSE
-Copyright:          (c) A.V.H. McPhail 2010, 2011, 2012, 2013
+Copyright:          (c) A.V.H. McPhail 2010, 2011, 2012, 2013, 2014
 Author:             Vivian McPhail
 Maintainer:         haskell.vivian.mcphail <at> gmail <dot> com
 Stability:          provisional
@@ -16,7 +16,7 @@
      .
      Feature requests, suggestions, and bug fixes welcome.
 Category:           Math, Statistics
-tested-with:        GHC ==7.6.3
+tested-with:        GHC ==7.10.2
 
 cabal-version:      >=1.8
 
@@ -30,8 +30,8 @@
     Build-Depends:      base >= 4 && < 5,
                         array, random,
                         vector,
-                        hmatrix >= 0.10.0.0,
-                        hmatrix-gsl-stats >= 0.1.2.9
+                        hmatrix >= 0.17,
+                        hmatrix-gsl-stats >= 0.4
 
     Extensions:         
 
@@ -51,6 +51,7 @@
     ghc-options:        -Wall -fno-warn-missing-signatures
                               -fno-warn-orphans
                               -fno-warn-unused-binds
+                        -O2
 
 source-repository head
     type:     git
diff --git a/lib/Numeric/Statistics.hs b/lib/Numeric/Statistics.hs
--- a/lib/Numeric/Statistics.hs
+++ b/lib/Numeric/Statistics.hs
@@ -2,8 +2,8 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics
--- Copyright   :  (c) A. V. H. McPhail 2010, 2012
--- License     :  BSD
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2012, 2014
+-- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
 -- Stability   :  provisional
@@ -31,14 +31,14 @@
 
 -----------------------------------------------------------------------------
 
---import Numeric.Vector
---import Numeric.Matrix
---import Numeric.Container
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra hiding(range)
+--import Numeric.LinearAlgebra.Data hiding(range)
+--import Numeric.LinearAlgebra.Devel
 
 import qualified Data.Array.IArray as I 
 import qualified Data.List as DL
 import qualified Data.Vector.Generic as GV
+--import qualified Data.Vector.Storable as SV
 
 import Foreign.Storable
 
@@ -68,36 +68,36 @@
 -----------------------------------------------------------------------------
 
 -- | the mean of a list of vectors
-meanList :: (Container Vector a, Num (Vector a)) => [Sample a] -> Sample a
+meanList :: (Container Vector a, Num (Vector a), Fractional a) => [Sample a] -> Sample a
 meanList []     = error "meanVectors: empty list"
 meanList [s]    = s
 meanList (s:ss) = let ln = fromIntegral $ length ss + 1
                   in scale (recip ln) $ foldl (+) s ss
 
 -- | the mean of an array of vectors
-meanArray :: (Container Vector a, Num (Vector a)) => Samples a -> Sample a
+meanArray :: (Container Vector a, Num (Vector a), Fractional a) => Samples a -> Sample a
 meanArray a = meanList $ I.elems a
 
 -- | the mean of a matrix with data series in rows
-meanMatrix :: (Container Vector a, Num (Vector a), Element a) => Matrix a -> Sample a
+meanMatrix :: (Container Vector a, Num (Vector a), Element a, Fractional a) => Matrix a -> Sample a
 meanMatrix a = meanList $ toRows a
 
 -----------------------------------------------------------------------------
 
 -- | the variance of a list of vectors
-varianceList :: (Container Vector a, Floating (Vector a)) => [Sample a] -> Sample a
+varianceList :: (Container Vector a, Floating (Vector a), Num a, Fractional a) => [Sample a] -> Sample a
 varianceList []  = error "varianceList: empty list"
-varianceList [s] = constant 0 (dim s)
+varianceList [s] = konst 0 (size s)
 varianceList l   = let mxs = meanList (map (** 2) l)
                        msx = (meanList l) ** 2
                    in mxs - msx
 
 -- | the variance of an array of vectors
-varianceArray :: (Container Vector a, Floating (Vector a)) => Samples a -> Sample a
+varianceArray :: (Container Vector a, Floating (Vector a), Num a, Fractional a) => Samples a -> Sample a
 varianceArray a = varianceList $ I.elems a
 
 -- | the variance of a matrix with data series in rows
-varianceMatrix :: (Container Vector a, Floating (Vector a), Element a) => Matrix a -> Sample a
+varianceMatrix :: (Container Vector a, Floating (Vector a), Element a, Num a, Fractional a) => Matrix a -> Sample a
 varianceMatrix a = varianceList $ toRows a
 
 -----------------------------------------------------------------------------
@@ -128,11 +128,11 @@
     -> Vector Double -- ^ intervals
     -> Vector Int    -- ^ data indexed by bin
 cut v c  = let c' = sort c
-           in mapVector (\x -> cut_helper 0 x c') v 
+           in GV.map (\x -> cut_helper 0 x c') v 
     where
       cut_helper j x d 
-          | j >= dim d                       = error "Numeric.Statistics: cut: data point not within interval"
-          | x >= (d @> j) && x <= (d @> (j+1)) = j
+          | j >= size d                       = error "Numeric.Statistics: cut: data point not within interval"
+          | x >= (d `atIndex` j) && x <= (d `atIndex` (j+1)) = j
           | otherwise                       = cut_helper (j + 1) x d
 
 -----------------------------------------------------------------------------
@@ -142,19 +142,19 @@
 --ranks :: Vector Double -> Vector Double
 ranks :: (Fractional b, Storable b) => Vector Double -> Vector b
 ranks v = let v' = sort v
-          in mapVector (\x -> 1 + rank_helper x v') v
+          in GV.map (\x -> 1 + rank_helper x v') v
               where rank_helper x v' = let is = GV.elemIndices x v'
-                                       in (realToFrac (GV.foldl (+) 0 is)) / (fromIntegral $ dim is)
+                                       in (realToFrac (GV.foldl (+) 0 is)) / (fromIntegral $ GV.length is)
 
 -----------------------------------------------------------------------------
 
 -- | kendall's rank correlation τ
 kendall :: Vector Double -> Vector Double -> Matrix Double
-kendall x y = let ln = dim x
+kendall x y = let ln = size x
                   rx = ranks x
                   ry = ranks y
                   r = fromColumns [rx,ry]
-                  m = signum $ (kronecker r (asColumn $ constant 1.0 ln)) - (kronecker (asRow $ constant 1.0 ln) r)
+                  m = signum $ (kronecker r (asColumn $ konst 1.0 ln)) - (kronecker (asRow $ konst 1.0 ln) r)
                   c = rows m - 1
               in correlationCoefficientMatrix $ I.listArray (0,c) (toColumns m)
 
@@ -164,7 +164,7 @@
 --logit :: Vector Double -> Vector Double
 logit :: (Floating b, Storable b)
         => Vector b -> Vector b
-logit v =  mapVector (\x -> - (log ((1 / x) - 1))) v
+logit v =  GV.map (\x -> - (log ((1 / x) - 1))) v
 
 -----------------------------------------------------------------------------
 
@@ -181,16 +181,16 @@
                                  Just m  -> m
                       xm     = fromRows $ map ((-) xu) $ toRows $ fromColumns xl
                       --um     = asColumn xu
-                      --w      = ((trans xm) <> xm + (trans um) <> um)/(fromIntegral $ xr - 1)
+                      --w      = ((tr' xm) <> xm + (tr' um) <> um)/(fromIntegral $ xr - 1)
                       --w'     = inv w
-                  in ((xm <> s' <> (trans xm)) @@> (0,0)) 
+                  in ((xm <> s' <> (tr' xm)) `atIndex` (0,0)) 
 
 -----------------------------------------------------------------------------
 
 -- | a list of element frequencies
 mode :: Vector Double -> [(Double,Integer)]
 mode v = let w = sort v
-         in DL.sortBy (\(_,n) (_,n') -> compare n' n) $ foldVector freqs [] w
+         in DL.sortBy (\(_,n) (_,n') -> compare n' n) $ GV.foldr freqs [] w
             where freqs x []          = [(x,1)]
                   freqs x ((f,n):fns)
                       | f == x         = ((f,n+1):fns) 
@@ -211,7 +211,7 @@
 --    | p == 2     = variance v -- gives sample variance
     | otherwise = let u = if c then centre v else v
                       w = if a then abs u else u
-                      x = mapVector (** (fromIntegral p)) w
+                      x = GV.map (** (fromIntegral p)) w
                   in mean x
 
 -----------------------------------------------------------------------------
@@ -227,13 +227,13 @@
     | rows x /= rows y = error "Numeric.Statistics: ols: incorrect matrix dimensions"
     | otherwise       = let (xr,xc) = (rows x,cols x)
                             (yr,yc) = (rows y,cols y)
-                            z = (trans x) <> x
+                            z = (tr' x) <> x
                             r = rank z
                             beta = if r == xc 
-                                      then (inv z) <> (trans x) <> y
+                                      then (inv z) <> (tr' x) <> y
                                       else (pinv x) <> y
                             rr = y - x <> beta
-                            sigma = ((trans rr) <> rr) / (fromIntegral $ xr - r)
+                            sigma = ((tr' rr) <> rr) / (fromIntegral $ xr - r)
                         in (beta,rr,sigma)
 
 -----------------------------------------------------------------------------
@@ -247,18 +247,18 @@
 -----------------------------------------------------------------------------
 
 -- | the difference between the maximum and minimum of the input
-range :: Container c e => c e -> e
+range :: (Container c e, Num e) => c e -> e
 range v = maxElement v - minElement v
 
 -----------------------------------------------------------------------------
 
 -- | count the number of runs greater than or equal to @n@ in the data
-run_count :: (Num a, Num t, Ord b, Ord a, Storable b) 
+run_count :: (Num a, Num t, Ord b, Ord a, Storable b, Container Vector b) 
             => a             -- ^ longest run to count
           -> Vector b        -- ^ data
           -> [(a, t)]        -- ^ [(run length,count)]
-run_count n v = let w = subVector 1 (dim v - 1) v
-                    x = foldVector run_count' [(1,v @> 0)] w
+run_count n v = let w = subVector 1 (size v - 1) v
+                    x = GV.foldr run_count' [(1,v `atIndex` 0)] w
                     y = map fst x
                     z = takeWhile (<= n) $  DL.sort y
                 in foldr count [] z
diff --git a/lib/Numeric/Statistics/Histogram.hs b/lib/Numeric/Statistics/Histogram.hs
--- a/lib/Numeric/Statistics/Histogram.hs
+++ b/lib/Numeric/Statistics/Histogram.hs
@@ -1,8 +1,8 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.Histogram
--- Copyright   :  (c) A. V. H. McPhail 2010
--- License     :  GPL-style
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
+-- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
 -- Stability   :  provisional
@@ -20,7 +20,8 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
+--import qualified Data.Packed.Vector as V
+import qualified Data.Vector.Storable as V
 
 import qualified Numeric.GSL.Histogram as H
 --import qualified Numeric.GSL.Histogram2D as H2
@@ -32,7 +33,7 @@
 
 -----------------------------------------------------------------------------
 
-vectorToTuples = toTuples . toList
+vectorToTuples = toTuples . V.toList
     where toTuples []         = error "need a minimum of two elements"
           toTuples [_]        = error "need a minimum of two elements"
           toTuples [x1,x2]    = [(x1,x2)]
@@ -40,14 +41,14 @@
 
 -----------------------------------------------------------------------------
 
-cumulativeToHistogram :: (Double -> Double)  -- ^ the cumulative distribution function D(x <= X)
-                   -> Vector Double          -- ^ the bins
+cumulativeToHistogram :: (Double -> Double)   -- ^ the cumulative distribution function D(x <= X)
+                   -> V.Vector Double        -- ^ the bins
                    -> H.Histogram            -- ^ the resulting histogram
 cumulativeToHistogram f v = H.addListWeighted (H.emptyRanges v) $ map (\(x1,x2) -> ((x1 + x2) / 2.0,f x2 - f x1)) (vectorToTuples v)
 
-gaussianHistogram :: Double              -- ^ mean
+gaussianHistogram :: Double               -- ^ mean
                   -> Double              -- ^ standard deviation
-                  -> Vector Double       -- ^ the bins
+                  -> V.Vector Double     -- ^ the bins
                   -> H.Histogram         -- ^ the resulting histogram
 gaussianHistogram u s = cumulativeToHistogram (\x -> C.density_1p C.Gaussian C.Lower s (x-u))
 
diff --git a/lib/Numeric/Statistics/ICA.hs b/lib/Numeric/Statistics/ICA.hs
--- a/lib/Numeric/Statistics/ICA.hs
+++ b/lib/Numeric/Statistics/ICA.hs
@@ -1,7 +1,8 @@
+{-# LANGUAGE FlexibleContexts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.ICA
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -23,7 +24,9 @@
 module Numeric.Statistics.ICA (
                                sigmoid, sigmoid',
                                demean, whiten,
-                               ica, icaDefaults
+                               ica, icaDefaults,
+                               --
+                               NormType(..)
                           ) where
 
 
@@ -33,6 +36,8 @@
 
 import Numeric.LinearAlgebra
 
+import qualified Data.Vector.Generic as GV
+
 import Numeric.GSL.Statistics
 
 import Numeric.Statistics
@@ -41,6 +46,16 @@
 
 -----------------------------------------------------------------------------
 
+data NormType = NormZero | NormOne | NormTwo | NormInf
+
+pnorm :: Normed (Vector a) => NormType -> Vector a -> R
+pnorm NormZero = norm_0
+pnorm NormOne  = norm_1
+pnorm NormTwo  = norm_2
+pnorm NormInf  = norm_Inf
+
+-----------------------------------------------------------------------------
+
 -- | sigmoid transfer function
 sigmoid :: Double -> Double
 sigmoid u = u * exp((-u**2)/2)
@@ -76,7 +91,7 @@
        -> Double                                      -- ^ eigenvalue threshold
        -> (I.Array Int (Vector Double),Matrix Double) -- ^ (whitened data,transform)
 whiten d q = let cv = covarianceMatrix d
-                 (val',vec') = eigSH cv           -- the covariance matrix is real symmetric
+                 (val',vec') = eigSH $ trustSym cv  -- the covariance matrix is real symmetric
                  val = toList val'
                  vec = toColumns vec'
                  v' = zip val vec
@@ -85,7 +100,7 @@
                  dd = diag $ (** (-0.5)) $ fromList dd'  -- square root of eigenvalues diagonalised
                  e = fromColumns e'
                  x = fromRows $ I.elems d
-                 t = e <> dd <> trans e          -- the actual mathematics
+                 t = e <> dd <> tr' e          -- the actual mathematics
                  x' = t <> x                     -- the actual mathematics
                  d' = I.listArray (I.bounds d) (toRows x') 
              in (d',t)
@@ -125,23 +140,23 @@
 update :: (Double -> Double) -> (Double -> Double) -> Matrix Double -> Matrix Double -> Matrix Double
 update g g' w x = let y = w <> x
                       ys = toRows y
-                      bis = map (\y' -> - mean (y' * (mapVector g y'))) ys
-                      ais = zipWith (\b y' -> -1 / (b - mean (mapVector g y'))) bis ys
+                      bis = map (\y' -> - mean (y' * (GV.map g y'))) ys
+                      ais = zipWith (\b y' -> -1 / (b - mean (GV.map g y'))) bis ys
                       r = rows y
                       ix = ((1,1),(r,r))
-                      cov = fromArray2D $ I.listArray ix $ map (\(m,n) -> covariance (mapVector g' (ys!!(m-1))) (ys!!(n-1))) $ I.range ix
+                      cov = fromArray2D $ I.listArray ix $ map (\(m,n) -> covariance (GV.map g' (ys!!(m-1))) (ys!!(n-1))) $ I.range ix
                   in w + (diag $ fromList ais) <> ((diag $ fromList bis) + cov) <> w  
 
 decorrelate :: Matrix Double -> Matrix Double
 decorrelate m = let (d',v') = eig m
                     d = fst $ fromComplex d'
                     v = fst $ fromComplex v'
-                in v <> (diag (d ** (-0.5))) <> trans v <> m
-{-decorrelate n t w = let w' = w / (scalar $ sqrt $ pnorm n (w <> trans w))
+                in v <> (diag (d ** (-0.5))) <> tr' v <> m
+{-decorrelate n t w = let w' = w / (scalar $ sqrt $ pnorm n (w <> tr' w))
                     in decorrelate' t w w'
     where decorrelate' t' m m' 
               | converged t' m m' = m'
-              | otherwise         = decorrelate' t' m' ((scale 1.5 m') - (scale 0.5 (m' <> trans m' <> m')))
+              | otherwise         = decorrelate' t' m' ((scale 1.5 m') - (scale 0.5 (m' <> tr' m' <> m')))
 -}
 
 normalise :: NormType -> Matrix Double -> Matrix Double
@@ -192,7 +207,7 @@
             -> I.Array Int (Vector Double) -- ^ data
             -> (I.Array Int (Vector Double),Matrix Double) -- ^ transformed data, ica transform
 icaDefaults r a = let c = I.rangeSize $ I.bounds a
-                      s = (dim $ (a I.! 1)) `div` 16
-                  in ica r sigmoid sigmoid' Infinity 0.0000001 s a
+                      s = (GV.length $ (a I.! 1)) `div` 16
+                  in ica r sigmoid sigmoid' NormInf 0.0000001 s a
 
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/Statistics/Information.hs b/lib/Numeric/Statistics/Information.hs
--- a/lib/Numeric/Statistics/Information.hs
+++ b/lib/Numeric/Statistics/Information.hs
@@ -2,7 +2,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.Information
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -23,15 +23,21 @@
 
 import Numeric.Statistics.PDF
 
-import Numeric.LinearAlgebra
+import qualified Numeric.LinearAlgebra as LA
+--import Numeric.LinearAlgebra.Data hiding(Vector)
 
+--import qualified Data.Vector as DV
+import qualified Data.Vector.Storable as V
+
+import Prelude hiding(map,zip)
+
 -----------------------------------------------------------------------------
 
 zeroToOne x
     | x == 0.0  = 1.0
     | otherwise = x
 
-logE = mapVector (log . zeroToOne)
+logE = V.map (log . zeroToOne)
 
 
 -----------------------------------------------------------------------------
@@ -39,21 +45,21 @@
 -- | the entropy \sum p_i l\ln{p_i} of a sequence
 entropy :: PDF a Double 
         => a                       -- ^ the underlying distribution
-        -> Vector Double           -- ^ the sequence
+        -> LA.Vector Double           -- ^ the sequence
         -> Double                  -- ^ the entropy
 entropy p x = let ps = probability p x
-              in negate $ (dot ps (logE ps))
+              in negate $ (LA.dot ps (logE ps))
 
 -- | the mutual information \sum_x \sum_y p(x,y) \ln{\frac{p(x,y)}{p(x)p(y)}}
 mutual_information :: (PDF a Double, PDF b (Double,Double)) 
                    => b                                          -- ^ the underlying distribution
                    -> a                                          -- ^ the first dimension distribution
                    -> a                                          -- ^ the second dimension distribution
-                   -> (Vector Double, Vector Double)             -- ^ the sequence
+                   -> (LA.Vector Double, LA.Vector Double)             -- ^ the sequence
                    -> Double         -- ^ the mutual information
-mutual_information p px py (x,y) = let ps = probability p $ zipVector x y
+mutual_information p px py (x,y) = let ps = probability p $ V.zipWith (,) x y
                                        xs = probability px x
                                        ys = probability py y
-                                   in (dot ps (logE ps - logE (xs*ys)))
+                                   in (LA.dot ps (logE ps - logE (xs*ys)))
 
 -----------------------------------------------------------------------------
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
@@ -1,7 +1,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.PCA
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -38,7 +38,7 @@
     -> Matrix Double
 pca d q = 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',vec') = eigSH $ trustSym cv -- the covariance matrix is real symmetric
               val = toList val'
               vec = toColumns vec'
               v' = zip val vec
@@ -52,7 +52,7 @@
      -> 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',vec') = eigSH $ trustSym cv  -- the covariance matrix is real symmetric
                val = toList val'
                vec = toColumns vec'
                v' = zip val vec
@@ -65,7 +65,7 @@
              -> Matrix Double                  -- ^ the principal components
              -> I.Array Int (Vector Double)    -- ^ the transformed data
 pcaTransform d m = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
-                   in I.listArray (1,cols m) $ toRows $ (trans m) <> (fromRows $ I.elems d')
+                   in I.listArray (1,cols m) $ toRows $ (tr' m) <> (fromRows $ I.elems d')
 
 -- | perform a dimension-reducing PCA modification, 
 --     using an eigenvalue threshhold
@@ -75,13 +75,13 @@
 pcaReduce d q = let u = fmap (scalar . mean) d
                     d' = zipWith (-) (I.elems d) (I.elems u)
                     cv = covarianceMatrix $ I.listArray (I.bounds d) d'
-                    (val',vec') = eigSH cv           -- the covariance matrix is real symmetric
+                    (val',vec') = eigSH $ trustSym cv -- the covariance matrix is real symmetric
                     val = toList val'
                     vec = toColumns vec'
                     v' = zip val vec
                     v = filter (\(x,_) -> x > q) v'  -- keep only eigens > than parameter
                     m = fromColumns $ snd $ unzip v
-                 in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (trans m) <> fromRows d') (I.elems u) 
+                 in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (tr' m) <> fromRows d') (I.elems u) 
 
 -- | perform a dimension-reducing PCA modification, using N components
 pcaReduceN :: I.Array Int (Vector Double)      -- ^ the data
@@ -90,12 +90,12 @@
 pcaReduceN d n = let u = fmap (scalar . mean) d
                      d' = zipWith (-) (I.elems d) (I.elems u)
                      cv = covarianceMatrix $ I.listArray (I.bounds d) d'
-                     (val',vec') = eigSH cv           -- the covariance matrix is real symmetric
+                     (val',vec') = eigSH $ trustSym 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'
                      m = fromColumns $ snd $ unzip v
-                  in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (trans m) <> fromRows d') (I.elems u) 
+                  in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (tr' m) <> fromRows d') (I.elems u) 
 
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/Statistics/PDF.hs b/lib/Numeric/Statistics/PDF.hs
--- a/lib/Numeric/Statistics/PDF.hs
+++ b/lib/Numeric/Statistics/PDF.hs
@@ -3,7 +3,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.PDF
--- Copyright   :  (c) A. V. H. McPhail 2010
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -25,7 +25,8 @@
 
 -----------------------------------------------------------------------------
 
-import qualified Data.Packed.Vector as V
+--import qualified Data.Packed.Vector as V
+import qualified Data.Vector.Storable as V
 
 import qualified Numeric.GSL.Histogram as H
 import qualified Numeric.GSL.Histogram2D as H2
@@ -46,7 +47,7 @@
     probability :: b -> V.Vector a -> V.Vector Double
 
 instance Storable b => PDF (PDFFunction b) b where
-    probability (P_Func f) = V.mapVector f
+    probability (P_Func f) = V.map f
 
 instance PDF H.Histogram Double where
     probability = H.prob
diff --git a/lib/Numeric/Statistics/Surrogate.hs b/lib/Numeric/Statistics/Surrogate.hs
--- a/lib/Numeric/Statistics/Surrogate.hs
+++ b/lib/Numeric/Statistics/Surrogate.hs
@@ -1,7 +1,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.Surrogate
--- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010
+-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010, 2014
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -19,8 +19,8 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Vector
---import Data.Packed.Matrix
+--import qualified Data.Packed.Vector as V
+import qualified Data.Vector.Storable as V
 
 import qualified Data.Array.IArray as I 
 
@@ -33,13 +33,13 @@
 -- | perform an analysis using surrogate data
 surrogate :: Int                          -- ^ random seed
           -> Int                          -- ^ number of repetitions
-          -> (I.Array Int (Vector Double) -> a) -- ^ the evaluation function
-          -> I.Array Int (Vector Double)  -- ^ the data
+          -> (I.Array Int (V.Vector Double) -> a) -- ^ the evaluation function
+          -> I.Array Int (V.Vector Double)  -- ^ the data
           -> I.Array Int a                -- ^ the results, with the evaluated real data in position 1 and the rest of the array containing the evaluated surrogate data
                                                                                                               
 surrogate r n f d = I.listArray (1,n+1) $ (f d) : (surrogate' (mkStdGen r) n f d)
 
-surrogate' :: StdGen -> Int-> (I.Array Int (Vector Double) -> a) -> I.Array Int (Vector Double) -> [a]
+surrogate' :: StdGen -> Int-> (I.Array Int (V.Vector Double) -> a) -> I.Array Int (V.Vector Double) -> [a]
 surrogate' _ 0 _ _ = []
 surrogate' g n f d = let (g',g'') = split g
                          d' = permute_data g' d
@@ -50,10 +50,10 @@
 randomList g n = let (r,g') = random g
                  in r : (randomList g' (n-1))
 
-permute_data :: StdGen -> I.Array Int (Vector Double) -> I.Array Int (Vector Double)
+permute_data :: StdGen -> I.Array Int (V.Vector Double) -> I.Array Int (V.Vector Double)
 permute_data g d = let s = I.rangeSize $ I.bounds d
                        rs = randomList g s
                        ds = zip rs $ I.elems d
-                   in I.listArray (I.bounds d) $ map (\(r,v) -> permute (random_permute r (dim v)) v) ds
+                   in I.listArray (I.bounds d) $ map (\(r,v) -> permute (random_permute r (V.length v)) v) ds
  
 -----------------------------------------------------------------------------