diff --git a/Statistics/Autocorrelation.hs b/Statistics/Autocorrelation.hs
--- a/Statistics/Autocorrelation.hs
+++ b/Statistics/Autocorrelation.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module    : Statistics.Autocorrelation
 -- Copyright : (c) 2009 Bryan O'Sullivan
@@ -16,31 +17,31 @@
     , autocorrelation
     ) where
 
-import Statistics.Sample (Sample, mean)
-import qualified Data.Vector.Unboxed as U
+import Statistics.Sample (mean)
+import qualified Data.Vector.Generic as G
 
 -- | Compute the autocovariance of a sample, i.e. the covariance of
 -- the sample against a shifted version of itself.
-autocovariance :: Sample -> U.Vector Double
-autocovariance a = U.map f . U.enumFromTo 0 $ l-2
+autocovariance :: (G.Vector v Double, G.Vector v Int) => v Double -> v Double
+autocovariance a = G.map f . G.enumFromTo 0 $ l-2
   where
-    f k = U.sum (U.zipWith (*) (U.take (l-k) c) (U.slice k (l-k) c))
+    f k = G.sum (G.zipWith (*) (G.take (l-k) c) (G.slice k (l-k) c))
           / fromIntegral l
-    c   = U.map (subtract (mean a)) a
-    l   = U.length a
+    c   = G.map (subtract (mean a)) a
+    l   = G.length a
 
 -- | Compute the autocorrelation function of a sample, and the upper
 -- and lower bounds of confidence intervals for each element.
 --
 -- /Note/: The calculation of the 95% confidence interval assumes a
 -- stationary Gaussian process.
-autocorrelation :: Sample -> (U.Vector Double, U.Vector Double, U.Vector Double)
+autocorrelation :: (G.Vector v Double, G.Vector v Int) => v Double -> (v Double, v Double, v Double)
 autocorrelation a = (r, ci (-), ci (+))
   where
-    r           = U.map (/ U.head c) c
+    r           = G.map (/ G.head c) c
       where c   = autocovariance a
-    dllse       = U.map f . U.scanl1 (+) . U.map square $ r
+    dllse       = G.map f . G.scanl1 (+) . G.map square $ r
       where f v = 1.96 * sqrt ((v * 2 + 1) / l)
-    l           = fromIntegral (U.length a)
-    ci f        = U.cons 1 . U.tail . U.map (f (-1/l)) $ dllse
+    l           = fromIntegral (G.length a)
+    ci f        = G.cons 1 . G.tail . G.map (f (-1/l)) $ dllse
     square x    = x * x
diff --git a/Statistics/Function.hs b/Statistics/Function.hs
--- a/Statistics/Function.hs
+++ b/Statistics/Function.hs
@@ -1,4 +1,5 @@
-{-# LANGUAGE Rank2Types, TypeOperators #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module    : Statistics.Function
 -- Copyright : (c) 2009, 2010 Bryan O'Sullivan
@@ -26,38 +27,38 @@
 import Data.Vector.Algorithms.Combinators (apply)
 import Data.Vector.Generic (unsafeFreeze)
 import qualified Data.Vector.Algorithms.Intro as I
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Unboxed.Mutable  as MU
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as M
 
 -- | Sort a vector.
-sort :: (U.Unbox e, Ord e) => U.Vector e -> U.Vector e
+sort :: (Ord e, G.Vector v e) => v e -> v e
 sort = apply I.sort
 {-# INLINE sort #-}
 
 -- | Partially sort a vector, such that the least /k/ elements will be
 -- at the front.
-partialSort :: (U.Unbox e, Ord e) =>
-               Int              -- ^ The number /k/ of least elements.
-            -> U.Vector e
-            -> U.Vector e
+partialSort :: (G.Vector v e, Ord e) =>
+               Int -- ^ The number /k/ of least elements.
+            -> v e
+            -> v e
 partialSort k = apply (\a -> I.partialSort a k)
 {-# INLINE partialSort #-}
 
 -- | Return the indices of a vector.
-indices :: (U.Unbox a) => U.Vector a -> U.Vector Int
-indices a = U.enumFromTo 0 (U.length a - 1)
+indices :: (G.Vector v a, G.Vector v Int) => v a -> v Int
+indices a = G.enumFromTo 0 (G.length a - 1)
 {-# INLINE indices #-}
 
 -- | Zip a vector with its indices.
-indexed :: U.Unbox e => U.Vector e -> U.Vector (Int,e)
-indexed a = U.zip (indices a) a
+indexed :: (G.Vector v e, G.Vector v Int, G.Vector v (Int,e)) => v e -> v (Int,e)
+indexed a = G.zip (indices a) a
 {-# INLINE indexed #-}
 
 data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
 -- | Compute the minimum and maximum of a vector in one pass.
-minMax :: U.Vector Double -> (Double , Double)
-minMax = fini . U.foldl go (MM (1/0) (-1/0))
+minMax :: (G.Vector v Double) => v Double -> (Double, Double)
+minMax = fini . G.foldl' go (MM (1/0) (-1/0))
   where
     go (MM lo hi) k = MM (min lo k) (max hi k)
     fini (MM lo hi) = (lo, hi)
@@ -65,12 +66,12 @@
 
 -- | Create a vector, using the given action to populate each
 -- element.
-create :: (PrimMonad m, U.Unbox e) => Int -> (Int -> m e) -> m (U.Vector e)
+create :: (PrimMonad m, G.Vector v e) => Int -> (Int -> m e) -> m (v e)
 create size itemAt = assert (size >= 0) $
-    MU.new size >>= loop 0
+    M.new size >>= loop 0
   where
     loop k arr | k >= size = unsafeFreeze arr
                | otherwise = do r <- itemAt k
-                                MU.write arr k r
+                                M.write arr k r
                                 loop (k+1) arr
 {-# INLINE create #-}
diff --git a/Statistics/KernelDensity.hs b/Statistics/KernelDensity.hs
--- a/Statistics/KernelDensity.hs
+++ b/Statistics/KernelDensity.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module    : Statistics.KernelDensity
 -- Copyright : (c) 2009 Bryan O'Sullivan
@@ -40,8 +41,8 @@
 import Statistics.Constants (m_1_sqrt_2, m_2_sqrt_pi)
 import Statistics.Function (minMax)
 import Statistics.Sample (stdDev)
-import Statistics.Types (Sample)
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic as G
 
 -- | Points from the range of a 'Sample'.
 newtype Points = Points {
@@ -61,10 +62,11 @@
 
 -- | Compute the optimal bandwidth from the observed data for the given
 -- kernel.
-bandwidth :: (Double -> Bandwidth)
-          -> Sample
+bandwidth :: G.Vector v Double =>
+             (Double -> Bandwidth)
+          -> v Double
           -> Bandwidth
-bandwidth kern values = stdDev values * kern (fromIntegral $ U.length values)
+bandwidth kern values = stdDev values * kern (fromIntegral $ G.length values)
 
 -- | Choose a uniform range of points at which to estimate a sample's
 -- probability density function.
@@ -74,9 +76,10 @@
 --
 -- If this function is passed an empty vector, it returns values of
 -- positive and negative infinity.
-choosePoints :: Int             -- ^ Number of points to select, /n/
+choosePoints :: G.Vector v Double =>
+                Int             -- ^ Number of points to select, /n/
              -> Double          -- ^ Sample bandwidth, /h/
-             -> Sample          -- ^ Input data
+             -> v Double        -- ^ Input data
              -> Points
 choosePoints n h sample = Points . U.map f $ U.enumFromTo 0 n'
   where lo     = a - h
@@ -116,27 +119,29 @@
 
 -- | Kernel density estimator, providing a non-parametric way of
 -- estimating the PDF of a random variable.
-estimatePDF :: Kernel           -- ^ Kernel function
+estimatePDF :: G.Vector v Double =>
+               Kernel           -- ^ Kernel function
             -> Bandwidth        -- ^ Bandwidth, /h/
-            -> Sample           -- ^ Sample data
+            -> v Double         -- ^ Sample data
             -> Points           -- ^ Points at which to estimate
             -> U.Vector Double
 estimatePDF kernel h sample
     | n < 2     = errorShort "estimatePDF"
     | otherwise = U.map k . fromPoints
   where
-    k p = U.sum . U.map (kernel f h p) $ sample
+    k p = G.sum . G.map (kernel f h p) $ sample
     f   = 1 / (h * fromIntegral n)
-    n   = U.length sample
+    n   = G.length sample
 {-# INLINE estimatePDF #-}
 
 -- | A helper for creating a simple kernel density estimation function
 -- with automatically chosen bandwidth and estimation points.
-simplePDF :: (Double -> Double) -- ^ Bandwidth function
+simplePDF :: G.Vector v Double =>
+             (Double -> Double) -- ^ Bandwidth function
           -> Kernel             -- ^ Kernel function
           -> Double             -- ^ Bandwidth scaling factor (3 for a Gaussian kernel, 1 for all others)
           -> Int                -- ^ Number of points at which to estimate
-          -> Sample             -- ^ Sample data
+          -> v Double           -- ^ sample data
           -> (Points, U.Vector Double)
 simplePDF fbw fpdf k numPoints sample =
     (points, estimatePDF fpdf bw sample points)
@@ -147,16 +152,18 @@
 -- | Simple Epanechnikov kernel density estimator.  Returns the
 -- uniformly spaced points from the sample range at which the density
 -- function was estimated, and the estimates at those points.
-epanechnikovPDF :: Int          -- ^ Number of points at which to estimate
-                -> Sample
+epanechnikovPDF :: G.Vector v Double =>
+                   Int          -- ^ Number of points at which to estimate
+                -> v Double     -- ^ Data sample
                 -> (Points, U.Vector Double)
 epanechnikovPDF = simplePDF epanechnikovBW epanechnikovKernel 1
 
 -- | Simple Gaussian kernel density estimator.  Returns the uniformly
 -- spaced points from the sample range at which the density function
 -- was estimated, and the estimates at those points.
-gaussianPDF :: Int              -- ^ Number of points at which to estimate
-            -> Sample
+gaussianPDF :: G.Vector v Double =>
+               Int              -- ^ Number of points at which to estimate
+            -> v Double         -- ^ Data sample
             -> (Points, U.Vector Double)
 gaussianPDF = simplePDF gaussianBW gaussianKernel 3
 
diff --git a/Statistics/Math.hs b/Statistics/Math.hs
--- a/Statistics/Math.hs
+++ b/Statistics/Math.hs
@@ -1,4 +1,5 @@
-{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE BangPatterns     #-}
+{-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module    : Statistics.Math
 -- Copyright : (c) 2009 Bryan O'Sullivan
@@ -26,25 +27,28 @@
     -- $references
     ) where
 
-import Data.Vector.Unboxed ((!))
+import Data.Vector.Generic ((!))
 import Data.Word (Word64)
 import Statistics.Constants (m_sqrt_2_pi)
 import Statistics.Distribution (cumulative)
 import Statistics.Distribution.Normal (standard)
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic as G
 
 data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
 -- | Evaluate a series of Chebyshev polynomials. Uses Clenshaw's
 -- algorithm.
-chebyshev :: Double             -- ^ Parameter of each function.
-          -> U.Vector Double    -- ^ Coefficients of each polynomial
-                                --   term, in increasing order.
+chebyshev :: (G.Vector v Double) =>
+             Double      -- ^ Parameter of each function.
+          -> v Double    -- ^ Coefficients of each polynomial term, in increasing order.
           -> Double
-chebyshev x a = fini . U.foldl step (C 0 0) $ U.enumFromStepN (U.length a - 1) (-1) (U.length a - 1)
+chebyshev x a = fini . U.foldl' step (C 0 0) $ U.enumFromStepN (len - 1) (-1) (len - 1)
     where step (C b1 b2) k = C ((a ! k) + x2 * b1 - b2) b1
           fini (C b1 b2)   = (a ! 0) + x * b1 - b2
-          x2                 = x * 2
+          x2               = x * 2
+          len              = G.length a
+{-# INLINE chebyshev #-}
 
 -- | The binomial coefficient.
 --
@@ -52,7 +56,7 @@
 choose :: Int -> Int -> Double
 n `choose` k
     | k > n     = 0
-    | k < 30    = U.foldl go 1 . U.enumFromTo 1 $ k'
+    | k < 30    = U.foldl' go 1 . U.enumFromTo 1 $ k'
     | otherwise = exp $ lg (n+1) - lg (k+1) - lg (n-k+1)
     where go a i = a * (nk + j) / j
               where j = fromIntegral i :: Double
@@ -71,8 +75,8 @@
 factorial n
     | n < 0     = error "Statistics.Math.factorial: negative input"
     | n <= 1    = 0
-    | n <= 14   = fini . U.foldl goLong (F 1 1) $ ns
-    | otherwise = U.foldl goDouble 1 $ ns
+    | n <= 14   = fini . U.foldl' goLong (F 1 1) $ ns
+    | otherwise = U.foldl' goDouble 1 $ ns
     where goDouble t k = t * fromIntegral k
           goLong (F z x) _ = F (z * x') x'
               where x' = x + 1
@@ -204,7 +208,7 @@
 logGammaL :: Double -> Double
 logGammaL x
     | x <= 0    = 1/0
-    | otherwise = fini . U.foldl go (L 0 (x+7)) $ a
+    | otherwise = fini . U.foldl' go (L 0 (x+7)) $ a
     where fini (L l _) = log (l+a0) + log m_sqrt_2_pi - x65 + (x-0.5) * log x65
           go (L l t) k = L (l + k / t) (t-1)
           x65 = x + 6.5
diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs
--- a/Statistics/Quantile.hs
+++ b/Statistics/Quantile.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module    : Statistics.Quantile
 -- Copyright : (c) 2009 Bryan O'Sullivan
@@ -38,27 +38,27 @@
     ) where
 
 import Control.Exception (assert)
-import Data.Vector.Unboxed ((!))
+import Data.Vector.Generic ((!))
 import Statistics.Constants (m_epsilon)
 import Statistics.Function (partialSort)
-import Statistics.Types (Sample)
-import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic as G
 
 -- | O(/n/ log /n/). Estimate the /k/th /q/-quantile of a sample,
 -- using the weighted average method.
-weightedAvg :: Int              -- ^ /k/, the desired quantile.
-            -> Int              -- ^ /q/, the number of quantiles.
-            -> Sample           -- ^ /x/, the sample data.
+weightedAvg :: G.Vector v Double => 
+               Int        -- ^ /k/, the desired quantile.
+            -> Int        -- ^ /q/, the number of quantiles.
+            -> v Double   -- ^ /x/, the sample data.
             -> Double
 weightedAvg k q x =
     assert (q >= 2) .
     assert (k >= 0) .
     assert (k < q) .
-    assert (U.all (not . isNaN) x) $
+    assert (G.all (not . isNaN) x) $
     xj + g * (xj1 - xj)
   where
     j   = floor idx
-    idx = fromIntegral (U.length x - 1) * fromIntegral k / fromIntegral q
+    idx = fromIntegral (G.length x - 1) * fromIntegral k / fromIntegral q
     g   = idx - fromIntegral j
     xj  = sx ! j
     xj1 = sx ! (j+1)
@@ -72,16 +72,17 @@
 -- using the continuous sample method with the given parameters.  This
 -- is the method used by most statistical software, such as R,
 -- Mathematica, SPSS, and S.
-continuousBy :: ContParam       -- ^ Parameters /a/ and /b/.
-             -> Int             -- ^ /k/, the desired quantile.
-             -> Int             -- ^ /q/, the number of quantiles.
-             -> Sample          -- ^ /x/, the sample data.
+continuousBy :: G.Vector v Double =>
+                ContParam  -- ^ Parameters /a/ and /b/.
+             -> Int        -- ^ /k/, the desired quantile.
+             -> Int        -- ^ /q/, the number of quantiles.
+             -> v Double   -- ^ /x/, the sample data.
              -> Double
 continuousBy (ContParam a b) k q x =
     assert (q >= 2) .
     assert (k >= 0) .
     assert (k <= q) .
-    assert (U.all (not . isNaN) x) $
+    assert (G.all (not . isNaN) x) $
     (1-h) * item (j-1) + h * item j
   where
     j               = floor (t + eps)
@@ -91,7 +92,7 @@
       | otherwise   = r
       where r       = t - fromIntegral j
     eps             = m_epsilon * 4
-    n               = U.length x
+    n               = G.length x
     item            = (sx !) . bracket
     sx              = partialSort (bracket j + 1) x
     bracket m       = min (max m 0) (n - 1)
@@ -104,14 +105,15 @@
 -- For instance, the interquartile range (IQR) can be estimated as
 -- follows:
 --
--- > midspread medianUnbiased 4 (U.to [1,1,2,2,3])
+-- > midspread medianUnbiased 4 (U.fromList [1,1,2,2,3])
 -- > ==> 1.333333
-midspread :: ContParam       -- ^ Parameters /a/ and /b/.
-          -> Int             -- ^ /q/, the number of quantiles.
-          -> Sample          -- ^ /x/, the sample data.
+midspread :: G.Vector v Double =>
+             ContParam  -- ^ Parameters /a/ and /b/.
+          -> Int        -- ^ /q/, the number of quantiles.
+          -> v Double   -- ^ /x/, the sample data.
           -> Double
 midspread (ContParam a b) k x =
-    assert (U.all (not . isNaN) x) .
+    assert (G.all (not . isNaN) x) .
     assert (k > 0) $
     quantile (1-frac) - quantile frac
   where
@@ -122,7 +124,7 @@
         | otherwise   = r
         where r       = t i - fromIntegral (j i)
     eps               = m_epsilon * 4
-    n                 = U.length x
+    n                 = G.length x
     item              = (sx !) . bracket
     sx                = partialSort (bracket (j (1-frac)) + 1) x
     bracket m         = min (max m 0) (n - 1)
diff --git a/Statistics/Resampling/Bootstrap.hs b/Statistics/Resampling/Bootstrap.hs
--- a/Statistics/Resampling/Bootstrap.hs
+++ b/Statistics/Resampling/Bootstrap.hs
@@ -83,7 +83,7 @@
         ni    = U.length resample
         n     = fromIntegral ni
         accel = sumCubes / (6 * (sumSquares ** 1.5))
-          where (sumSquares :< sumCubes) = U.foldl f (0 :< 0) jack
+          where (sumSquares :< sumCubes) = U.foldl' f (0 :< 0) jack
                 f (s :< c) j = s + d2 :< c + d2 * d
                     where d  = jackMean - j
                           d2 = d * d
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module    : Statistics.Sample
 -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan
@@ -15,6 +15,7 @@
     (
     -- * Types
       Sample
+    , WeightedSample
     -- * Descriptive functions
     , range
 
@@ -52,18 +53,20 @@
 
 import Statistics.Function (minMax)
 import Statistics.Types (Sample,WeightedSample)
-import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic as G
 
+-- Operator ^ will be overriden
+import Prelude hiding ((^))
 
-range :: Sample -> Double
+range :: (G.Vector v Double) => v Double -> Double
 range s = hi - lo
     where (lo , hi) = minMax s
 {-# INLINE range #-}
 
 -- | Arithmetic mean.  This uses Welford's algorithm to provide
 -- numerical stability, using a single pass over the sample data.
-mean :: Sample -> Double
-mean = fini . U.foldl go (T 0 0)
+mean :: (G.Vector v Double) => v Double -> Double
+mean = fini . G.foldl' go (T 0 0)
   where
     fini (T a _) = a
     go (T m n) x = T m' n'
@@ -73,8 +76,8 @@
 
 -- | Arithmetic mean for weighted sample. It uses algorithm analogous
 --   to one in 'mean'
-meanWeighted :: WeightedSample -> Double
-meanWeighted = fini . U.foldl go (V 0 0)
+meanWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double
+meanWeighted = fini . G.foldl' go (V 0 0)
     where
       fini (V a _) = a
       go (V m w) (x,xw) = V m' w'
@@ -85,16 +88,16 @@
 
 -- | Harmonic mean.  This algorithm performs a single pass over the
 -- sample.
-harmonicMean :: Sample -> Double
-harmonicMean = fini . U.foldl go (T 0 0)
+harmonicMean :: (G.Vector v Double) => v Double -> Double
+harmonicMean = fini . G.foldl' go (T 0 0)
   where
     fini (T b a) = fromIntegral a / b
     go (T x y) n = T (x + (1/n)) (y+1)
 {-# INLINE harmonicMean #-}
 
 -- | Geometric mean of a sample containing no negative values.
-geometricMean :: Sample -> Double
-geometricMean = fini . U.foldl go (T 1 0)
+geometricMean :: (G.Vector v Double) => v Double -> Double
+geometricMean = fini . G.foldl' go (T 1 0)
   where
     fini (T p n) = p ** (1 / fromIntegral n)
     go (T p n) a = T (p * a) (n + 1)
@@ -108,12 +111,12 @@
 --
 -- For samples containing many values very close to the mean, this
 -- function is subject to inaccuracy due to catastrophic cancellation.
-centralMoment :: Int -> Sample -> Double
+centralMoment :: (G.Vector v Double) => Int -> v Double -> Double
 centralMoment a xs
     | a < 0  = error "Statistics.Sample.centralMoment: negative input"
     | a == 0 = 1
     | a == 1 = 0
-    | otherwise = U.sum (U.map go xs) / fromIntegral (U.length xs)
+    | otherwise = G.sum (G.map go xs) / fromIntegral (G.length xs)
   where
     go x = (x-m) ^ a
     m    = mean xs
@@ -126,15 +129,15 @@
 --
 -- For samples containing many values very close to the mean, this
 -- function is subject to inaccuracy due to catastrophic cancellation.
-centralMoments :: Int -> Int -> Sample -> (Double, Double)
+centralMoments :: (G.Vector v Double) => Int -> Int -> v Double -> (Double, Double)
 centralMoments a b xs
     | a < 2 || b < 2 = (centralMoment a xs , centralMoment b xs)
-    | otherwise      = fini . U.foldl go (V 0 0) $ xs
+    | otherwise      = fini . G.foldl' go (V 0 0) $ xs
   where go (V i j) x = V (i + d^a) (j + d^b)
             where d  = x - m
         fini (V i j) = (i / n , j / n)
         m            = mean xs
-        n            = fromIntegral (U.length xs)
+        n            = fromIntegral (G.length xs)
 {-# INLINE centralMoments #-}
 
 -- | Compute the skewness of a sample. This is a measure of the
@@ -159,7 +162,7 @@
 --
 -- For samples containing many values very close to the mean, this
 -- function is subject to inaccuracy due to catastrophic cancellation.
-skewness :: Sample -> Double
+skewness :: (G.Vector v Double) => v Double -> Double
 skewness xs = c3 * c2 ** (-1.5)
     where (c3 , c2) = centralMoments 3 2 xs
 {-# INLINE skewness #-}
@@ -177,7 +180,7 @@
 --
 -- For samples containing many values very close to the mean, this
 -- function is subject to inaccuracy due to catastrophic cancellation.
-kurtosis :: Sample -> Double
+kurtosis :: (G.Vector v Double) => v Double -> Double
 kurtosis xs = c4 / (c2 * c2) - 3
     where (c4 , c2) = centralMoments 4 2 xs
 {-# INLINE kurtosis #-}
@@ -201,48 +204,50 @@
 sqr :: Double -> Double
 sqr x = x * x
 
-robustSumVar :: Sample -> Double
-robustSumVar samp = U.sum . U.map (sqr . subtract m) $ samp
+robustSumVar :: (G.Vector v Double) => v Double -> Double
+robustSumVar samp = G.sum . G.map (sqr . subtract m) $ samp
     where
       m = mean samp
+{-# INLINE robustSumVar #-}
 
 -- | Maximum likelihood estimate of a sample's variance.  Also known
 -- as the population variance, where the denominator is /n/.
-variance :: Sample -> Double
+variance :: (G.Vector v Double) => v Double -> Double
 variance samp
     | n > 1     = robustSumVar samp / fromIntegral n
     | otherwise = 0
     where
-      n = U.length samp
+      n = G.length samp
 {-# INLINE variance #-}
 
 -- | Unbiased estimate of a sample's variance.  Also known as the
 -- sample variance, where the denominator is /n/-1.
-varianceUnbiased :: Sample -> Double
+varianceUnbiased :: (G.Vector v Double) => v Double -> Double
 varianceUnbiased samp
     | n > 1     = robustSumVar samp / fromIntegral (n-1)
     | otherwise = 0
     where
-      n = U.length samp
+      n = G.length samp
 {-# INLINE varianceUnbiased #-}
 
 -- | Standard deviation.  This is simply the square root of the
 -- unbiased estimate of the variance.
-stdDev :: Sample -> Double
+stdDev :: (G.Vector v Double) => v Double -> Double
 stdDev = sqrt . varianceUnbiased
-
+{-# INLINE stdDev #-}
 
-robustSumVarWeighted :: WeightedSample -> V
-robustSumVarWeighted samp = U.foldl go (V 0 0) samp
+robustSumVarWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> V
+robustSumVarWeighted samp = G.foldl' go (V 0 0) samp
     where
       go (V s w) (x,xw) = V (s + xw*d*d) (w + xw)
           where d = x - m
       m = meanWeighted samp
+{-# INLINE robustSumVarWeighted #-}
 
 -- | Weighted variance. This is biased estimation.
-varianceWeighted :: WeightedSample -> Double
+varianceWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double
 varianceWeighted samp
-    | U.length samp > 1 = fini $ robustSumVarWeighted samp
+    | G.length samp > 1 = fini $ robustSumVarWeighted samp
     | otherwise         = 0
     where
       fini (V s w) = s / w
@@ -259,8 +264,8 @@
 -- mean, Knuth's algorithm gives inaccurate results due to
 -- catastrophic cancellation.
 
-fastVar :: Sample -> T1
-fastVar = U.foldl go (T1 0 0 0)
+fastVar :: (G.Vector v Double) => v Double -> T1
+fastVar = G.foldl' go (T1 0 0 0)
   where
     go (T1 n m s) x = T1 n' m' s'
       where n' = n + 1
@@ -269,7 +274,7 @@
             d  = x - m
 
 -- | Maximum likelihood estimate of a sample's variance.
-fastVariance :: Sample -> Double
+fastVariance :: (G.Vector v Double) => v Double -> Double
 fastVariance = fini . fastVar
   where fini (T1 n _m s)
           | n > 1     = s / fromIntegral n
@@ -277,7 +282,7 @@
 {-# INLINE fastVariance #-}
 
 -- | Unbiased estimate of a sample's variance.
-fastVarianceUnbiased :: Sample -> Double
+fastVarianceUnbiased :: (G.Vector v Double) => v Double -> Double
 fastVarianceUnbiased = fini . fastVar
   where fini (T1 n _m s)
           | n > 1     = s / fromIntegral (n - 1)
@@ -286,12 +291,18 @@
 
 -- | Standard deviation.  This is simply the square root of the
 -- maximum likelihood estimate of the variance.
-fastStdDev :: Sample -> Double
+fastStdDev :: (G.Vector v Double) => v Double -> Double
 fastStdDev = sqrt . fastVariance
 {-# INLINE fastStdDev #-}
 
 ------------------------------------------------------------------------
 -- Helper code. Monomorphic unpacked accumulators.
+
+-- (^) operator from Prelude is just slow.
+(^) :: Double -> Int -> Double
+x ^ 1 = x
+x ^ n = x * (x ^ (n-1))
+{-# INLINE (^) #-}
 
 -- don't support polymorphism, as we can't get unboxed returns if we use it.
 data T = T {-# UNPACK #-}!Double {-# UNPACK #-}!Int
diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs
--- a/Statistics/Sample/Powers.hs
+++ b/Statistics/Sample/Powers.hs
@@ -1,4 +1,5 @@
-{-# LANGUAGE BangPatterns, TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE BangPatterns     #-}
 -- |
 -- Module    : Statistics.Sample.Powers
 -- Copyright : (c) 2009, 2010 Bryan O'Sullivan
@@ -19,8 +20,7 @@
 module Statistics.Sample.Powers
     (
     -- * Types
-      Sample
-    , Powers
+      Powers
 
     -- * Constructor
     , powers
@@ -53,9 +53,9 @@
 import Statistics.Function (indexed)
 import Statistics.Internal (inlinePerformIO)
 import Statistics.Math (choose)
-import Statistics.Types (Sample)
 import System.IO.Unsafe (unsafePerformIO)
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed.Mutable as MU
 
 newtype Powers = Powers (U.Vector Double)
@@ -76,12 +76,13 @@
 -- * For 'kurtosis', at least 4 simple powers are required.
 --
 -- This function is subject to stream fusion.
-powers :: Int                   -- ^ /n/, the number of powers, where /n/ >= 2.
-       -> Sample
+powers :: G.Vector v Double =>
+          Int                   -- ^ /n/, the number of powers, where /n/ >= 2.
+       -> v Double
        -> Powers
 powers k
     | k < 2     = error "Statistics.Sample.powers: too few powers"
-    | otherwise = fini . U.foldl go (unsafePerformIO $ create)
+    | otherwise = fini . G.foldl' go (unsafePerformIO $ MU.newWith l 0)
   where
     go ms x = inlinePerformIO $ loop 0 1
         where loop !i !xk | i == l = return ms
@@ -89,18 +90,15 @@
                 MU.read ms i >>= MU.write ms i . (+ xk)
                 loop (i+1) (xk*x)
     fini = Powers . unsafePerformIO . unsafeFreeze
-    create = MU.new l >>= fill 0
-        where fill !i ms | i == l    = return ms
-                         | otherwise = MU.write ms i 0 >> fill (i+1) ms
-    l = k + 1
+    l    = k + 1
 {-# INLINE powers #-}
 
--- | The order (number) of simple powers collected from a 'Sample'.
+-- | The order (number) of simple powers collected from a 'sample'.
 order :: Powers -> Int
 order (Powers pa) = U.length pa - 1
 {-# INLINE order #-}
 
--- | Compute the /k/th central moment of a 'Sample'.  The central
+-- | Compute the /k/th central moment of a sample.  The central
 -- moment is also known as the moment about the mean.
 centralMoment :: Int -> Powers -> Double
 centralMoment k p@(Powers pa)
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.5.1.2
+version:        0.6.0.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful