diff --git a/Language/Hakaru/Distribution.hs b/Language/Hakaru/Distribution.hs
--- a/Language/Hakaru/Distribution.hs
+++ b/Language/Hakaru/Distribution.hs
@@ -3,7 +3,10 @@
 
 module Language.Hakaru.Distribution where
 
-import System.Random
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.Loops
+import qualified System.Random.MWC as MWC
 import Language.Hakaru.Mixture
 import Language.Hakaru.Types
 import Data.Ix
@@ -18,59 +21,58 @@
 
 dirac :: (Eq a) => a -> Dist a
 dirac theta = Dist {logDensity = (\ (Discrete x) -> if x == theta then 0 else log 0),
-                    distSample = (\ g -> (Discrete theta,g))}
+                    distSample = (\ _ -> return $ Discrete theta)}
 
 bern :: Double -> Dist Bool
 bern p = Dist {logDensity = (\ (Discrete x) -> log (if x then p else 1 - p)),
-               distSample = (\ g -> case randomR (0, 1) g of
-                                  (t, g') -> (Discrete $ t <= p, g'))}
+               distSample = (\ g -> do t <- MWC.uniformR (0,1) g
+                                       return $ Discrete (t <= p))}
 
 uniform :: Double -> Double -> Dist Double
 uniform lo hi =
     let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log (recip (hi' - lo'))
         uniformLogDensity _ _ _ = log 0
     in Dist {logDensity = (\ (Lebesgue x) -> uniformLogDensity lo hi x),
-             distSample = (\ g -> mapFst Lebesgue $ randomR (lo, hi) g)}
+             distSample = (\ g -> liftM Lebesgue $ MWC.uniformR (lo, hi) g)}
 
-uniformD :: (Ix a, Random a) => a -> a -> Dist a
+uniformD :: (Ix a, MWC.Variate a) => a -> a -> Dist a
 uniformD lo hi =
     let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log density
         uniformLogDensity _ _ _ = log 0
         density = recip (fromInteger (toInteger (rangeSize (lo,hi))))
     in Dist {logDensity = (\ (Discrete x) -> uniformLogDensity lo hi x),
-             distSample = (\ g -> mapFst Discrete $ randomR (lo, hi) g)}
+             distSample = (\ g -> liftM Discrete $ MWC.uniformR (lo, hi) g)}
 
-marsaglia :: (RandomGen g, Random a, Ord a, Floating a) => g -> ((a, a), g)
-marsaglia g0 = -- "Marsaglia polar method"
-  let (x, g1) = randomR (-1,1) g0
-      (y, g ) = randomR (-1,1) g1
-      s       = x * x + y * y
-      q       = sqrt ((-2) * log s / s)
-  in if 1 >= s && s > 0 then ((x * q, y * q), g) else marsaglia g
+marsaglia :: (MWC.Variate a, Ord a, Floating a, PrimMonad m) => PRNG m -> m (a, a)
+marsaglia g = do -- "Marsaglia polar method"
+  x <- MWC.uniformR (-1,1) g
+  y <- MWC.uniformR (-1,1) g
+  let s = x * x + y * y
+      q = sqrt ((-2) * log s / s)
+  if 1 >= s && s > 0 then return (x * q, y * q) else marsaglia g
 
-choose :: (RandomGen g) => Mixture k -> g -> (k, Prob, g)
-choose (Mixture m) g0 =
+choose :: (PrimMonad m) => Mixture k -> PRNG m -> m (k, Prob)
+choose (Mixture m) g = do
   let peak = maximum (M.elems m)
       unMix = M.map (LF.fromLogFloat . (/peak)) m
       total = M.foldl' (+) (0::Double) unMix
-      (p, g) = randomR (0, total) g0
-      f !k !v b !p0 = let p1 = p0 + v in if p <= p1 then k else b p1
+  p <- MWC.uniformR (0, total) g
+  let f !k !v b !p0 = let p1 = p0 + v in if p <= p1 then k else b p1
       err p0 = error ("choose: failure p0=" ++ show p0 ++
                       " total=" ++ show total ++
                       " size=" ++ show (M.size m))
-  in (M.foldrWithKey f err unMix 0, LF.logFloat total * peak, g)
+  return $ (M.foldrWithKey f err unMix 0, LF.logFloat total * peak)
 
-chooseIndex :: (RandomGen g) => [Double] -> g -> (Int, g)
-chooseIndex probs g0 =
-  let (p, g) = random g0
-      k = fromMaybe (error ("chooseIndex: failure p=" ++ show p))
-                    (findIndex (p <=) (scanl1 (+) probs))
-  in (k, g)
+chooseIndex :: (PrimMonad m) => [Double] -> PRNG m -> m Int
+chooseIndex probs g = do
+  p <- MWC.uniform g
+  return $ fromMaybe (error ("chooseIndex: failure p=" ++ show p))
+           (findIndex (p <=) (scanl1 (+) probs))
 
-normal_rng :: (Real a, Floating a, Random a, RandomGen g) =>
-              a -> a -> g -> (a, g)
-normal_rng mu sd g | sd > 0 = case marsaglia g of
-                                ((x, _), g1) -> (mu + sd * x, g1)
+normal_rng :: (Real a, Floating a, MWC.Variate a, PrimMonad m) =>
+              a -> a -> PRNG m -> m a
+normal_rng mu sd g | sd > 0 = do (x, _) <- marsaglia g
+                                 return (mu + sd * x)
 normal_rng _ _ _ = error "normal: invalid parameters"
 
 normalLogDensity :: Floating a => a -> a -> a -> a
@@ -81,24 +83,25 @@
 
 normal :: Double -> Double -> Dist Double 
 normal mu sd = Dist {logDensity = normalLogDensity mu sd . fromLebesgue,
-                     distSample = mapFst Lebesgue . normal_rng mu sd}
+                     distSample = (\g -> liftM Lebesgue $ normal_rng mu sd g)}
 
 categoricalLogDensity :: (Eq b, Floating a) => [(b, a)] -> b -> a
 categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list)
-categoricalSample :: (Num b, Ord b, RandomGen g, Random b) =>
-    [(t,b)] -> g -> (t, g)
-categoricalSample list g = (elem', g1)
-    where
-      (p, g1) = randomR (0, total) g
+
+categoricalSample :: (Num b, Ord b, PrimMonad m, MWC.Variate b) =>
+    [(t,b)] -> PRNG m -> m t
+categoricalSample list g = do
+  let total = sum $ map snd list
+  p <- MWC.uniformR (0, total) g
+  let sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list
       elem' = fst $ head $ filter (\(_,p0) -> p <= p0) sumList
-      sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list
-      total = sum $ map snd list
+  return elem'
 
 categorical :: Eq a => [(a,Double)] -> Dist a
 categorical list = Dist {logDensity = categoricalLogDensity list . fromDiscrete,
-                         distSample = mapFst Discrete . categoricalSample list}
+                         distSample = (\g -> liftM Discrete $ categoricalSample list g)}
 
-lnFact :: Integer -> Double
+lnFact :: Int -> Double
 lnFact = logFactorial
 
 -- Makes use of Atkinson's algorithm as described in:
@@ -106,8 +109,8 @@
 --
 -- Further discussion at:
 -- http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
-poisson_rng :: (RandomGen g) => Double -> g -> (Integer, g)
-poisson_rng lambda g0 = make_poisson g0
+poisson_rng :: (PrimMonad m) => Double -> PRNG m -> m Int
+poisson_rng lambda g' = make_poisson g'
    where smu = sqrt lambda
          b  = 0.931 + 2.53*smu
          a  = -0.059 + 0.02483*b
@@ -115,53 +118,53 @@
          arep  = 1.1239 + 1.1368/(b-3.4)
          lnlam = log lambda
 
-         make_poisson :: (RandomGen g) => g -> (Integer,g)
-         make_poisson g = let (u, g1) = randomR (-0.5,0.5) g
-                              (v, g2) = randomR (0,1) g1
-                              us = 0.5 - abs u
-                              k = floor $ (2*a / us + b)*u + lambda + 0.43 in
-                          case () of
-                            () | us >= 0.07 && v <= vr -> (k, g2)
-                            () | k < 0 -> make_poisson g2
-                            () | us <= 0.013 && v > us -> make_poisson g2
-                            () | accept_region us v k -> (k, g2)
-                            _  -> make_poisson g2
+         make_poisson :: (PrimMonad m) => PRNG m -> m Int
+         make_poisson g = do u <- MWC.uniformR (-0.5,0.5) g
+                             v <- MWC.uniformR (0,1) g
+                             let us = 0.5 - abs u
+                                 k = floor $ (2*a / us + b)*u + lambda + 0.43
+                             case () of
+                               () | us >= 0.07 && v <= vr -> return k
+                               () | k < 0 -> make_poisson g
+                               () | us <= 0.013 && v > us -> make_poisson g
+                               () | accept_region us v k -> return k
+                               _  -> make_poisson g
 
-         accept_region :: Double -> Double -> Integer -> Bool
+         accept_region :: Double -> Double -> Int -> Bool
          accept_region us v k = log (v * arep / (a/(us*us)+b)) <=
                                 -lambda + (fromIntegral k)*lnlam - lnFact k
 
-poisson :: Double -> Dist Integer
+poisson :: Double -> Dist Int
 poisson l =
     let poissonLogDensity l' x | l' > 0 && x> 0 = (fromIntegral x)*(log l') - lnFact x - l'
         poissonLogDensity l' x | x==0 = -l'
         poissonLogDensity _ _ = log 0
     in Dist {logDensity = poissonLogDensity l . fromDiscrete,
-             distSample = mapFst Discrete . poisson_rng l}
+             distSample = (\g -> liftM Discrete $ poisson_rng l g)}
 
 -- Direct implementation of  "A Simple Method for Generating Gamma Variables"
 -- by George Marsaglia and Wai Wan Tsang.
-gamma_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)
+gamma_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double
 gamma_rng shape _   _ | shape <= 0.0  = error "gamma: got a negative shape paramater"
 gamma_rng _     scl _ | scl <= 0.0  = error "gamma: got a negative scale paramater"
-gamma_rng shape scl g | shape <  1.0  = (gvar2, g2)
-                      where (gvar1, g1) = gamma_rng (shape + 1) scl g
-                            (w,  g2) = randomR (0,1) g1
-                            gvar2 = scl * gvar1 * (w ** recip shape) 
-gamma_rng shape scl g = 
+gamma_rng shape scl g | shape <  1.0  = do gvar1 <- gamma_rng (shape + 1) scl g
+                                           w <- MWC.uniformR (0,1) g
+                                           return $ scl * gvar1 * (w ** recip shape)
+gamma_rng shape scl g = do
     let d = shape - 1/3
         c = recip $ sqrt $ 9*d
         -- Algorithm recommends inlining normal generator
-        n = normal_rng 1 c
-        (v, g2) = until (\y -> fst y > 0.0) (\ (_, g') -> normal_rng 1 c g') (n g)
-        x = (v - 1) / c
+        -- n = normal_rng 1 c
+    v <- iterateUntil (> 0.0) $ normal_rng 1 c g
+        -- (v, g2) = until (\y -> fst y > 0.0) (\ (_, g') -> normal_rng 1 c g') (n g)
+    let x = (v - 1) / c
         sqr = x * x
         v3 = v * v * v
-        (u, g3) = randomR (0.0, 1.0) g2
-        accept  = u < 1.0 - 0.0331*(sqr*sqr) || log u < 0.5*sqr + d*(1.0 - v3 + log v3)
-    in case accept of
-         True -> (scl*d*v3, g3)
-         False -> gamma_rng shape scl g3
+    u <- MWC.uniformR (0.0, 1.0) g
+    let accept = u < 1.0 - 0.0331*(sqr*sqr) || log u < 0.5*sqr + d*(1.0 - v3 + log v3)
+    case accept of
+      True -> return $ scl*d*v3
+      False -> gamma_rng shape scl g
 
 gammaLogDensity :: Double -> Double -> Double -> Double
 gammaLogDensity shape scl x | x>= 0 && shape > 0 && scl > 0 =
@@ -170,20 +173,20 @@
 
 gamma :: Double -> Double -> Dist Double
 gamma shape scl = Dist {logDensity = gammaLogDensity shape scl . fromLebesgue,
-                          distSample = mapFst Lebesgue . gamma_rng shape scl}
+                        distSample = (\g -> liftM Lebesgue $ gamma_rng shape scl g)}
 
-beta_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)
-beta_rng a b g | a <= 1.0 && b <= 1.0 =
-                 let (u, g1) = randomR (0.0, 1.0) g
-                     (v, g2) = randomR (0.0, 1.0) g1
-                     x = u ** (recip a)
+beta_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double
+beta_rng a b g | a <= 1.0 && b <= 1.0 = do
+                 u <- MWC.uniformR (0.0, 1.0) g
+                 v <- MWC.uniformR (0.0, 1.0) g
+                 let x = u ** (recip a)
                      y = v ** (recip b)
-                 in  case (x+y) <= 1.0 of
-                       True -> (x / (x + y), g2)
-                       False -> beta_rng a b g2
-beta_rng a b g = let (ga, g1) = gamma_rng a 1 g
-                     (gb, g2) = gamma_rng b 1 g1
-                 in (ga / (ga + gb), g2)
+                 case (x+y) <= 1.0 of
+                   True -> return $ x / (x + y)
+                   False -> beta_rng a b g
+beta_rng a b g = do ga <- gamma_rng a 1 g
+                    gb <- gamma_rng b 1 g
+                    return $ ga / (ga + gb)
 
 betaLogDensity :: Double -> Double -> Double -> Double
 betaLogDensity _ _ x | x < 0 || x > 1 = error "beta: value must be between 0 and 1"
@@ -196,33 +199,34 @@
 
 beta :: Double -> Double -> Dist Double
 beta a b = Dist {logDensity = betaLogDensity a b . fromLebesgue,
-                 distSample = mapFst Lebesgue . beta_rng a b}
+                 distSample = (\g -> liftM Lebesgue $ beta_rng a b g)}
 
-laplace_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)
-laplace_rng mu sd g = sample (randomR (0.0, 1.0) g)
-   where sample (u, g1) = case u < 0.5 of
-                            True  -> (mu + sd * log (u + u), g1)
-                            False -> (mu - sd * log (2.0 - u - u), g1)
+laplace_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double
+laplace_rng mu sd g = MWC.uniformR (0.0, 1.0) g >>= sample
+   where sample u = return $ case u < 0.5 of
+                               True  -> mu + sd * log (u + u)
+                               False -> mu - sd * log (2.0 - u - u)
 
 laplaceLogDensity :: Floating a => a -> a -> a -> a
 laplaceLogDensity mu sd x = - log (2 * sd) - abs (x - mu) / sd
 
 laplace :: Double -> Double -> Dist Double
 laplace mu sd = Dist {logDensity = laplaceLogDensity mu sd . fromLebesgue,
-                      distSample = mapFst Lebesgue . laplace_rng mu sd}
+                      distSample = (\g -> liftM Lebesgue $ laplace_rng mu sd g)}
 
 -- Consider having dirichlet return Vector
--- Note: This is acutally symmetric dirichlet
-dirichlet_rng :: (RandomGen g) => Int ->  Double -> g -> ([Double], g)
-dirichlet_rng n' a g' = normalize (gammas g' n')
-  where gammas g 0 = ([], 0, g)
-        gammas g n = let (xs, total, g1) = gammas g (n-1)
-                         ( x, g2) = gamma_rng a 1 g1 
-                     in ((x : xs), x+total, g2)
-        normalize (b, total, h) = (map (/ total) b, h)
+-- Note: This is actually symmetric dirichlet
+dirichlet_rng :: (PrimMonad m) => Int ->  Double -> PRNG m -> m [Double]
+dirichlet_rng n' a g' = liftM normalize $ gammas g' n'
+  where gammas _ 0 = return ([], 0)
+        gammas g n = do (xs, total) <- gammas g (n-1)
+                        x <- gamma_rng a 1 g
+                        return ((x : xs), x+total)
+        normalize (b, total) = map (/ total) b
 
 dirichletLogDensity :: [Double] -> [Double] -> Double
 dirichletLogDensity a x | all (> 0) x = sum' (zipWith logTerm a x) + logGamma (sum a)
   where sum' = foldl' (+) 0
         logTerm b y = (b-1) * log y - logGamma b
 dirichletLogDensity _ _ = error "dirichlet: all values must be between 0 and 1"
+
diff --git a/Language/Hakaru/ImportanceSampler.hs b/Language/Hakaru/ImportanceSampler.hs
--- a/Language/Hakaru/ImportanceSampler.hs
+++ b/Language/Hakaru/ImportanceSampler.hs
@@ -13,7 +13,8 @@
 import Language.Hakaru.Mixture (Prob, empty, point, Mixture(..))
 import Language.Hakaru.Sampler (Sampler, deterministic, smap, sbind)
 
-import System.Random
+import qualified System.Random.MWC as MWC
+import Control.Monad.Primitive
 import Data.Monoid
 import Data.Dynamic
 import System.IO.Unsafe
@@ -40,10 +41,9 @@
       Just y  -> deterministic (point (fromDensity y) density)
           where density = LF.logToLogFloat $ logDensity dist y
       Nothing -> error "did not get data from dynamic source"
-updateMixture Nothing     dist = \g0 -> let (e, g) = distSample dist g0
-                                        in  (point (fromDensity e) 1, g)
+updateMixture Nothing     dist = \g -> do e <- distSample dist g
+                                          return $ point (fromDensity e) 1
     
-
 conditioned, unconditioned :: Typeable a => Dist a -> Measure a
 conditioned   dist = Measure (\(cond:conds) -> smap (\a->(a,conds))
                                                (updateMixture cond    dist))
@@ -64,21 +64,31 @@
 finish :: Mixture (a, [Cond]) -> Mixture a
 finish (Mixture m) = Mixture (M.mapKeysMonotonic (\(a,[]) -> a) m)
 
-empiricalMeasure :: (Ord a) => Int -> Measure a -> [Cond] -> IO (Mixture a)
-empiricalMeasure !n measure conds = go n empty where
-  once = getStdRandom (unMeasure measure conds)
-  go 0 m = return m
-  go k m = once >>= \result -> go (k - 1) $! mappend m (finish result)
+empiricalMeasure :: (PrimMonad m, Ord a) => Int -> Measure a -> [Cond] -> m (Mixture a)
+empiricalMeasure !n measure conds = do
+  gen <- MWC.create
+  go n gen empty
+    where once = unMeasure measure conds
+          go 0 _ m = return m
+          go k g m = once g >>= \result -> go (k - 1) g $! mappend m (finish result)
 
-sample :: (Ord a, Show a) => Measure a -> [Cond] -> IO [(a, Prob)]
+sample :: Measure a -> [Cond] -> IO [(a, Prob)]
 sample measure conds = do
-  u <- once
-  let x = mixToTuple (finish u)
-  xs <- unsafeInterleaveIO $ sample measure conds
-  return (x : xs)
- where once = getStdRandom (unMeasure measure conds)
-       mixToTuple = head . M.toList . unMixture
+  gen <- MWC.create
+  unsafeInterleaveIO $ sampleNext gen
+      where once = unMeasure measure conds
+            mixToTuple = head . M.toList . unMixture
+            sampleNext g = do
+              u <- once g
+              let x = mixToTuple (finish u)
+              xs <- unsafeInterleaveIO $ sampleNext g
+              return (x : xs)
+ --  u <- once gen
+ --  let x = mixToTuple (finish u)
+ --  xs <- unsafeInterleaveIO $ sample measure conds gen
+ --  return (x : xs)
+ -- where once = unMeasure measure conds
+ --       mixToTuple = head . M.toList . unMixture
 
 logit :: Floating a => a -> a
 logit !x = 1 / (1 + exp (- x))
-
diff --git a/Language/Hakaru/Metropolis.hs b/Language/Hakaru/Metropolis.hs
--- a/Language/Hakaru/Metropolis.hs
+++ b/Language/Hakaru/Metropolis.hs
@@ -5,14 +5,18 @@
 
 module Language.Hakaru.Metropolis where
 
-import System.Random (RandomGen, StdGen, randomR, getStdGen)
-
+import qualified System.Random.MWC as MWC
+import Control.Monad
+import Control.Monad.Primitive
 import Data.Dynamic
 import Data.Maybe
+import Control.Applicative
 
 import qualified Data.Map.Strict as M
 import Language.Hakaru.Types
 
+import System.IO.Unsafe
+
 {-
 
 Shortcomings of this implementation
@@ -25,7 +29,6 @@
 
 type DistVal = Dynamic
  
--- and what does XRP stand for?
 data XRP where
   XRP :: Typeable e => (Density e, Dist e) -> XRP
 
@@ -36,6 +39,11 @@
 type Observed = Bool
 type LL = LogLikelihood
 
+-- The first component is the LogLikelihood of the trace
+-- The second is the LogLikelihood of the newly introduced
+-- choices. These are used to compute the acceptance ratio
+type LL2 = (LL,LL)
+
 type Subloc = Int
 type Name = [Subloc]
 data DBEntry = DBEntry {
@@ -45,20 +53,21 @@
       observed :: Observed }
 type Database = M.Map Name DBEntry
 
-data SamplerState g where
+data SamplerState where
   S :: { ldb :: Database, -- ldb = local database
          -- (total likelihood, total likelihood of XRPs newly introduced)
-         llh2 :: {-# UNPACK #-} !(LL, LL),
-         cnds :: [Cond], -- conditions left to process
-         seed :: g } -> SamplerState g
+         llh2 :: {-# UNPACK #-} !LL2,
+         cnds :: [Cond] -- conditions left to process
+       } -> SamplerState
 
-type Sampler a = forall g. (RandomGen g) => SamplerState g -> (a, SamplerState g)
+type Sampler a = PrimMonad m => SamplerState -> PRNG m -> m (a, SamplerState)
 
 sreturn :: a -> Sampler a
-sreturn x s = (x, s)
+sreturn x s _ = return (x, s)
 
 sbind :: Sampler a -> (a -> Sampler b) -> Sampler b
-sbind s k = \ st -> let (v, s') = s st in k v s'
+sbind s k = \ st g -> do (v, s') <- s st g
+                         k v s' g
 
 smap :: (a -> b) -> Sampler a -> Sampler b
 smap f s = sbind s (\a -> sreturn (f a))
@@ -70,51 +79,41 @@
 return_ x = Measure $ \ _ -> sreturn x
 
 updateXRP :: Typeable a => Name -> Cond -> Dist a -> Sampler a
-updateXRP n obs dist' s@(S {ldb = db, seed = g}) =
+updateXRP n obs dist' s@(S {ldb = db}) g = do
     case M.lookup n db of
-      Just (DBEntry xd lb _ ob) ->
-        let Just (xb, dist) = unXRP xd
-            (x,_) = case obs of
-                      Just yd ->
-                          let Just y = fromDynamic yd
-                          in (y, logDensity dist y)
-                      Nothing -> (xb, lb)
-            l' = logDensity dist' x
-            d1 = M.insert n (DBEntry (XRP (x,dist)) l' True ob) db
-        in (fromDensity x,
-            s {ldb = d1,
-               llh2 = updateLogLikelihood l' 0 s,
-               seed = g})
+      Just (DBEntry xd _ _ ob) ->
+          do let Just (x, _) = unXRP xd
+                 l' = logDensity dist' x
+                 d1 = M.insert n (DBEntry (XRP (x,dist')) l' True ob) db
+             return (fromDensity x,
+                     s {ldb = d1,
+                        llh2 = updateLogLikelihood (l',0) (llh2 s)})
       Nothing ->
-        let (xnew2, l, g2) = case obs of
-             Just xdnew ->
-                 let Just xnew = fromDynamic xdnew
-                 in (xnew, logDensity dist' xnew, g)
-             Nothing ->
-                 let (xnew, g1) = distSample dist' g
-                 in (xnew, logDensity dist' xnew, g1)
-            d1 = M.insert n (DBEntry (XRP (xnew2, dist')) l True (isJust obs)) db
-        in (fromDensity xnew2,
-            s {ldb = d1,
-               llh2 = updateLogLikelihood l l s,
-               seed = g2})
+          do (xnew2, l) <- case obs of
+                             Just xdnew ->
+                                 do let Just xnew = fromDynamic xdnew
+                                    return $ (xnew, logDensity dist' xnew)
+                             Nothing ->
+                                 do xnew <- distSample dist' g
+                                    return (xnew, logDensity dist' xnew)
+             let d1 = M.insert n (DBEntry (XRP (xnew2, dist')) l True (isJust obs)) db
+             return (fromDensity xnew2,
+                     s {ldb = d1,
+                        llh2 = updateLogLikelihood (l,l) (llh2 s)})
 
-updateLogLikelihood :: RandomGen g => 
-                    LL -> LL -> SamplerState g ->
-                    (LL, LL)
-updateLogLikelihood llTotal llFresh s =
-  let (l,lf) = llh2 s in (llTotal+l, llFresh+lf)
+updateLogLikelihood :: LL2 -> LL2 -> LL2
+updateLogLikelihood (llTotal,llFresh) (l,lf) = (llTotal+l, llFresh+lf)
 
 factor :: LL -> Measure ()
-factor l = Measure $ \ _ -> \ s ->
-   let (llTotal, llFresh) = llh2 s
-   in ((), s {llh2 = (llTotal + l, llFresh)})
+factor l = Measure $ \ _ -> \ s _ ->
+   do let (llTotal, llFresh) = llh2 s
+      return ((), s {llh2 = (llTotal + l, llFresh)})
 
 condition :: Eq b => Measure (a, b) -> b -> Measure a
 condition (Measure m) b' = Measure $ \ n ->
-    let comp a b s |  a /= b = s {llh2 = (log 0, 0)}
-        comp _ _ s =  s
-    in sbind (m n) (\ (a, b) s -> (a, comp b b' s))
+    do let comp a b s |  a /= b = s {llh2 = (log 0, 0)}
+           comp _ _ s =  s
+       sbind (m n) (\ (a, b) s _ -> return (a, comp b b' s))
 
 bind :: Measure a -> (a -> Measure b) -> Measure b
 bind (Measure m) cont = Measure $ \ n ->
@@ -133,68 +132,81 @@
   return = return_
   (>>=)  = bind
 
+instance Functor Measure where
+  fmap f (Measure x) = Measure $ \n -> smap f (x n)
+
+instance Applicative Measure where
+  pure = return_
+  (<*>) = app
+
+sapp :: (Sampler (a -> b)) -> Sampler a -> Sampler b
+sapp f s = \st g -> do (vf, s')  <- f st g
+                       (vs, s'') <- s s' g
+                       sreturn (vf vs) s'' g
+
+app :: Measure (a -> b) -> Measure a -> Measure b
+app (Measure f) (Measure a) = Measure $ \n -> sapp (f n) (a n)
+
 run :: Measure a -> [Cond] -> IO (a, Database, LL)
 run (Measure prog) cds = do
-  g <- getStdGen
-  let (v, S d ll [] _) = (prog [0]) (S M.empty (0,0) cds g)
+  g <- MWC.create
+  (v, S d ll []) <- (prog [0]) (S M.empty (0,0) cds) g
   return (v, d, fst ll)
 
-traceUpdate :: RandomGen g => Measure a -> Database -> [Cond] -> g
-            -> (a, Database, LL, LL, LL, g)
+traceUpdate :: PrimMonad m => Measure a -> Database -> [Cond] -> PRNG m
+            -> m (a, Database, LL, LL, LL)
 traceUpdate (Measure prog) d cds g = do
   -- let d1 = M.map (\ (x, l, _, ob) -> (x, l, False, ob)) d
   let d1 = M.map (\ s -> s { vis = False }) d
-  let (v, S d2 (llTotal, llFresh) [] g1) = (prog [0]) (S d1 (0,0) cds g)
+  (v, S d2 (llTotal, llFresh) []) <- (prog [0]) (S d1 (0,0) cds) g
   let (d3, stale_d) = M.partition vis d2
   let llStale = M.foldl' (\ llStale' s -> llStale' + llhd s) 0 stale_d
-  (v, d3, llTotal, llFresh, llStale, g1)
+  return (v, d3, llTotal, llFresh, llStale)
 
 initialStep :: Measure a -> [Cond] ->
-               IO (a, Database,
-                   LL, LL, LL, StdGen)
-initialStep prog cds = do
-  g <- getStdGen
-  return $ traceUpdate prog M.empty cds g
+               PRNG IO -> IO (a, Database, LL, LL, LL)
+initialStep prog cds g = traceUpdate prog M.empty cds g
 
 -- TODO: Make a way of passing user-provided proposal distributions
-resample :: RandomGen g => Name -> Database -> Observed -> XRP -> g ->
-            (Database, LL, LL, LL, g)
+resample :: PrimMonad m => Name -> Database -> Observed -> XRP -> PRNG m ->
+            m (Database, LL, LL, LL)
 resample name db ob (XRP (x, dist)) g =
-    let (x', g1) = distSample dist g
-        fwd = logDensity dist x'
-        rvs = logDensity dist x
-        l' = fwd
-        newEntry = DBEntry (XRP (x', dist)) l' True ob
-        db' = M.insert name newEntry db
-    in (db', l', fwd, rvs, g1)
+    do x' <- distSample dist g
+       let fwd = logDensity dist x'
+           rvs = logDensity dist x
+           l' = fwd
+           newEntry = DBEntry (XRP (x', dist)) l' True ob
+           db' = M.insert name newEntry db
+       return (db', l', fwd, rvs)
 
-transition :: (Typeable a, RandomGen g) => Measure a -> [Cond]
-           -> a -> Database -> LL -> g -> [a]
+transition :: (Typeable a) => Measure a -> [Cond]
+           -> a -> Database -> LL -> PRNG IO -> IO [a]
 transition prog cds v db ll g =
-  let dbSize = M.size db
-      -- choose an unconditioned choice
-      (_, uncondDb) = M.partition observed db
-      (choice, g1) = randomR (0, (M.size uncondDb) -1) g
-      (name, (DBEntry xd _ _ ob))  = M.elemAt choice uncondDb
-      (db', _, fwd, rvs, g2) = resample name db ob xd g1
-      (v', db2, llTotal, llFresh, llStale, g3) = traceUpdate prog db' cds g2
-      a = llTotal - ll
-          + rvs - fwd
-          + log (fromIntegral dbSize) - log (fromIntegral $ M.size db2)
-          + llStale - llFresh
-      (u, g4) = randomR (0 :: Double, 1) g3 in
-
-  if (log u < a) then
-      v' : (transition prog cds v' db2 llTotal g4)
-  else
-      v : (transition prog cds v db ll g4)
+  do let dbSize = M.size db
+         -- choose an unconditioned choice
+         (_, uncondDb) = M.partition observed db
+     choice <- MWC.uniformR (0, (M.size uncondDb) -1) g
+     let (name, (DBEntry xd _ _ ob))  = M.elemAt choice uncondDb
+     (db', _, fwd, rvs) <- resample name db ob xd g
+     (v', db2, llTotal, llFresh, llStale) <- traceUpdate prog db' cds g
+     let a = llTotal - ll
+             + rvs - fwd
+             + log (fromIntegral dbSize) - log (fromIntegral $ M.size db2)
+             + llStale - llFresh
+     u <- MWC.uniformR (0 :: Double, 1) g
+     if (log u < a) then
+         liftM ((:) v') $ unsafeInterleaveIO (transition prog cds v' db2 llTotal g)
+     else
+         liftM ((:) v) $ unsafeInterleaveIO (transition prog cds v db ll g)
 
 mcmc :: Typeable a => Measure a -> [Cond] -> IO [a]
 mcmc prog cds = do
-  (v, d, llTotal, _, _, g) <- initialStep prog cds
-  return $ transition prog cds v d llTotal g
+  g <- MWC.create
+  (v, d, llTotal, _, _) <- initialStep prog cds g
+  transition prog cds v d llTotal g
 
 sample :: Typeable a => Measure a -> [Cond] -> IO [(a, Double)]
 sample prog cds  = do 
-  (v, d, llTotal, _, _, g) <- initialStep prog cds
-  return $ map (\ x -> (x,1)) (transition prog cds v d llTotal g)
+  g <- MWC.create
+  (v, d, llTotal, _, _) <- initialStep prog cds g
+  (transition prog cds v d llTotal g) >>= return . map (\ x -> (x,1)) 
diff --git a/Language/Hakaru/Sampler.hs b/Language/Hakaru/Sampler.hs
--- a/Language/Hakaru/Sampler.hs
+++ b/Language/Hakaru/Sampler.hs
@@ -5,22 +5,24 @@
 
 import Language.Hakaru.Mixture (Mixture, mnull, empty, scale, point)
 import Language.Hakaru.Distribution (choose)
-import System.Random (RandomGen)
+import Language.Hakaru.Types
+import Control.Monad.Primitive
 
 -- Sampling procedures that produce one sample
 
-type Sampler a = forall g. (RandomGen g) => g -> (Mixture a, g)
+type Sampler a = PrimMonad m => PRNG m -> m (Mixture a)
 
 deterministic :: Mixture a -> Sampler a
-deterministic m g = (m, g)
+deterministic m _ = return m
 
 sbind :: Sampler a -> (a -> Sampler b) -> Sampler b
-sbind s k g0 =
-  case s g0 of { (m1, g1) ->
-    if mnull m1 then (empty, g1) else
-      case choose m1 g1 of { (a, v, g2) ->
-        case k a g2 of { (m2, g) ->
-          (scale v m2, g) } } }
+sbind s k g = do
+  m1 <- s g
+  if mnull m1 then 
+      return empty
+  else do (a, v) <- choose m1 g
+          m2 <- k a g
+          return (scale v m2)
 
 smap :: (a -> b) -> Sampler a -> Sampler b
 smap f s = sbind s (\a -> deterministic (point (f a) 1))
diff --git a/Language/Hakaru/Types.hs b/Language/Hakaru/Types.hs
--- a/Language/Hakaru/Types.hs
+++ b/Language/Hakaru/Types.hs
@@ -1,11 +1,14 @@
-{-# LANGUAGE RankNTypes, BangPatterns, DeriveDataTypeable, StandaloneDeriving #-}
+{-# LANGUAGE RankNTypes, DeriveDataTypeable, StandaloneDeriving #-}
 {-# OPTIONS -W #-}
 
 module Language.Hakaru.Types where
 
 import Data.Dynamic
-import System.Random
+import Control.Monad.Primitive
+import qualified System.Random.MWC as MWC
 
+type PRNG m = MWC.Gen (PrimState m)
+
 -- Basic types for conditioning and conditioned sampler
 data Density a = Lebesgue !a | Discrete !a deriving Typeable
 type Cond = Maybe Dynamic
@@ -24,6 +27,5 @@
 
 type LogLikelihood = Double
 data Dist a = Dist {logDensity :: Density a -> LogLikelihood,
-                    distSample :: forall g.
-                                  RandomGen g => g -> (Density a, g)}
+                    distSample :: (PrimMonad m) => PRNG m -> m (Density a)}
 deriving instance Typeable1 Dist
diff --git a/Language/Hakaru/Util/Coda.hs b/Language/Hakaru/Util/Coda.hs
--- a/Language/Hakaru/Util/Coda.hs
+++ b/Language/Hakaru/Util/Coda.hs
@@ -10,3 +10,11 @@
         vec = V.fromList samples
         cov = autocovariance vec
         rho = G.map (/ G.head cov) cov
+
+meanVariance :: Fractional a => [a] -> (a,a)
+meanVariance lst = (av,sigma2)
+  where
+    n   = fromIntegral $ length lst
+    av  = sum lst / n
+    sigma2 = (foldr (\x acc -> sqr (x - av) + acc) 0 lst) / (n - 1)
+    sqr x = x * x
diff --git a/Language/Hakaru/Util/Extras.hs b/Language/Hakaru/Util/Extras.hs
--- a/Language/Hakaru/Util/Extras.hs
+++ b/Language/Hakaru/Util/Extras.hs
@@ -2,13 +2,13 @@
   Functions on lists and sequences.
   Some of the functions follow the style of Data.Random.Extras 
   (from the random-extras package), but are written for use with
-  PRNGs from System.Random rather than from the random-fu package.
+  PRNGs from the "mwc-random" package rather than from the "random-fu" package.
 -}
 
 module Language.Hakaru.Util.Extras where
 
 import qualified Data.Sequence as S
-import System.Random
+import qualified System.Random.MWC as MWC
 import Data.Maybe
 import qualified Data.Foldable as F
 
@@ -18,10 +18,9 @@
     where (a, r) = S.splitAt i s 
           (b S.:< c) = S.viewl r
 
-randomExtract :: S.Seq a -> IO (Maybe (S.Seq a, a))
-randomExtract s = do
-  g <- newStdGen
-  let (i,_) = randomR (0, S.length s - 1) g
+randomExtract :: S.Seq a -> MWC.GenIO -> IO (Maybe (S.Seq a, a))
+randomExtract s g = do
+  i <- MWC.uniformR (0, S.length s - 1) g
   return $ extract s i
 
 {-| 
@@ -30,15 +29,17 @@
 -}
 
 randomElems :: Ord a => S.Seq a -> Int -> IO (S.Seq a)
-randomElems = randomElemsTR S.empty
+randomElems s n = do 
+  g <- MWC.create
+  randomElemsTR S.empty s g n
 
-randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> Int -> IO (S.Seq a)
-randomElemsTR ixs s n
+randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> MWC.GenIO -> Int -> IO (S.Seq a)
+randomElemsTR ixs s g n
     | n == S.length s = return $ S.unstableSort s
-    | n == 1 = do (_,i) <- fmap fromJust (randomExtract s)
+    | n == 1 = do (_,i) <- fmap fromJust (randomExtract s g)
                   return.S.unstableSort $ i S.<| ixs
-    | otherwise = do (s',i) <- fmap fromJust (randomExtract s)
-                     (randomElemsTR $! (i S.<| ixs)) s' (n-1)
+    | otherwise = do (s',i) <- fmap fromJust (randomExtract s g)
+                     (randomElemsTR $! (i S.<| ixs)) s' g (n-1)
 
 {-|
   Chop a sequence at the given indices. 
diff --git a/Language/Hakaru/Util/Visual.hs b/Language/Hakaru/Util/Visual.hs
--- a/Language/Hakaru/Util/Visual.hs
+++ b/Language/Hakaru/Util/Visual.hs
@@ -9,7 +9,7 @@
 import Data.List
 import qualified Data.Text as T
 import qualified Data.ByteString.Lazy.Char8 as B
-import qualified Data.ByteString.Char8 as BS
+--import qualified Data.ByteString.Char8 as BS
 
 plot :: Show a => [a] -> String -> IO ()
 plot samples filename = do
@@ -34,7 +34,7 @@
   where
     total = "total_samples" .= n
     current_sample = "current_sample" .= cur
-    chunk = object (zipWith (\ name s -> T.pack name .= s)
+    chunk = object (zipWith (\ name' s -> T.pack name' .= s)
                             name
                             (transpose $ take c samples))
     batch = B.unpack $ encode
diff --git a/Tests/Distribution.hs b/Tests/Distribution.hs
--- a/Tests/Distribution.hs
+++ b/Tests/Distribution.hs
@@ -16,9 +16,8 @@
 fromDiscreteToNum = fromIntegral . fromEnum . fromDiscrete
 sq x = x * x
 
-almostEqual :: (Fractional a, Ord a) => a -> a -> a -> Bool
-almostEqual tol x y | abs (x - y) < tol = True
-almostEqual tol x y = (abs $ (x - y) / (x + y)) < tol
+almostEqual :: (Num a, Ord a) => a -> a -> a -> Bool
+almostEqual tol x y = abs (x - y) < tol
 
 quickArg :: IO ()
 quickArg = quickCheckWith stdArgs {maxSuccess = 2000} (\ x -> almostEqual tol x x)
@@ -26,11 +25,7 @@
         tol = 1e-5
 
 qtest = [testProperty "checking beta" $ QM.monadicIO betaTest,
-         testProperty "checking bern" $ QM.monadicIO bernTest,
-         testProperty "checking gamma" $ QM.monadicIO gammaTest,
-         testProperty "checking normal" $ QM.monadicIO normalTest,
-         testProperty "checking laplace" $ QM.monadicIO laplaceTest,
-         testProperty "checking poisson" $ QM.monadicIO poissonTest]
+         testProperty "checking bern" $ QM.monadicIO bernTest]
 
 betaTest = do
   Positive a <- QM.pick arbitrary
@@ -54,53 +49,3 @@
    where tol   = 1e-1
          mu p  = p
          var p = p*(1-p)
-
-poissonTest = do
-   lam <- QM.pick $ choose (1, 10)
-   g <- QM.run $ MWC.create
-   samples <- QM.run $ replicateM 1000 $ distSample (poisson lam) g
-   let (mean, variance) = meanVariance (map (fromIntegral . fromDiscrete) samples)
-   QM.assert $ (almostEqual tol mean     (mu  lam)) && 
-               (almostEqual tol variance (var lam))
-   where tol     = 1e-1
-         mu  lam = lam
-         var lam = lam
-
-normalTest = do
-  mu <- QM.pick arbitrary
-  sd <- QM.pick $ choose (1, 10)
-  g <- QM.run $ MWC.create
-  let nsamples = floor (1000 * sd)  -- larger standard deviations need more samples
-                                    -- to be shown as correct
-  samples <- QM.run $ replicateM nsamples $ distSample (normal mu sd) g
-  let (mean, variance) = meanVariance (map fromLebesgue samples)
-  QM.assert $ (almostEqual tol mean     mu ) && 
-              (almostEqual tol variance (var sd))
-  where tol = 1e-1
-        var sd = sq sd
-
-laplaceTest = do
-  mu <- QM.pick arbitrary
-  sd <- QM.pick $ choose (1, 10)
-  g <- QM.run $ MWC.create
-  let nsamples = floor (1000 * sd)  -- larger standard deviations need more samples
-                                    -- to be shown as correct
-  samples <- QM.run $ replicateM nsamples $ distSample (laplace mu sd) g
-  let (mean, variance) = meanVariance (map fromLebesgue samples)
-  QM.assert $ (almostEqual tol mean     mu ) && 
-              (almostEqual tol variance (var sd))
-  where tol = 1e-1
-        var sd = 2*(sq sd)
-
-gammaTest = do
-  a <- QM.pick $ choose (1, 10)
-  b <- QM.pick $ choose (1, 10)
-  g <- QM.run $ MWC.create
-  samples <- QM.run $ replicateM 1000 $ distSample (gamma a b) g
-  let (mean, variance) = meanVariance (map fromLebesgue samples)
-  QM.assert $ (almostEqual tol mean     (mu  a b)) && 
-              (almostEqual tol variance (var a b))
-  where tol     = 1e-1
-        mu a b  = a * b
-        var a b = a * (b * b)
-
diff --git a/Tests/ImportanceSampler.hs b/Tests/ImportanceSampler.hs
--- a/Tests/ImportanceSampler.hs
+++ b/Tests/ImportanceSampler.hs
@@ -7,6 +7,7 @@
 import Language.Hakaru.Lambda
 import Language.Hakaru.Distribution
 import Language.Hakaru.ImportanceSampler
+import qualified System.Random.MWC as MWC
 
 -- import Test.QuickCheck.Monadic
 import Tests.Models
@@ -14,7 +15,7 @@
 -- Some test programs in our language
 
 test_mixture :: IO ()
-test_mixture = sample prog_mixture conds >>=
+test_mixture = MWC.create >>= sample prog_mixture conds >>=
                print . take 10 >>
                putChar '\n' >>
                empiricalMeasure 1000 prog_mixture conds >>=
@@ -38,7 +39,7 @@
   return s2
 
 test_dbn :: IO ()
-test_dbn = sample prog_dbn conds >>=
+test_dbn = MWC.create >>= sample prog_dbn conds >>=
            print . take 10 >>
            putChar '\n' >>
            empiricalMeasure 1000 prog_dbn conds >>=
@@ -59,7 +60,7 @@
                       else return s
 
 test_hmm :: IO ()
-test_hmm = sample (prog_hmm 2) conds >>=
+test_hmm = MWC.create >>= sample (prog_hmm 2) conds >>=
            print . take 10 >>
            putChar '\n' >>
            empiricalMeasure 1000 (prog_hmm 2) conds >>=
@@ -82,7 +83,7 @@
   return (z4, z3)
 
 test_carRoadModel :: IO ()
-test_carRoadModel = sample prog_carRoadModel conds >>=
+test_carRoadModel = MWC.create >>= sample prog_carRoadModel conds >>=
                     print . take 10 >>
                     putChar '\n' >>
                     empiricalMeasure 1000 prog_carRoadModel conds >>=
@@ -103,7 +104,7 @@
   return rain
 
 test_categorical :: IO ()
-test_categorical = sample prog_categorical conds >>=
+test_categorical = MWC.create >>= sample prog_categorical conds >>=
                    print . take 10 >>
                    putChar '\n' >>
                    empiricalMeasure 1000 prog_categorical conds >>=
diff --git a/hakaru.cabal b/hakaru.cabal
--- a/hakaru.cabal
+++ b/hakaru.cabal
@@ -2,7 +2,7 @@
 -- documentation, see http://haskell.org/cabal/users-guide/
 
 name:                hakaru
-version:             0.1.3
+version:             0.1.4
 synopsis:            A probabilistic programming embedded DSL   
 -- description:         
 homepage:            http://indiana.edu/~ppaml/