diff --git a/Numeric/MCMC/Metropolis.hs b/Numeric/MCMC/Metropolis.hs
--- a/Numeric/MCMC/Metropolis.hs
+++ b/Numeric/MCMC/Metropolis.hs
@@ -21,6 +21,7 @@
 
 module Numeric.MCMC.Metropolis (
     mcmc
+  , chain
   , metropolis
 
   -- * Re-exported
@@ -31,7 +32,8 @@
   , MWC.asGenIO
   ) where
 
-import Control.Monad (when)
+import Control.Monad (when, replicateM)
+import Control.Monad.Codensity (lowerCodensity)
 import Control.Monad.Primitive (PrimMonad, PrimState)
 import Control.Monad.Trans.Class (lift)
 import Control.Monad.Trans.State.Strict (execStateT, get, put)
@@ -40,7 +42,7 @@
 #if __GLASGOW_HASKELL__ < 710
 import Data.Traversable (Traversable, traverse)
 #endif
-import Pipes (Producer, yield, (>->), runEffect)
+import Pipes (Producer, Consumer, yield, (>->), runEffect, await)
 import qualified Pipes.Prelude as Pipes (mapM_, take)
 import System.Random.MWC.Probability (Gen, Prob)
 import qualified System.Random.MWC.Probability as MWC
@@ -69,19 +71,52 @@
   accept <- lift (MWC.bernoulli acceptProb)
   when accept (put (Chain chainTarget proposalScore proposal chainTunables))
 
--- A Markov chain driven by the Metropolis transition operator.
-chain
+-- Drive a Markov chain via the Metropolis transition operator.
+drive
   :: (Traversable f, PrimMonad m)
   => Double
   -> Chain (f Double) b
   -> Gen (PrimState m)
-  -> Producer (Chain (f Double) b) m ()
-chain radial = loop where
+  -> Producer (Chain (f Double) b) m c
+drive radial = loop where
   loop state prng = do
     next <- lift (MWC.sample (execStateT (metropolis radial) state) prng)
     yield next
     loop next prng
 
+-- | Trace 'n' iterations of a Markov chain and collect the results in a list.
+--
+-- >>> let rosenbrock [x0, x1] = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+-- >>> results <- withSystemRandom . asGenIO $ chain 3 1 [0, 0] rosenbrock
+-- >>> mapM_ print results
+-- 0.0,0.0
+-- 1.4754117657794871e-2,0.5033208261760778
+-- 3.8379699517007895e-3,0.24627131099479127
+chain
+  :: (PrimMonad m, Traversable f)
+  => Int
+  -> Double
+  -> f Double
+  -> (f Double -> Double)
+  -> Gen (PrimState m)
+  -> m [Chain (f Double) b]
+chain n radial position target gen = runEffect $
+        drive radial origin gen
+    >-> collect n
+  where
+    ctarget = Target target Nothing
+
+    origin = Chain {
+        chainScore    = lTarget ctarget position
+      , chainTunables = Nothing
+      , chainTarget   = ctarget
+      , chainPosition = position
+      }
+
+    collect :: Monad m => Int -> Consumer a m [a]
+    collect size = lowerCodensity $
+      replicateM size (lift Pipes.await)
+
 -- | Trace 'n' iterations of a Markov chain and stream them to stdout.
 --
 -- >>> let rosenbrock [x0, x1] = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
@@ -98,7 +133,7 @@
   -> Gen (PrimState m)
   -> m ()
 mcmc n radial chainPosition target gen = runEffect $
-        chain radial Chain {..} gen
+        drive radial Chain {..} gen
     >-> Pipes.take n
     >-> Pipes.mapM_ (liftIO . print)
   where
diff --git a/mighty-metropolis.cabal b/mighty-metropolis.cabal
--- a/mighty-metropolis.cabal
+++ b/mighty-metropolis.cabal
@@ -1,5 +1,5 @@
 name:                mighty-metropolis
-version:             1.1.0
+version:             1.2.0
 synopsis:            The Metropolis algorithm.
 homepage:            http://github.com/jtobin/mighty-metropolis
 license:             MIT
@@ -15,8 +15,9 @@
   Wander around parameter space according to a simple spherical Gaussian
   distribution.
   .
-  Exports a 'mcmc' function that prints a trace to stdout, as well as a
-  'metropolis' transition operator that can be used more generally.
+  Exports a 'mcmc' function that prints a trace to stdout, a 'chain' function
+  for collecting results in-memory, and a 'metropolis' transition operator that
+  can be used more generally.
   .
   > import Numeric.MCMC.Metropolis
   >
@@ -38,6 +39,7 @@
       Numeric.MCMC.Metropolis
   build-depends:
       base             >= 4 && < 6
+    , kan-extensions   >= 5 && < 6
     , pipes            >= 4 && < 5
     , primitive        >= 0.6 && < 1.0
     , mcmc-types       >= 1.0.1
diff --git a/test/BNN.hs b/test/BNN.hs
--- a/test/BNN.hs
+++ b/test/BNN.hs
@@ -11,5 +11,7 @@
   x1 = index xs 1
 
 main :: IO ()
-main = withSystemRandom . asGenIO $ mcmc 100 1 (fromList [0, 0]) bnn
+main = withSystemRandom . asGenIO $ \gen -> do
+  _ <- chain 50 1 (fromList [0, 0]) bnn gen
+  mcmc 50 1 (fromList [0, 0]) bnn gen
 
diff --git a/test/Rosenbrock.hs b/test/Rosenbrock.hs
--- a/test/Rosenbrock.hs
+++ b/test/Rosenbrock.hs
@@ -8,5 +8,7 @@
 rosenbrock [x0, x1] = negate (5  *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
 
 main :: IO ()
-main = withSystemRandom . asGenIO $ mcmc 100 1 [0, 0] rosenbrock
+main = withSystemRandom . asGenIO $ \gen -> do
+  _ <- chain 50 1 [0, 0] rosenbrock gen
+  mcmc 50 1 [0, 0] rosenbrock gen