diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2009, Bryan O'Sullivan
+Copyright (c) 2009, 2010 Bryan O'Sullivan
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
diff --git a/Statistics/Autocorrelation.hs b/Statistics/Autocorrelation.hs
--- a/Statistics/Autocorrelation.hs
+++ b/Statistics/Autocorrelation.hs
@@ -16,8 +16,8 @@
     , autocorrelation
     ) where
 
-import qualified Data.Vector.Unboxed as U
 import Statistics.Sample (Sample, mean)
+import qualified Data.Vector.Unboxed as U
 
 -- | Compute the autocovariance of a sample, i.e. the covariance of
 -- the sample against a shifted version of itself.
diff --git a/Statistics/Function.hs b/Statistics/Function.hs
--- a/Statistics/Function.hs
+++ b/Statistics/Function.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE Rank2Types, TypeOperators #-}
 -- |
 -- Module    : Statistics.Function
--- Copyright : (c) 2009 Bryan O'Sullivan
+-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
 -- License   : BSD3
 --
 -- Maintainer  : bos@serpentine.com
@@ -15,26 +15,26 @@
       minMax
     , sort
     , partialSort
+    , indexed
     , indices
-    -- * Array setup
-    , createU
-    , createIO
+    -- * Vector setup
+    , create
     ) where
 
 import Control.Exception (assert)
-import Control.Monad.ST (ST, unsafeIOToST, unsafeSTToIO)
+import Control.Monad.Primitive (PrimMonad)
 import Data.Vector.Algorithms.Combinators (apply)
-import qualified Data.Vector.Unboxed as U
 import Data.Vector.Generic (unsafeFreeze)
-import qualified Data.Vector.Unboxed.Mutable  as MU
 import qualified Data.Vector.Algorithms.Intro as I
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable  as MU
 
--- | Sort an array.
+-- | Sort a vector.
 sort :: (U.Unbox e, Ord e) => U.Vector e -> U.Vector e
 sort = apply I.sort
 {-# INLINE sort #-}
 
--- | Partially sort an array, such that the least /k/ elements will be
+-- | 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.
@@ -43,36 +43,34 @@
 partialSort k = apply (\a -> I.partialSort a k)
 {-# INLINE partialSort #-}
 
--- | Return the indices of an array.
+-- | 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)
 {-# 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
+{-# INLINE indexed #-}
+
 data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
--- | Compute the minimum and maximum of an array in one pass.
+-- | 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))
   where
     go (MM lo hi) k = MM (min lo k) (max hi k)
-    fini (MM lo hi) = (lo , hi)
+    fini (MM lo hi) = (lo, hi)
 {-# INLINE minMax #-}
 
--- | Create an array, using the given 'ST' action to populate each
+-- | Create a vector, using the given action to populate each
 -- element.
-createU :: (U.Unbox e) => forall s. Int -> (Int -> ST s e) -> ST s (U.Vector e)
-createU size itemAt = assert (size >= 0) $
+create :: (PrimMonad m, U.Unbox e) => Int -> (Int -> m e) -> m (U.Vector e)
+create size itemAt = assert (size >= 0) $
     MU.new size >>= loop 0
   where
     loop k arr | k >= size = unsafeFreeze arr
                | otherwise = do r <- itemAt k
                                 MU.write arr k r
                                 loop (k+1) arr
-{-# INLINE createU #-}
-
--- | Create an array, using the given 'IO' action to populate each
--- element.
-createIO :: (U.Unbox e) => Int -> (Int -> IO e) -> IO (U.Vector e)
-createIO size itemAt =
-    unsafeSTToIO $ createU size (unsafeIOToST . itemAt)
-{-# INLINE createIO #-}
+{-# INLINE create #-}
diff --git a/Statistics/KernelDensity.hs b/Statistics/KernelDensity.hs
--- a/Statistics/KernelDensity.hs
+++ b/Statistics/KernelDensity.hs
@@ -37,11 +37,11 @@
     , simplePDF
     ) where
 
-import qualified Data.Vector.Unboxed as U
+import Statistics.Constants (m_1_sqrt_2, m_2_sqrt_pi)
 import Statistics.Function (minMax)
 import Statistics.Sample (stdDev)
-import Statistics.Constants (m_1_sqrt_2, m_2_sqrt_pi)
 import Statistics.Types (Sample)
+import qualified Data.Vector.Unboxed as U
 
 -- | Points from the range of a 'Sample'.
 newtype Points = Points {
diff --git a/Statistics/Math.hs b/Statistics/Math.hs
--- a/Statistics/Math.hs
+++ b/Statistics/Math.hs
@@ -26,12 +26,12 @@
     -- $references
     ) where
 
-import qualified Data.Vector.Unboxed as U
 import Data.Vector.Unboxed ((!))
 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
 
 data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs
--- a/Statistics/Quantile.hs
+++ b/Statistics/Quantile.hs
@@ -38,11 +38,11 @@
     ) where
 
 import Control.Exception (assert)
-import qualified Data.Vector.Unboxed as U
 import Data.Vector.Unboxed ((!))
 import Statistics.Constants (m_epsilon)
 import Statistics.Function (partialSort)
 import Statistics.Types (Sample)
+import qualified Data.Vector.Unboxed as U
 
 -- | O(/n/ log /n/). Estimate the /k/th /q/-quantile of a sample,
 -- using the weighted average method.
diff --git a/Statistics/RandomVariate.hs b/Statistics/RandomVariate.hs
deleted file mode 100644
--- a/Statistics/RandomVariate.hs
+++ /dev/null
@@ -1,6 +0,0 @@
-module Statistics.RandomVariate
-    (
-      module System.Random.MWC
-    ) where
-
-import System.Random.MWC
diff --git a/Statistics/Resampling.hs b/Statistics/Resampling.hs
--- a/Statistics/Resampling.hs
+++ b/Statistics/Resampling.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module    : Statistics.Resampling
--- Copyright : (c) 2009 Bryan O'Sullivan
+-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
 -- License   : BSD3
 --
 -- Maintainer  : bos@serpentine.com
@@ -16,16 +16,16 @@
     , resample
     ) where
 
-import Control.Monad (forM_)
-import Control.Monad.ST (ST)
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Unboxed.Mutable as MU
-import Data.Vector.Unboxed ((!))
-import Data.Vector.Generic (unsafeFreeze)
+import Control.Monad (forM_, liftM)
+import Control.Monad.Primitive (PrimMonad, PrimState)
 import Data.Vector.Algorithms.Intro (sort)
-import Statistics.Function (createU, indices)
-import System.Random.MWC (Gen, uniform)
+import Data.Vector.Generic (unsafeFreeze)
+import Data.Vector.Unboxed ((!))
+import Statistics.Function (create, indexed, indices)
 import Statistics.Types (Estimator, Sample)
+import System.Random.MWC (Gen, uniform)
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as MU
 
 -- | A resample drawn randomly, with replacement, from a set of data
 -- points.  Distinct from a normal array to make it harder for your
@@ -36,16 +36,16 @@
 
 -- | Resample a data set repeatedly, with replacement, computing each
 -- estimate over the resampled data.
-resample :: Gen s -> [Estimator] -> Int -> Sample -> ST s [Resample]
+resample :: (PrimMonad m) => Gen (PrimState m) -> [Estimator] -> Int -> Sample -> m [Resample]
 resample gen ests numResamples samples = do
   results <- mapM (const (MU.new numResamples)) $ ests
   loop 0 (zip ests results)
   mapM_ sort results
-  mapM (fmap Resample . unsafeFreeze) results
+  mapM (liftM Resample . unsafeFreeze) results
  where
   loop k ers | k >= numResamples = return ()
              | otherwise = do
-    re <- createU n $ \_ -> do
+    re <- create n $ \_ -> do
             r <- uniform gen
             return (samples ! (abs r `mod` n))
     forM_ ers $ \(est,arr) ->
@@ -59,10 +59,6 @@
 jackknife est sample = U.map f . indices $ sample
     where f i = est (dropAt i sample)
 {-# INLINE jackknife #-}
-
--- Reimplementation of indexed
-indexed :: U.Unbox e => U.Vector e -> U.Vector (Int,e)
-indexed a = U.zip (U.enumFromN 0 (U.length a)) a
 
 -- | Drop the /k/th element of a vector.
 dropAt :: U.Unbox e => Int -> U.Vector e -> U.Vector e
diff --git a/Statistics/Resampling/Bootstrap.hs b/Statistics/Resampling/Bootstrap.hs
--- a/Statistics/Resampling/Bootstrap.hs
+++ b/Statistics/Resampling/Bootstrap.hs
@@ -18,13 +18,13 @@
     ) where
 
 import Control.Exception (assert)
-import qualified Data.Vector.Unboxed as U
 import Data.Vector.Unboxed ((!))
-import Statistics.Distribution.Normal
 import Statistics.Distribution (cumulative, quantile)
+import Statistics.Distribution.Normal
 import Statistics.Resampling (Resample(..), jackknife)
 import Statistics.Sample (mean)
 import Statistics.Types (Estimator, Sample)
+import qualified Data.Vector.Unboxed as U
 
 -- | A point and interval estimate computed via an 'Estimator'.
 data Estimate = Estimate {
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -50,9 +50,9 @@
     -- $references
     ) where
 
-import qualified Data.Vector.Unboxed as U
 import Statistics.Function (minMax)
 import Statistics.Types (Sample,WeightedSample)
+import qualified Data.Vector.Unboxed as U
 
 
 range :: Sample -> Double
diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs
--- a/Statistics/Sample/Powers.hs
+++ b/Statistics/Sample/Powers.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE BangPatterns, TypeOperators #-}
 -- |
 -- Module    : Statistics.Sample.Powers
--- Copyright : (c) 2009 Bryan O'Sullivan
+-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
 -- License   : BSD3
 --
 -- Maintainer  : bos@serpentine.com
@@ -47,16 +47,16 @@
     -- $references
     ) where
 
-import Control.Monad.ST (unsafeSTToIO)
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Unboxed.Mutable as MU
 import Data.Vector.Generic (unsafeFreeze)
 import Data.Vector.Unboxed ((!))
 import Prelude hiding (sum)
+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.Unboxed.Mutable as MU
 
 newtype Powers = Powers (U.Vector Double)
     deriving (Eq, Show)
@@ -81,14 +81,14 @@
        -> Powers
 powers k
     | k < 2     = error "Statistics.Sample.powers: too few powers"
-    | otherwise = fini . U.foldl go (unsafePerformIO . unsafeSTToIO $ create)
+    | otherwise = fini . U.foldl go (unsafePerformIO $ create)
   where
-    go ms x = inlinePerformIO . unsafeSTToIO $ loop 0 1
+    go ms x = inlinePerformIO $ loop 0 1
         where loop !i !xk | i == l = return ms
                           | otherwise = do
                 MU.read ms i >>= MU.write ms i . (+ xk)
                 loop (i+1) (xk*x)
-    fini = Powers . unsafePerformIO . unsafeSTToIO . unsafeFreeze
+    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
@@ -99,10 +99,6 @@
 order :: Powers -> Int
 order (Powers pa) = U.length pa - 1
 {-# INLINE order #-}
-
--- Reimplementation of indexed
-indexed :: U.Unbox e => U.Vector e -> U.Vector (Int,e)
-indexed a = U.zip (U.enumFromN 0 (U.length a)) a
 
 -- | Compute the /k/th central moment of a 'Sample'.  The central
 -- moment is also known as the moment about the mean.
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.5.0.0
+version:        0.5.1.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -22,7 +22,7 @@
 homepage:       http://darcs.serpentine.com/statistics
 author:         Bryan O'Sullivan <bos@serpentine.com>
 maintainer:     Bryan O'Sullivan <bos@serpentine.com>
-copyright:      2009 Bryan O'Sullivan
+copyright:      2009, 2010 Bryan O'Sullivan
 category:       Math, Statistics
 build-type:     Simple
 cabal-version:  >= 1.2
@@ -44,7 +44,6 @@
     Statistics.KernelDensity
     Statistics.Math
     Statistics.Quantile
-    Statistics.RandomVariate
     Statistics.Resampling
     Statistics.Resampling.Bootstrap
     Statistics.Sample
@@ -56,6 +55,7 @@
     base < 5,
     erf,
     mwc-random >= 0.5.0.0,
+    primitive,
     time,
     vector >= 0.5,
     vector-algorithms >= 0.3