diff --git a/Statistics/Regression.hs b/Statistics/Regression.hs
--- a/Statistics/Regression.hs
+++ b/Statistics/Regression.hs
@@ -10,15 +10,25 @@
       olsRegress
     , ols
     , rSquare
+    , bootstrapRegress
     ) where
 
 import Control.Applicative ((<$>))
+import Control.Concurrent (forkIO)
+import Control.Concurrent.Chan (newChan, readChan, writeChan)
+import Control.DeepSeq (rnf)
+import Control.Monad (forM_, replicateM)
+import GHC.Conc (getNumCapabilities)
 import Prelude hiding (pred, sum)
 import Statistics.Function as F
 import Statistics.Matrix hiding (map)
 import Statistics.Matrix.Algorithms (qr)
+import Statistics.Resampling (splitGen)
+import Statistics.Resampling.Bootstrap (Estimate(..))
 import Statistics.Sample (mean)
 import Statistics.Sample.Internal (sum)
+import System.Random.MWC (GenIO, uniformR)
+import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as M
@@ -97,3 +107,50 @@
     r   = sum $ flip U.imap resp $ \i x -> square (x - p i)
     t   = sum $ flip U.map resp $ \x -> square (x - mean resp)
     p i = sum . flip U.imap coeff $ \j -> (* unsafeIndex pred i j)
+
+-- | Bootstrap a regression function.  Returns both the results of the
+-- regression and the requested confidence interval values.
+bootstrapRegress :: GenIO
+                 -> Int         -- ^ Number of resamples to compute.
+                 -> Double      -- ^ Confidence interval.
+                 -> ([Vector] -> Vector -> (Vector, Double))
+                 -- ^ Regression function.
+                 -> [Vector]    -- ^ Predictor vectors.
+                 -> Vector      -- ^ Responder vector.
+                 -> IO (V.Vector Estimate, Estimate)
+bootstrapRegress gen0 numResamples ci rgrss preds0 resp0
+  | numResamples < 1   = error $ "bootstrapRegress: number of resamples " ++
+                                 "must be positive"
+  | ci <= 0 || ci >= 1 = error $ "bootstrapRegress: confidence interval " ++
+                                 "must lie between 0 and 1"
+  | otherwise = do
+  caps <- getNumCapabilities
+  gens <- splitGen caps gen0
+  done <- newChan
+  forM_ (zip gens (balance caps numResamples)) $ \(gen,count) -> do
+    forkIO $ do
+      v <- V.replicateM count $ do
+           let n = U.length resp0
+           ixs <- U.replicateM n $ uniformR (0,n-1) gen
+           let resp  = U.backpermute resp0 ixs
+               preds = map (flip U.backpermute ixs) preds0
+           return $ rgrss preds resp
+      rnf v `seq` writeChan done v
+  (coeffsv, r2v) <- (G.unzip . V.concat) <$> replicateM caps (readChan done)
+  let coeffs  = flip G.imap (G.convert coeffss) $ \i x ->
+                est x . U.generate numResamples $ \k -> ((coeffsv G.! k) G.! i)
+      r2      = est r2s (G.convert r2v)
+      (coeffss, r2s) = rgrss preds0 resp0
+      est s v = Estimate s (w G.! lo) (w G.! hi) ci
+        where w  = F.sort v
+              lo = round c
+              hi = truncate (n - c)
+              n  = fromIntegral numResamples
+              c  = n * ((1 - ci) / 2)
+  return (coeffs, r2)
+
+-- | Balance units of work across workers.
+balance :: Int -> Int -> [Int]
+balance numSlices numItems = zipWith (+) (replicate numSlices q)
+                                         (replicate r 1 ++ repeat 0)
+ where (q,r) = numItems `quotRem` numSlices
diff --git a/Statistics/Resampling.hs b/Statistics/Resampling.hs
--- a/Statistics/Resampling.hs
+++ b/Statistics/Resampling.hs
@@ -21,12 +21,12 @@
     , jackknifeStdDev
     , resample
     , estimate
+    , splitGen
     ) where
 
 import Data.Aeson (FromJSON, ToJSON)
 import Control.Concurrent (forkIO, newChan, readChan, writeChan)
-import Control.Monad (forM_, liftM, replicateM_)
-import Control.Monad.Primitive (PrimState)
+import Control.Monad (forM_, liftM, replicateM, replicateM_)
 import Data.Binary (Binary(..))
 import Data.Data (Data, Typeable)
 import Data.Vector.Algorithms.Intro (sort)
@@ -39,7 +39,7 @@
 import Statistics.Function (indices)
 import Statistics.Sample (mean, stdDev, variance, varianceUnbiased)
 import Statistics.Types (Estimator(..), Sample)
-import System.Random.MWC (Gen, initialize, uniform, uniformVector)
+import System.Random.MWC (GenIO, initialize, uniform, uniformVector)
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
@@ -70,7 +70,7 @@
 -- available CPUs.  At least with GHC 7.0, parallel performance seems
 -- best if the parallel garbage collector is disabled (RTS option
 -- @-qg@).
-resample :: Gen (PrimState IO)
+resample :: GenIO
          -> [Estimator]         -- ^ Estimation functions.
          -> Int                 -- ^ Number of resamples to compute.
          -> Sample              -- ^ Original sample.
@@ -83,8 +83,8 @@
           where (q,r) = numResamples `quotRem` numCapabilities
   results <- mapM (const (MU.new numResamples)) ests
   done <- newChan
-  forM_ (zip ixs (tail ixs)) $ \ (start,!end) -> do
-    gen' <- initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))
+  gens <- splitGen numCapabilities gen
+  forM_ (zip3 ixs (tail ixs) gens) $ \ (start,!end,gen') -> do
     forkIO $ do
       let loop k ers | k >= end = writeChan done ()
                      | otherwise = do
@@ -175,3 +175,11 @@
 singletonErr :: String -> a
 singletonErr func = error $
                     "Statistics.Resampling." ++ func ++ ": singleton input"
+
+-- | Split a generator into several that can run independently.
+splitGen :: Int -> GenIO -> IO [GenIO]
+splitGen n gen
+  | n <= 0    = return []
+  | otherwise =
+  fmap (gen:) . replicateM (n-1) $
+  initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.13.1.1
+version:        0.13.2.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
