diff --git a/mcmc.cabal b/mcmc.cabal
--- a/mcmc.cabal
+++ b/mcmc.cabal
@@ -1,6 +1,6 @@
 cabal-version:      3.0
 name:               mcmc
-version:            0.8.1.0
+version:            0.8.2.0
 synopsis:           Sample from a posterior using Markov chain Monte Carlo
 description:
   Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme>
diff --git a/src/Mcmc/MarginalLikelihood.hs b/src/Mcmc/MarginalLikelihood.hs
--- a/src/Mcmc/MarginalLikelihood.hs
+++ b/src/Mcmc/MarginalLikelihood.hs
@@ -21,6 +21,7 @@
   )
 where
 
+import Control.Concurrent (getNumCapabilities)
 import Control.Concurrent.Async hiding (link)
 import Control.Monad
 import Control.Monad.IO.Class
@@ -87,13 +88,15 @@
   { mlAnalysisName :: AnalysisName,
     mlAlgorithm :: MLAlgorithm,
     mlNPoints :: NPoints,
-    -- | Initial burn in at the starting point of the path.
+    -- | Initial burn in at the starting point of the path (or each segment if
+    -- running in parallel).
     mlInitialBurnIn :: BurnInSettings,
     -- | Repetitive burn in at each point on the path.
     mlPointBurnIn :: BurnInSettings,
     -- | The number of iterations performed at each point.
     mlIterations :: Iterations,
     mlExecutionMode :: ExecutionMode,
+    mlParallelizationMode :: ParallelizationMode,
     mlLogMode :: LogMode,
     mlVerbosity :: Verbosity
   }
@@ -139,23 +142,27 @@
 
 sampleAtPoint ::
   ToJSON a =>
+  Bool ->
   Point ->
   Settings ->
   LikelihoodFunction a ->
   MHG a ->
   ML (MHG a)
-sampleAtPoint x ss lhf a = do
+sampleAtPoint isInitialBurnIn x ss lhf a = do
   a'' <- liftIO $ mcmc ss' a'
   let ch'' = fromMHG a''
       ac = acceptances ch''
       mAr = sequence $ acceptanceRates ac
   logDebugB "sampleAtPoint: Summarize cycle."
   logDebugB $ summarizeCycle AllProposals ac $ cycle ch''
-  case mAr of
-    Nothing -> logWarnB "Some acceptance rates are unavailable."
-    Just ar -> do
-      unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1."
-      unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9."
+  unless
+    isInitialBurnIn
+    ( case mAr of
+        Nothing -> logWarnB "Some acceptance rates are unavailable."
+        Just ar -> do
+          unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1."
+          unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9."
+    )
   return a''
   where
     -- For debugging set a proper analysis name.
@@ -181,19 +188,18 @@
 
 traversePoints ::
   ToJSON a =>
-  -- Current point.
-  Int ->
   NPoints ->
-  [Point] ->
+  [(Int, Point)] ->
   Settings ->
   LikelihoodFunction a ->
   MHG a ->
   -- For each point a vector of obtained likelihoods stored in the log domain.
   ML [VU.Vector Likelihood]
-traversePoints _ _ [] _ _ _ = return []
-traversePoints i k (b : bs) ss lhf a = do
-  logInfoS $ "Point " <> show i <> " of " <> show k' <> ": " <> show b <> "."
-  a' <- sampleAtPoint b ss lhf a
+traversePoints _ [] _ _ _ = return []
+traversePoints k ((idb, b) : bs) ss lhf a = do
+  let msg = printf "Point %4d of %4d: %.12f." idb k' b
+  logInfoS msg
+  a' <- sampleAtPoint False b ss lhf a
   -- Get the links samples at this point.
   ls <- liftIO $ takeT n $ trace $ fromMHG a'
   -- Extract the likelihoods.
@@ -205,16 +211,67 @@
   -- resulting in a severe memory leak.
   let !lhs = VU.convert $ VB.map (lhf . state) ls
   -- Sample the other points.
-  lhss <- traversePoints (i + 1) k bs ss lhf a'
+  lhss <- traversePoints k bs ss lhf a'
   return $ lhs : lhss
   where
     n = fromIterations $ sIterations ss
     (NPoints k') = k
 
+nChunks :: Int -> [a] -> [[a]]
+nChunks k xs = chop (chunks k l) xs
+  where
+    l = length xs
+
+chunks :: Int -> Int -> [Int]
+chunks c n = filter (> 0) ns
+  where
+    n' = n `div` c
+    r = n `mod` c
+    ns = replicate r (n' + 1) ++ replicate (c - r) n'
+
+chop :: [Int] -> [a] -> [[a]]
+chop [] [] = []
+chop (n : ns) xs
+  | n > 0 = take n xs : chop ns (drop n xs)
+  | otherwise = error "chop: n negative or zero"
+chop _ _ = error "chop: not all list elements handled"
+
+mlRunPar ::
+  ToJSON a =>
+  ParallelizationMode ->
+  NPoints ->
+  [(Int, Point)] ->
+  ExecutionMode ->
+  Verbosity ->
+  PriorFunction a ->
+  LikelihoodFunction a ->
+  Cycle a ->
+  Monitor a ->
+  a ->
+  StdGen ->
+  ML [VU.Vector Likelihood]
+mlRunPar pm k xs em vb prf lhf cc mn i0 g = do
+  nThreads <- case pm of
+    Sequential -> do
+      logInfoB "Sequential execution."
+      pure 1
+    Parallel -> do
+      n <- liftIO getNumCapabilities
+      logInfoS $ "Parallel execution with " <> show n <> " cores."
+      pure n
+  let xsChunks = nChunks nThreads xs
+  r <- ask
+  xss <-
+    liftIO $
+      mapConcurrently
+        (\thesePoints -> runReaderT (mlRun k thesePoints em vb prf lhf cc mn i0 g) r)
+        xsChunks
+  pure $ concat xss
+
 mlRun ::
   ToJSON a =>
   NPoints ->
-  [Point] ->
+  [(Int, Point)] ->
   ExecutionMode ->
   Verbosity ->
   PriorFunction a ->
@@ -241,12 +298,13 @@
       ssP = Settings nm biP is trLen em Sequential NoSave LogFileOnly vb'
   logDebugB "mlRun: Initialize MHG algorithm."
   a0 <- liftIO $ mhg ssI prf lhf cc mn i0 g
-  logDebugS $ "mlRun: Perform initial burn in at first point " <> show x0 <> "."
-  a1 <- sampleAtPoint x0 ssI lhf a0
+  let msg = printf "Initial burn in at point %.12f with ID %4d." x0 id0
+  logInfoS msg
+  a1 <- sampleAtPoint True x0 ssI lhf a0
   logDebugB "mlRun: Traverse points."
-  traversePoints 1 k xs ssP lhf a1
+  traversePoints k xs ssP lhf a1
   where
-    x0 = head xs
+    (id0, x0) = head xs
 
 -- Use lists since the number of points is expected to be low.
 integrateSimpsonTriangle ::
@@ -280,8 +338,8 @@
   (lhssForward, lhssBackward) <-
     lift $
       concurrently
-        (runReaderT (mlRun k bsForward em vb prf lhf cc mn i0 g0) r)
-        (runReaderT (mlRun k bsBackward em vb prf lhf cc mn i0 g1) r)
+        (runReaderT (mlRunPar pm k (zip [1 ..] bsForward) em vb prf lhf cc mn i0 g0) r)
+        (runReaderT (mlRunPar pm k (zip [1 ..] bsBackward) em vb prf lhf cc mn i0 g1) r)
   logInfoEndTime
 
   logDebugB "tiWrapper: Calculate mean log likelihoods."
@@ -300,6 +358,7 @@
     bsForward = getPoints k
     bsBackward = reverse bsForward
     em = mlExecutionMode s
+    pm = mlParallelizationMode s
     vb = mlVerbosity s
 
 -- Helper function to exponentiate log domain values with a double value.
@@ -356,7 +415,7 @@
   ML MarginalLikelihood
 sssWrapper s prf lhf cc mn i0 g = do
   logInfoB "Stepping stone sampling."
-  logLhss <- mlRun k bsForward' em vb prf lhf cc mn i0 g
+  logLhss <- mlRunPar pm k (zip [1 ..] bsForward') em vb prf lhf cc mn i0 g
   logInfoB "The last point does not need to be sampled with stepping stone sampling."
   logDebugB "sssWrapper: Calculate marginal likelihood."
   return $ sssCalculateMarginalLikelihood bsForward logLhss
@@ -365,6 +424,7 @@
     bsForward = getPoints k
     bsForward' = init bsForward
     em = mlExecutionMode s
+    pm = mlParallelizationMode s
     vb = mlVerbosity s
 
 -- | Estimate the marginal likelihood.
@@ -382,9 +442,9 @@
   -- Initialize.
   e <- initializeEnvironment s
 
-  when (mlVerbosity s == Debug) $ do
-    let n = fromAnalysisName $ mlAnalysisName s
-    createDirectoryIfMissing True n
+  -- Create marginal likelihood analysis directory.
+  let n = fromAnalysisName $ mlAnalysisName s
+  createDirectoryIfMissing True n
 
   -- Run.
   runReaderT