Move resampling code to the ST monad.

Statistics/Function.hs

@@ -1,4 +1,4 @@
-{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE Rank2Types, TypeOperators #-}
 -- |
 -- Module    : Statistics.Function
 -- Copyright : (c) 2009 Bryan O'Sullivan
@@ -19,7 +19,7 @@ module Statistics.Function
     ) where
 
 import Control.Exception (assert)
-import Control.Monad.ST (unsafeSTToIO)
+import Control.Monad.ST (ST)
 import Data.Array.Vector.Algorithms.Combinators (apply)
 import Data.Array.Vector
 import qualified Data.Array.Vector.Algorithms.Intro as I
@@ -49,12 +49,12 @@ minMax = fini . foldlU go (MM (1/0) (-1/0))
 {-# INLINE minMax #-}
 
 -- | Create an array, using the given action to populate each element.
-createU :: (UA e) => Int -> (Int -> IO e) -> IO (UArr e)
+createU :: (UA e) => forall s. Int -> (Int -> ST s e) -> ST s (UArr e)
 createU size itemAt = assert (size >= 0) $
-    unsafeSTToIO (newMU size) >>= loop 0
+    newMU size >>= loop 0
   where
-    loop k arr | k >= size = unsafeSTToIO (unsafeFreezeAllMU arr)
+    loop k arr | k >= size = unsafeFreezeAllMU arr
                | otherwise = do
       r <- itemAt k
-      unsafeSTToIO (writeMU arr k r)
+      writeMU arr k r
       loop (k+1) arr

Statistics/Resampling.hs

@@ -17,12 +17,12 @@ module Statistics.Resampling
     ) where
 
 import Control.Monad (forM_)
-import Control.Monad.ST (unsafeSTToIO)
+import Control.Monad.ST (ST)
 import Data.Array.Vector
 import Data.Array.Vector.Algorithms.Intro (sort)
 import Statistics.Function (createU)
+import Statistics.RandomVariate (Gen, uniform)
 import Statistics.Types (Estimator, Sample)
-import System.Random.Mersenne (MTGen, random)
 
 -- | A resample drawn randomly, with replacement, from a set of data
 -- points.  Distinct from a normal array to make it harder for your
@@ -33,20 +33,19 @@ newtype Resample = Resample {
 
 -- | Resample a data set repeatedly, with replacement, computing each
 -- estimate over the resampled data.
-resample :: MTGen -> [Estimator] -> Int -> Sample -> IO [Resample]
+resample :: Gen s -> [Estimator] -> Int -> Sample -> ST s [Resample]
 resample gen ests numResamples samples = do
-  results <- unsafeSTToIO . mapM (const (newMU numResamples)) $ ests
+  results <- mapM (const (newMU numResamples)) $ ests
   loop 0 (zip ests results)
-  unsafeSTToIO $ do
-    mapM_ sort results
-    mapM (fmap Resample . unsafeFreezeAllMU) results
+  mapM_ sort results
+  mapM (fmap Resample . unsafeFreezeAllMU) results
  where
   loop k ers | k >= numResamples = return ()
              | otherwise = do
     re <- createU n $ \_ -> do
-            r <- random gen
+            r <- uniform gen
             return (indexU samples (abs r `mod` n))
-    unsafeSTToIO . forM_ ers $ \(est,arr) ->
+    forM_ ers $ \(est,arr) ->
         writeMU arr k . est $ re
     loop (k+1) ers
   n = lengthU samples

statistics.cabal

@@ -55,7 +55,6 @@ library
   build-depends:
     base < 5,
     erf,
-    mersenne-random,
     time,
     uvector >= 0.1.0.4,
     uvector-algorithms >= 0.2