diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -5,6 +5,32 @@
 ## Unreleased changes
 
 
+## 0.3.0
+
+-   New shorter example/test for dating trees.
+-   `noData` allows running a chain without likelihood function.
+-   Give proposal parameters `PName`, `PDescription`, and `PWeight` newtype
+    wrappers.
+-   Give `Tune` a data type.
+-   Allow periodical cleansing of state (`Cleaner`).
+-   Add description string to proposals, so that they can be identified in an
+    easier way.
+-   Add simplices and proposals on simplices.
+-   `slideUniform` renamed to `slideUniformSymmetric`.
+-   Merge tools into internal.
+-   Do not export internal modules.
+
+
+## 0.2.4
+
+-   **Change order of arguments for proposals**.
+-   'slideStem' was renamed to 'slideBranch'.
+-   Change ProposalSimple from newtype to type.
+-   Contravariant instances of parameter and batch monitors. Use `(>$<)` instead
+    of `(@.)` and `(@#)`.
+-   Add `gammaDirichlet` prior for partitioned dating analyses.
+
+
 ## 0.2.3
 
 -   Contrary proposals.
diff --git a/bench/Normal.hs b/bench/Normal.hs
--- a/bench/Normal.hs
+++ b/bench/Normal.hs
@@ -37,7 +37,7 @@
 proposals :: Cycle Double
 proposals =
   fromList
-    [slideSymmetric 1.0 "medium" 1 True]
+    [slideSymmetric 1.0 (PName "Medium") (PWeight 1) Tune]
 
 mons :: [MonitorParameter Double]
 mons = [monitorDouble "mu"]
@@ -45,9 +45,6 @@
 monStd :: MonitorStdOut Double
 monStd = monitorStdOut mons 200
 
--- monFile :: MonitorFile Double
--- monFile = monitorFile "Mu" mons 200
-
 mon :: Monitor Double
 mon = Monitor monStd [] []
 
@@ -68,7 +65,7 @@
 proposalsBactrian :: Cycle Double
 proposalsBactrian =
   fromList
-    [slideBactrian 0.5 1.0 "bactrian" 1 True]
+    [slideBactrian 0.5 1.0 (PName "Bactrian") (PWeight 1) Tune]
 
 normalBactrianBench :: GenIO -> IO ()
 normalBactrianBench g = do
diff --git a/bench/Poisson.hs b/bench/Poisson.hs
--- a/bench/Poisson.hs
+++ b/bench/Poisson.hs
@@ -39,19 +39,19 @@
     m = sum ys / fromIntegral (length ys)
 
 f :: Int -> Double -> I -> Log Double
-f ft yr (alpha, beta) = Exp $ logProbability (poisson l) (fromIntegral ft)
+f ft yr (a, b) = Exp $ logProbability (poisson l) (fromIntegral ft)
   where
-    l = exp $ alpha + beta * yr
+    l = exp $ a + b * yr
 
 lh :: I -> Log Double
 lh x =
   product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]
 
 proposalAlpha :: Proposal I
-proposalAlpha = _1 @~ slideSymmetric 0.2 "alpha" 2 False
+proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (PWeight 1) NoTune
 
 proposalBeta :: Proposal I
-proposalBeta = _2 @~ slideSymmetric 0.2 "beta" 1 False
+proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (PWeight 1) NoTune
 
 proposals :: Cycle I
 proposals = fromList [proposalAlpha, proposalBeta]
diff --git a/mcmc.cabal b/mcmc.cabal
--- a/mcmc.cabal
+++ b/mcmc.cabal
@@ -1,6 +1,6 @@
 cabal-version:  2.2
 name:           mcmc
-version:        0.2.4
+version:        0.3.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>
 category:       Math, Statistics
@@ -23,7 +23,6 @@
 library
   exposed-modules:
       Mcmc
-      Mcmc.Internal.ByteString
       Mcmc.Item
       Mcmc.Mcmc
       Mcmc.Metropolis
@@ -38,12 +37,14 @@
       Mcmc.Proposal.Generic
       Mcmc.Proposal.Scale
       Mcmc.Proposal.Slide
+      Mcmc.Proposal.Simplex
       Mcmc.Save
       Mcmc.Status
-      Mcmc.Tools.Shuffle
       Mcmc.Trace
       Mcmc.Verbosity
   other-modules:
+      Mcmc.Internal.ByteString
+      Mcmc.Internal.Shuffle
       Paths_mcmc
   autogen-modules:
       Paths_mcmc
@@ -57,10 +58,13 @@
     , containers
     , data-default
     , directory
+    , dirichlet
     , double-conversion
     , log-domain
+    , math-functions
     , microlens
     , mwc-random
+    , primitive
     , statistics
     , time
     , transformers
@@ -77,7 +81,7 @@
       Paths_mcmc
   hs-source-dirs:
       test
-  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N
+  ghc-options: -Wall
   build-depends:
       QuickCheck
     , base >=4.7 && <5
@@ -100,7 +104,7 @@
       Paths_mcmc
   hs-source-dirs:
       bench
-  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N
+  ghc-options: -Wall
   build-depends:
       base >=4.7 && <5
     , criterion
diff --git a/src/Mcmc.hs b/src/Mcmc.hs
--- a/src/Mcmc.hs
+++ b/src/Mcmc.hs
@@ -82,17 +82,21 @@
     -- B. R., Huelsenbeck, J. P., …, Revbayes: bayesian phylogenetic inference using
     -- graphical models and an interactive model-specification language, Systematic
     -- Biology, 65(4), 726–736 (2016). http://dx.doi.org/10.1093/sysbio/syw021
+    PName (..),
+    PWeight (..),
     Proposal,
     (@~),
+    Tune (..),
     scale,
     scaleUnbiased,
     scaleContrarily,
     scaleBactrian,
     slide,
     slideSymmetric,
-    slideUniform,
+    slideUniformSymmetric,
     slideContrarily,
     slideBactrian,
+    module Mcmc.Proposal.Simplex,
     Cycle,
     fromList,
     Order (..),
@@ -113,10 +117,13 @@
     -- space (see above) and to monitor the MCMC run, as well as some auxiliary
     -- information.
     status,
+    Cleaner (..),
+    cleanWith,
     saveWith,
     force,
     quiet,
     debug,
+    noData,
 
     -- * Monitor
 
@@ -160,6 +167,7 @@
 import Mcmc.Proposal
 import Mcmc.Proposal.Bactrian
 import Mcmc.Proposal.Scale
+import Mcmc.Proposal.Simplex
 import Mcmc.Proposal.Slide
 import Mcmc.Save
 import Mcmc.Status
diff --git a/src/Mcmc/Internal/Shuffle.hs b/src/Mcmc/Internal/Shuffle.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Internal/Shuffle.hs
@@ -0,0 +1,64 @@
+-- |
+-- Module      :  Mcmc.Internal.Shuffle
+-- Description :  Shuffle a list
+-- Copyright   :  (c) Dominik Schrempf 2020
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  unstable
+-- Portability :  portable
+--
+-- Creation date: Wed May 20 14:37:09 2020.
+--
+-- From https://wiki.haskell.org/Random_shuffle.
+module Mcmc.Internal.Shuffle
+  ( shuffle,
+    shuffleN,
+    grabble,
+  )
+where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Vector (Vector)
+import qualified Data.Vector as V
+import qualified Data.Vector.Mutable as M
+import System.Random.MWC
+  ( GenIO,
+    uniformR,
+  )
+
+-- | Shuffle a list.
+shuffle :: [a] -> GenIO -> IO [a]
+shuffle xs g = head <$> grabble xs 1 (length xs) g
+
+-- | Shuffle a list @n@ times.
+shuffleN :: [a] -> Int -> GenIO -> IO [[a]]
+shuffleN xs n = grabble xs n (length xs)
+
+-- -- Using System.Random.Shuffle. Speed is the same, so stay without additional dependency.
+-- -- | Shuffle a list @n@ times.
+-- shuffleN :: [a] -> Int -> GenIO -> IO [[a]]
+-- shuffleN xs n g = replicateM n $ fmap (shuffle xs) (rseqM (length xs - 1) g)
+--   where
+--     rseqM :: Int -> GenIO -> IO [Int]
+--     rseqM 0 _ = return []
+--     rseqM i gen = liftM2 (:) (uniformR (0, i) gen) (rseqM (i - 1) gen)
+
+-- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@
+-- elements from @xs@, without replacement, and that @m@ times.
+grabble :: [a] -> Int -> Int -> GenIO -> IO [[a]]
+grabble xs m n gen = do
+  swapss <- replicateM m $
+    forM [0 .. min (l - 1) n] $ \i -> do
+      j <- uniformR (i, l) gen
+      return (i, j)
+  return $ map (V.toList . V.take n . swapElems (V.fromList xs)) swapss
+  where
+    l = length xs - 1
+
+swapElems :: Vector a -> [(Int, Int)] -> Vector a
+swapElems xs swaps = runST $ do
+  mxs <- V.unsafeThaw xs
+  mapM_ (uncurry $ M.unsafeSwap mxs) swaps
+  V.unsafeFreeze mxs
diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs
--- a/src/Mcmc/Mcmc.hs
+++ b/src/Mcmc/Mcmc.hs
@@ -15,15 +15,16 @@
 -- Functions to work with the 'Mcmc' state transformer.
 module Mcmc.Mcmc
   ( Mcmc,
-    mcmcOutT,
+    mcmcOutB,
     mcmcOutS,
-    mcmcWarnT,
+    mcmcWarnB,
     mcmcWarnS,
-    mcmcInfoT,
+    mcmcInfoB,
     mcmcInfoS,
-    mcmcDebugT,
+    mcmcDebugB,
     mcmcDebugS,
     mcmcAutotune,
+    mcmcClean,
     mcmcResetA,
     mcmcSummarizeCycle,
     mcmcReport,
@@ -40,12 +41,14 @@
 import Data.Maybe
 import Data.Time.Clock
 import Data.Time.Format
+import Mcmc.Item
 import Mcmc.Monitor
 import Mcmc.Monitor.Time
 import Mcmc.Proposal
 import Mcmc.Save
 import Mcmc.Status hiding (debug)
 import Mcmc.Verbosity
+import Numeric.Log
 import System.Directory
 import System.IO
 import Prelude hiding (cycle)
@@ -58,66 +61,96 @@
 msgPrepare c t = BL.cons c $ ": " <> t
 
 -- | Write to standard output and log file.
-mcmcOutT :: BL.ByteString -> Mcmc a ()
-mcmcOutT msg = do
+mcmcOutB :: BL.ByteString -> Mcmc a ()
+mcmcOutB msg = do
   h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> gets logHandle
   liftIO $ BL.putStrLn msg >> BL.hPutStrLn h msg
 
 -- | Write to standard output and log file.
 mcmcOutS :: String -> Mcmc a ()
-mcmcOutS = mcmcOutT . BL.pack
+mcmcOutS = mcmcOutB . BL.pack
 
 -- Perform warning action.
 mcmcWarnA :: Mcmc a () -> Mcmc a ()
 mcmcWarnA a = gets verbosity >>= \v -> info v a
 
 -- | Print warning message.
-mcmcWarnT :: BL.ByteString -> Mcmc a ()
-mcmcWarnT = mcmcWarnA . mcmcOutT . msgPrepare 'W'
+mcmcWarnB :: BL.ByteString -> Mcmc a ()
+mcmcWarnB = mcmcWarnA . mcmcOutB . msgPrepare 'W'
 
 -- | Print warning message.
 mcmcWarnS :: String -> Mcmc a ()
-mcmcWarnS = mcmcWarnT . BL.pack
+mcmcWarnS = mcmcWarnB . BL.pack
 
 -- Perform info action.
 mcmcInfoA :: Mcmc a () -> Mcmc a ()
 mcmcInfoA a = gets verbosity >>= \v -> info v a
 
 -- | Print info message.
-mcmcInfoT :: BL.ByteString -> Mcmc a ()
-mcmcInfoT = mcmcInfoA . mcmcOutT . msgPrepare 'I'
+mcmcInfoB :: BL.ByteString -> Mcmc a ()
+mcmcInfoB = mcmcInfoA . mcmcOutB . msgPrepare 'I'
 
 -- | Print info message.
 mcmcInfoS :: String -> Mcmc a ()
-mcmcInfoS = mcmcInfoT . BL.pack
+mcmcInfoS = mcmcInfoB . BL.pack
 
 -- Perform debug action.
 mcmcDebugA :: Mcmc a () -> Mcmc a ()
 mcmcDebugA a = gets verbosity >>= \v -> debug v a
 
 -- | Print debug message.
-mcmcDebugT :: BL.ByteString -> Mcmc a ()
-mcmcDebugT = mcmcDebugA . mcmcOutT . msgPrepare 'D'
+mcmcDebugB :: BL.ByteString -> Mcmc a ()
+mcmcDebugB = mcmcDebugA . mcmcOutB . msgPrepare 'D'
 
 -- | Print debug message.
 mcmcDebugS :: String -> Mcmc a ()
-mcmcDebugS = mcmcDebugT . BL.pack
+mcmcDebugS = mcmcDebugB . BL.pack
 
 -- | Auto tune the 'Proposal's in the 'Cycle' of the chain. Reset acceptance counts.
 -- See 'autotuneCycle'.
 mcmcAutotune :: Mcmc a ()
 mcmcAutotune = do
-  mcmcDebugT "Auto tune."
+  mcmcDebugB "Auto tune."
   s <- get
   let a = acceptance s
       c = cycle s
       c' = autotuneCycle a c
   put $ s {cycle = c'}
 
+-- | Clean the state.
+mcmcClean :: Mcmc a ()
+mcmcClean = do
+  s <- get
+  let cl = cleaner s
+      i = iteration s
+  case cl of
+    Just (Cleaner n f) | i `mod` n == 0 -> do
+      mcmcDebugB "Clean state."
+      let (Item st pr lh) = item s
+      mcmcDebugS $
+        "Old log prior and log likelihood: " ++ show (ln pr) ++ ", " ++ show (ln lh) ++ "."
+      let prF = priorF s
+          lhF = likelihoodF s
+          st' = f st
+          pr' = prF st'
+          lh' = lhF st'
+      mcmcDebugS $
+        "New log prior and log likelihood: " ++ show (ln pr') ++ ", " ++ show (ln lh') ++ "."
+      let dLogPr = abs $ ln pr - ln pr'
+          dLogLh = abs $ ln lh - ln lh'
+      when
+        (dLogPr > 0.01)
+        (mcmcWarnS $ "Log of old and new prior differ by " ++ show dLogPr ++ ".")
+      when
+        (dLogPr > 0.01)
+        (mcmcWarnS $ "Log of old and new likelihood differ by " ++ show dLogLh ++ ".")
+      put $ s {item = Item st' pr' lh'}
+    _ -> return ()
+
 -- | Reset acceptance counts.
 mcmcResetA :: Mcmc a ()
 mcmcResetA = do
-  mcmcDebugT "Reset acceptance ratios."
+  mcmcDebugB "Reset acceptance ratios."
   s <- get
   let a = acceptance s
   put $ s {acceptance = resetA a}
@@ -150,7 +183,7 @@
         (True, _, _) -> openFile lfn AppendMode
   put s {logHandle = mh}
   mcmcDebugS $ "Log file name: " ++ lfn ++ "."
-  mcmcDebugT "Log file opened."
+  mcmcDebugB "Log file opened."
 
 -- Set the total number of iterations, the current time and open the 'Monitor's
 -- of the chain. See 'mOpen'.
@@ -176,14 +209,18 @@
   let b = burnInIterations s
       t = autoTuningPeriod s
       n = iterations s
+      c = cleaner s
   case b of
     Just b' -> mcmcInfoS $ "Burn in for " <> show b' <> " iterations."
     Nothing -> return ()
   case t of
     Just t' -> mcmcInfoS $ "Auto tune every " <> show t' <> " iterations (during burn in only)."
     Nothing -> return ()
+  case c of
+    Just (Cleaner c' _) -> mcmcInfoS $ "Clean state every " <> show c' <> " iterations."
+    Nothing -> return ()
   mcmcInfoS $ "Run chain for " <> show n <> " iterations."
-  mcmcInfoT "Initial state."
+  mcmcInfoB "Initial state."
   mcmcMonitorExec
 
 -- Save the status of an MCMC run. See 'saveStatus'.
@@ -192,11 +229,11 @@
   s <- get
   case save s of
     Just n -> do
-      mcmcInfoT $ "Save Markov chain with trace of length " <> BL.pack (show n) <> "."
-      mcmcInfoT "For long traces, or complex objects, this may take a while."
+      mcmcInfoB $ "Save Markov chain with trace of length " <> BL.pack (show n) <> "."
+      mcmcInfoB "For long traces, or complex objects, this may take a while."
       liftIO $ saveStatus (name s <> ".mcmc") s
-      mcmcInfoT "Done saving Markov chain."
-    Nothing -> mcmcInfoT "Do not save the Markov chain."
+      mcmcInfoB "Done saving Markov chain."
+    Nothing -> mcmcInfoB "Do not save the Markov chain."
 
 -- | Execute the 'Monitor's of the chain. See 'mExec'.
 mcmcMonitorExec :: ToJSON a => Mcmc a ()
@@ -209,14 +246,14 @@
       tr = trace s
       vb = verbosity s
   mt <- liftIO $ mExec vb i ss st tr j m
-  forM_ mt mcmcOutT
+  forM_ mt mcmcOutB
 
 -- Close the 'Monitor's of the chain. See 'mClose'.
 mcmcClose :: ToJSON a => Mcmc a ()
 mcmcClose = do
   s <- get
-  mcmcSummarizeCycle >>= mcmcInfoT
-  mcmcInfoT "Metropolis-Hastings sampler finished."
+  mcmcSummarizeCycle >>= mcmcInfoB
+  mcmcInfoB "Metropolis-Hastings sampler finished."
   let m = monitor s
   m' <- liftIO $ mClose m
   put $ s {monitor = m'}
@@ -225,7 +262,7 @@
   let rt = case start s of
         Nothing -> error "mcmcClose: Start time not set."
         Just (_, st) -> t `diffUTCTime` st
-  mcmcInfoT $ "Wall clock run time: " <> renderDuration rt <> "."
+  mcmcInfoB $ "Wall clock run time: " <> renderDuration rt <> "."
   mcmcInfoS $ "End time: " <> fTime t
   case logHandle s of
     Just h -> liftIO $ hClose h
diff --git a/src/Mcmc/Metropolis.hs b/src/Mcmc/Metropolis.hs
--- a/src/Mcmc/Metropolis.hs
+++ b/src/Mcmc/Metropolis.hs
@@ -46,9 +46,10 @@
 -- Handbook of Markov Chain Monte Carlo), or (b) almost surely reject the
 -- proposal when either fY or q are zero (Chapter 1). Since I trust the author
 -- of Chapter 1 (Charles Geyer) I choose to follow option (b).
-mhRatio :: Log Double -> Log Double -> Log Double -> Log Double
--- q = qYX / qXY
-mhRatio fX fY q = fY * q / fX
+mhRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double
+-- q = qYX / qXY * jXY; see 'ProposalSimple'.
+-- j = Jacobian.
+mhRatio fX fY q j = fY / fX * q * j
 {-# INLINE mhRatio #-}
 
 mhPropose :: Proposal a -> Mcmc a ()
@@ -61,11 +62,11 @@
       a = acceptance s
       g = generator s
   -- 1. Sample new state.
-  (!y, !q) <- liftIO $ p x g
+  (!y, !q, !j) <- liftIO $ p x g
   -- 2. Calculate Metropolis-Hastings ratio.
   let !pY = pF y
       !lY = lF y
-      !r = mhRatio (pX * lX) (pY * lY) q
+      !r = mhRatio (pX * lX) (pY * lY) q j
   -- 3. Accept or reject.
   if ln r >= 0.0
     then put $ s {item = Item y pY lY, acceptance = pushA m True a}
@@ -86,6 +87,7 @@
       t = trace s
       n = iteration s
   put $ s {trace = pushT i t, iteration = succ n}
+  mcmcClean
   mcmcMonitorExec
 
 -- Run N iterations.
@@ -104,13 +106,13 @@
   | b > t = do
     mcmcResetA
     mhNIter t
-    mcmcSummarizeCycle >>= mcmcDebugT
+    mcmcSummarizeCycle >>= mcmcDebugB
     mcmcAutotune
     mhBurnInN (b - t) (Just t)
   | otherwise = do
     mcmcResetA
     mhNIter b
-    mcmcSummarizeCycle >>= mcmcInfoT
+    mcmcSummarizeCycle >>= mcmcInfoB
     mcmcInfoS $ "Acceptance ratios calculated over the last " <> show b <> " iterations."
 mhBurnInN b Nothing = mhNIter b
 
@@ -123,7 +125,7 @@
     mcmcInfoS $ "Burn in for " <> show b <> " cycles."
     mcmcDebugS $ "Auto tuning period is " <> show t <> "."
     mhBurnInN b t
-    mcmcInfoT "Burn in finished."
+    mcmcInfoB "Burn in finished."
 
 -- Run for given number of iterations.
 mhRun :: ToJSON a => Int -> Mcmc a ()
@@ -142,8 +144,8 @@
 
 mhT :: ToJSON a => Mcmc a ()
 mhT = do
-  mcmcInfoT "Metropolis-Hastings sampler."
-  mcmcSummarizeCycle >>= mcmcInfoT
+  mcmcInfoB "Metropolis-Hastings sampler."
+  mcmcSummarizeCycle >>= mcmcInfoB
   mcmcReport
   s <- get
   let b = fromMaybe 0 (burnInIterations s)
@@ -152,9 +154,9 @@
 
 mhContinueT :: ToJSON a => Int -> Mcmc a ()
 mhContinueT dn = do
-  mcmcInfoT "Continuation of Metropolis-Hastings sampler."
+  mcmcInfoB "Continuation of Metropolis-Hastings sampler."
   mcmcInfoS $ "Run chain for " <> show dn <> " additional iterations."
-  mcmcSummarizeCycle >>= mcmcInfoT
+  mcmcSummarizeCycle >>= mcmcInfoB
   mhRun dn
 
 -- | Continue a Markov chain for a given number of Metropolis-Hastings steps.
diff --git a/src/Mcmc/Monitor/Parameter.hs b/src/Mcmc/Monitor/Parameter.hs
--- a/src/Mcmc/Monitor/Parameter.hs
+++ b/src/Mcmc/Monitor/Parameter.hs
@@ -43,7 +43,7 @@
     mpFunc :: a -> BB.Builder
   }
 
-instance Contravariant (MonitorParameter) where
+instance Contravariant MonitorParameter where
   contramap f (MonitorParameter n m) = MonitorParameter n (m . f)
 
 -- | Convert a parameter monitor from one data type to another.
diff --git a/src/Mcmc/Prior.hs b/src/Mcmc/Prior.hs
--- a/src/Mcmc/Prior.hs
+++ b/src/Mcmc/Prior.hs
@@ -13,7 +13,9 @@
 -- Creation date: Thu Jul 23 13:26:14 2020.
 module Mcmc.Prior
   ( -- * Continuous priors
+    largerThan,
     positive,
+    lowerThan,
     negative,
     uniform,
     normal,
@@ -35,17 +37,25 @@
 import qualified Statistics.Distribution.Gamma as S
 import qualified Statistics.Distribution.Normal as S
 
--- | Improper uniform prior; larger than 0.
+-- | Improper uniform prior; strictly larger than a given value.
+largerThan :: Double -> Double -> Log Double
+largerThan a x
+  | x <= a = 0
+  | otherwise = 1
+
+-- | Improper uniform prior; strictly larger than zero.
 positive :: Double -> Log Double
-positive x
-  | x <= 0 = 0
+positive = largerThan 0
+
+-- | Improper uniform prior; strictly lower than a given value.
+lowerThan :: Double -> Double -> Log Double
+lowerThan b x
+  | x >= b = 0
   | otherwise = 1
 
--- | Improper uniform prior; lower than 0.
+-- | Improper uniform prior; strictly lower than zero.
 negative :: Double -> Log Double
-negative x
-  | x >= 0 = 0
-  | otherwise = 1
+negative = lowerThan 0
 
 -- | Uniform prior on [a, b].
 uniform ::
diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs
--- a/src/Mcmc/Proposal.hs
+++ b/src/Mcmc/Proposal.hs
@@ -2,8 +2,6 @@
 {-# LANGUAGE OverloadedStrings #-}
 {-# LANGUAGE RankNTypes #-}
 
--- TODO: Proposals on simplices: SimplexElementScale (?).
-
 -- |
 -- Module      :  Mcmc.Proposal
 -- Description :  Proposals and cycles
@@ -17,10 +15,14 @@
 -- Creation date: Wed May 20 13:42:53 2020.
 module Mcmc.Proposal
   ( -- * Proposal
+    PName (..),
+    PDescription (..),
+    PWeight (..),
     Proposal (..),
     (@~),
     ProposalSimple,
     Tuner (tParam, tFunc),
+    Tune (..),
     createProposal,
     tune,
 
@@ -57,36 +59,54 @@
 import Data.Maybe
 import Lens.Micro
 import Mcmc.Internal.ByteString
-import Mcmc.Tools.Shuffle
+import Mcmc.Internal.Shuffle
 import Numeric.Log hiding (sum)
 import System.Random.MWC
 
+-- | Proposal name.
+newtype PName = PName {fromPName :: String}
+  deriving (Show, Eq, Ord)
+
+-- | Proposal description.
+newtype PDescription = PDescription {fromPDescription :: String}
+  deriving (Show, Eq, Ord)
+
+-- | The weight determines how often a 'Proposal' is executed per iteration of
+-- the Markov chain.
+newtype PWeight = PWeight {fromPWeight :: Int}
+  deriving (Show, Eq, Ord)
+
 -- | A 'Proposal' is an instruction about how the Markov chain will traverse the
 -- state space @a@. Essentially, it is a probability mass or probability density
 -- conditioned on the current state (i.e., a kernel).
 --
 -- A 'Proposal' may be tuneable in that it contains information about how to enlarge
 -- or shrink the step size to tune the acceptance ratio.
+--
+-- No proposals with the same name and description are allowed in a 'Cycle'.
 data Proposal a = Proposal
-  { -- | Name (no proposals with the same name are allowed in a 'Cycle').
-    pName :: String,
+  { -- | Name of the affected variable.
+    pName :: PName,
+    -- | Description of the proposal type and parameters.
+    pDescription :: PDescription,
     -- | The weight determines how often a 'Proposal' is executed per iteration of
     -- the Markov chain.
-    pWeight :: Int,
+    pWeight :: PWeight,
     -- | Simple proposal without name, weight, and tuning information.
     pSimple :: ProposalSimple a,
     -- | Tuning is disabled if set to 'Nothing'.
     pTuner :: Maybe (Tuner a)
   }
 
+-- XXX: This should be removed.
 instance Show (Proposal a) where
-  show m = show $ pName m
+  show m = fromPName (pName m) <> " " <> fromPDescription (pDescription m) <> ", weight " <> show (fromPWeight $ pWeight m)
 
 instance Eq (Proposal a) where
-  m == n = pName m == pName n
+  m == n = pName m == pName n && pDescription m == pDescription n
 
 instance Ord (Proposal a) where
-  compare = compare `on` pName
+  compare = compare `on` (\p -> (pDescription p, pName p, pWeight p))
 
 -- | Convert a proposal from one data type to another using a lens.
 --
@@ -96,7 +116,7 @@
 -- scaleFirstEntryOfTuple = _1 @~ scale
 -- @
 (@~) :: Lens' b a -> Proposal a -> Proposal b
-(@~) l (Proposal n w s t) = Proposal n w (convertS l s) (convertT l <$> t)
+(@~) l (Proposal n d w s t) = Proposal n d w (convertS l s) (convertT l <$> t)
 
 -- | Simple proposal without tuning information.
 --
@@ -105,16 +125,26 @@
 --
 -- In order to calculate the Metropolis-Hastings ratio, we need to know the
 -- ratio of the backward to forward kernels (i.e., the probability masses or
--- probability densities). For unbiased proposals, this ratio is 1.0. For biased
--- proposals, the ratio is ?.
-type ProposalSimple a = a -> GenIO -> IO (a, Log Double)
+-- probability densities) and the absolute value of the determinant of the
+-- Jacobian matrix.
+--
+-- For unbiased proposals, these values are 1.0 such that
+--
+-- @
+-- proposalSimpleUnbiased x g = return (x', 1.0, 1.0)
+-- @
+--
+-- For biased proposals, the kernel ratio is qYX / qXY, where qXY is the
+-- probability density to move from X to Y, and the absolute value of the
+-- determinant of the Jacobian matrix differs from 1.0.
+type ProposalSimple a = a -> GenIO -> IO (a, Log Double, Log Double)
 
 convertS :: Lens' b a -> ProposalSimple a -> ProposalSimple b
 convertS l s = s'
   where
     s' v g = do
-      (x', r) <- s (v ^. l) g
-      return (set l x' v, r)
+      (x', r, j) <- s (v ^. l) g
+      return (set l x' v, r, j)
 
 -- | Tune the acceptance ratio of a 'Proposal'; see 'tune', or 'autotuneCycle'.
 data Tuner a = Tuner
@@ -127,21 +157,27 @@
   where
     f' x = convertS l $ f x
 
+-- | Tune the proposal?
+data Tune = Tune | NoTune
+  deriving (Show, Eq)
+
 -- | Create a possibly tuneable proposal.
 createProposal ::
+  -- | Description of the proposal type and parameters.
+  PDescription ->
   -- | Function creating a simple proposal for a given tuning parameter. The
   -- larger the tuning parameter, the larger the proposal (and the lower the
   -- expected acceptance ratio), and vice versa.
   (Double -> ProposalSimple a) ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Activate tuning?
-  Bool ->
+  Tune ->
   Proposal a
-createProposal f n w True = Proposal n w (f 1.0) (Just $ Tuner 1.0 f)
-createProposal f n w False = Proposal n w (f 1.0) Nothing
+createProposal d f n w Tune = Proposal n d w (f 1.0) (Just $ Tuner 1.0 f)
+createProposal d f n w NoTune = Proposal n d w (f 1.0) Nothing
 
 -- Minimal tuning parameter; subject to change.
 tuningParamMin :: Double
@@ -223,11 +259,9 @@
 fromList [] =
   error "fromList: Received an empty list but cannot create an empty Cycle."
 fromList xs =
-  if length (nub nms) == length nms
+  if length (nub xs) == length xs
     then Cycle xs def
-    else error "fromList: Proposals don't have unique names."
-  where
-    nms = map pName xs
+    else error "fromList: Proposals are not unique."
 
 -- | Set the order of 'Proposal's in a 'Cycle'.
 setOrder :: Order -> Cycle a -> Cycle a
@@ -243,7 +277,7 @@
     return [psR ++ reverse psR | psR <- psRs]
   SequentialReversibleO -> return $ replicate n $ ps ++ reverse ps
   where
-    !ps = concat [replicate (pWeight m) m | m <- xs]
+    !ps = concat [replicate (fromPWeight $ pWeight m) m | m <- xs]
 
 -- | Tune 'Proposal's in the 'Cycle'. See 'tune'.
 tuneCycle :: Map (Proposal a) Double -> Cycle a -> Cycle a
@@ -271,10 +305,12 @@
   BL.ByteString ->
   BL.ByteString ->
   BL.ByteString ->
+  BL.ByteString ->
   BL.ByteString
-renderRow name weight nAccept nReject acceptRatio tuneParam manualAdjustment = "   " <> nm <> wt <> na <> nr <> ra <> tp <> mt
+renderRow name ptype weight nAccept nReject acceptRatio tuneParam manualAdjustment = "   " <> nm <> pt <> wt <> na <> nr <> ra <> tp <> mt
   where
     nm = alignLeft 30 name
+    pt = alignLeft 50 ptype
     wt = alignRight 8 weight
     na = alignRight 15 nAccept
     nr = alignRight 15 nReject
@@ -284,12 +320,13 @@
 
 proposalHeader :: BL.ByteString
 proposalHeader =
-  renderRow "Proposal" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter" "Consider manual adjustment"
+  renderRow "Name" "Description" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter" "Consider manual adjustment"
 
 summarizeProposal :: Proposal a -> Maybe (Int, Int, Double) -> BL.ByteString
 summarizeProposal m r =
   renderRow
-    (BL.pack name)
+    (BL.pack $ fromPName $ pName m)
+    (BL.pack $ fromPDescription $ pDescription m)
     weight
     nAccept
     nReject
@@ -297,8 +334,7 @@
     tuneParamStr
     manualAdjustmentStr
   where
-    name = pName m
-    weight = BB.toLazyByteString $ BB.intDec $ pWeight m
+    weight = BB.toLazyByteString $ BB.intDec $ fromPWeight $ pWeight m
     nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) r
     nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) r
     acceptRatio = BL.fromStrict $ maybe "" (BC.toFixed 3 . (^. _3)) r
@@ -324,7 +360,7 @@
       ++ [hLine proposalHeader]
   where
     ps = ccProposals c
-    mpi = BB.toLazyByteString $ BB.intDec $ sum $ map pWeight ps
+    mpi = BB.toLazyByteString $ BB.intDec $ sum $ map (fromPWeight . pWeight) ps
     ar m = acceptanceRatio m a
 
 -- | For each key @k@, store the number of accepted and rejected proposals.
diff --git a/src/Mcmc/Proposal/Bactrian.hs b/src/Mcmc/Proposal/Bactrian.hs
--- a/src/Mcmc/Proposal/Bactrian.hs
+++ b/src/Mcmc/Proposal/Bactrian.hs
@@ -52,10 +52,10 @@
   Double ->
   Double ->
   GenIO ->
-  IO (Double, Log Double)
+  IO (Double, Log Double, Log Double)
 bactrianAdditive m s x g = do
   dx <- genBactrian m s g
-  return (x + dx, 1.0)
+  return (x + dx, 1.0, 1.0)
 
 -- bactrianSimple lens spike stdDev tune forwardOp backwardOp
 bactrianAdditiveSimple ::
@@ -75,9 +75,9 @@
 -- The [Bactrian
 -- kernel](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/figure/fig01) is
 -- a mixture of two symmetrically arranged normal distributions. The spike
--- parameter loosely determines the standard deviations of the individual humps
--- while the second parameter refers to the standard deviation of the complete
--- Bactrian kernel.
+-- parameter (0, 1) loosely determines the standard deviations (>0.0) of the
+-- individual humps while the second parameter refers to the standard deviation
+-- of the complete Bactrian kernel.
 --
 -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/.
 slideBactrian ::
@@ -86,13 +86,15 @@
   -- | Standard deviation.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-slideBactrian m s = createProposal (bactrianAdditiveSimple m s)
+slideBactrian m s = createProposal description (bactrianAdditiveSimple m s)
+  where
+    description = PDescription $ "Slide Bactrian; spike: " ++ show m ++ ", sd: " ++ show s
 
 -- We have:
 -- x  (1+dx ) = x'
@@ -108,12 +110,14 @@
   Double ->
   Double ->
   GenIO ->
-  IO (Double, Log Double)
+  IO (Double, Log Double, Log Double)
 bactrianMult m s x g = do
-  dx <- genBactrian m s g
-  let qXY = logDensityBactrian m s dx
-      qYX = logDensityBactrian m s (fInv dx)
-  return (x * (1 + dx), qYX / qXY)
+  du <- genBactrian m s g
+  let qXY = logDensityBactrian m s du
+      qYX = logDensityBactrian m s (fInv du)
+      u = 1.0 + du
+      jac = Exp $ log $ recip u
+  return (x * u, qYX / qXY, jac)
 
 bactrianMultSimple :: Double -> Double -> Double -> ProposalSimple Double
 bactrianMultSimple m s t
@@ -130,10 +134,12 @@
   -- | Standard deviation.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-scaleBactrian m s = createProposal (bactrianMultSimple m s)
+scaleBactrian m s = createProposal description (bactrianMultSimple m s)
+  where
+    description = PDescription $ "Scale Bactrian; spike: " ++ show m <> ", sd: " <> show s
diff --git a/src/Mcmc/Proposal/Generic.hs b/src/Mcmc/Proposal/Generic.hs
--- a/src/Mcmc/Proposal/Generic.hs
+++ b/src/Mcmc/Proposal/Generic.hs
@@ -18,58 +18,54 @@
 import Mcmc.Proposal
 import Numeric.Log
 import Statistics.Distribution
-import System.Random.MWC
 
-sampleCont ::
-  (ContDistr d, ContGen d) =>
-  d ->
-  (a -> Double -> a) ->
-  Maybe (Double -> Double) ->
-  a ->
-  GenIO ->
-  IO (a, Log Double)
-sampleCont d f mfInv x g = do
-  dx <- genContVar d g
-  let r = case mfInv of
-        Nothing -> 1.0
-        Just fInv ->
-          let qXY = Exp $ logDensity d dx
-              qYX = Exp $ logDensity d (fInv dx)
-           in qYX / qXY
-  return (x `f` dx, r)
-{-# INLINEABLE sampleCont #-}
-
 -- | Generic function to create proposals for continuous parameters ('Double').
 genericContinuous ::
   (ContDistr d, ContGen d) =>
   -- | Probability distribution
   d ->
-  -- | Forward operator, e.g. (+), so that x + dx = x'.
+  -- | Forward operator.
+  --
+  -- For example, for a multiplicative proposal on one variable the forward
+  -- operator is @(*)@, so that @x * u = y@.
   (a -> Double -> a) ->
-  -- | Inverse operator, e.g., 'negate', so that x' + (negate dx) = x. Only
-  -- required for biased proposals.
+  -- | Inverse operator.
+  --
+  -- For example, 'recip' for a multiplicative proposal on one variable, since
+  -- @y * (recip u) = x * u * (recip u) = x@.
+  --
+  -- Required for biased proposals.
   Maybe (Double -> Double) ->
+  -- | Function to compute the absolute value of the determinant of the Jacobian
+  -- matrix. For example, for a multiplicative proposal on one variable, we have
+  --
+  -- @
+  -- detJacobian _ u = Exp $ log $ recip u
+  -- @
+  --
+  -- That is, the determinant of the Jacobian matrix of multiplication is just
+  -- the reciprocal value of @u@ (with conversion to log domain).
+  --
+  -- Required for proposals for which absolute value of the determinant of the
+  -- Jacobian differs from 1.0.
+  --
+  -- Conversion to log domain is necessary, because some determinants of
+  -- Jacobians are very small (or large).
+  Maybe (a -> Double -> Log Double) ->
   ProposalSimple a
-genericContinuous d f fInv = sampleCont d f fInv
-
-sampleDiscrete ::
-  (DiscreteDistr d, DiscreteGen d) =>
-  d ->
-  (a -> Int -> a) ->
-  Maybe (Int -> Int) ->
-  a ->
-  GenIO ->
-  IO (a, Log Double)
-sampleDiscrete d f mfInv x g = do
-  dx <- genDiscreteVar d g
-  let r = case mfInv of
+genericContinuous d f mInv mJac x g = do
+  u <- genContVar d g
+  let r = case mInv of
         Nothing -> 1.0
         Just fInv ->
-          let qXY = Exp $ logProbability d dx
-              qYX = Exp $ logProbability d (fInv dx)
+          let qXY = Exp $ logDensity d u
+              qYX = Exp $ logDensity d (fInv u)
            in qYX / qXY
-  return (x `f` dx, r)
-{-# INLINEABLE sampleDiscrete #-}
+      j = case mJac of
+        Nothing -> 1.0
+        Just fJac -> fJac x u
+  return (x `f` u, r, j)
+{-# INLINEABLE genericContinuous #-}
 
 -- | Generic function to create proposals for discrete parameters ('Int').
 genericDiscrete ::
@@ -82,4 +78,13 @@
   -- required for biased proposals.
   Maybe (Int -> Int) ->
   ProposalSimple a
-genericDiscrete fd f fInv = sampleDiscrete fd f fInv
+genericDiscrete d f mfInv x g = do
+  u <- genDiscreteVar d g
+  let r = case mfInv of
+        Nothing -> 1.0
+        Just fInv ->
+          let qXY = Exp $ logProbability d u
+              qYX = Exp $ logProbability d (fInv u)
+           in qYX / qXY
+  return (x `f` u, r, 1.0)
+{-# INLINEABLE genericDiscrete #-}
diff --git a/src/Mcmc/Proposal/Scale.hs b/src/Mcmc/Proposal/Scale.hs
--- a/src/Mcmc/Proposal/Scale.hs
+++ b/src/Mcmc/Proposal/Scale.hs
@@ -20,12 +20,20 @@
 
 import Mcmc.Proposal
 import Mcmc.Proposal.Generic
+import Numeric.Log
 import Statistics.Distribution.Gamma
 
 -- The actual proposal with tuning parameter. The tuning parameter does not
 -- change the mean.
 scaleSimple :: Double -> Double -> Double -> ProposalSimple Double
-scaleSimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) (*) (Just recip)
+scaleSimple k th t =
+  genericContinuous
+    (gammaDistr (k / t) (th * t))
+    (*)
+    (Just recip)
+    (Just jac)
+  where
+    jac _ = Exp . log . recip
 
 -- | Multiplicative proposal with Gamma distributed kernel.
 scale ::
@@ -34,13 +42,15 @@
   -- | Scale.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-scale k th = createProposal (scaleSimple k th)
+scale k th = createProposal description (scaleSimple k th)
+  where
+    description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th
 
 -- | Multiplicative proposal with Gamma distributed kernel.
 --
@@ -50,20 +60,31 @@
   -- | Shape.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-scaleUnbiased k = createProposal (scaleSimple k (1 / k))
-
-contra :: (Double, Double) -> Double -> (Double, Double)
-contra (x, y) z = (x * z, y / z)
+scaleUnbiased k = createProposal description (scaleSimple k (1 / k))
+  where
+    description = PDescription $ "Scale unbiased; shape: " ++ show k
 
 scaleContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)
-scaleContrarilySimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) contra (Just recip)
+scaleContrarilySimple k th t =
+  genericContinuous
+    (gammaDistr (k / t) (th * t))
+    contra
+    (Just recip)
+    (Just jac)
+  where
+    contra (x, y) u = (x * u, y / u)
+    jac _ u = Exp $ log $ recip $ u * u
 
+-- -- Determinant of Jacobian matrix.
+-- contraJac :: (Double, Double) -> Double
+-- contraJac (x, y) = x * y
+
 -- | Multiplicative proposal with Gamma distributed kernel.
 --
 -- The two values are scaled contrarily so that their product stays constant.
@@ -74,10 +95,12 @@
   -- | Scale.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal (Double, Double)
-scaleContrarily k th = createProposal (scaleContrarilySimple k th)
+scaleContrarily k th = createProposal description (scaleContrarilySimple k th)
+  where
+    description = PDescription $ "Scale contrariliy; shape: " ++ show k ++ ", scale: " ++ show th
diff --git a/src/Mcmc/Proposal/Simplex.hs b/src/Mcmc/Proposal/Simplex.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Proposal/Simplex.hs
@@ -0,0 +1,203 @@
+{-# LANGUAGE TemplateHaskell #-}
+
+-- |
+-- Module      :  Mcmc.Proposal.Simplex
+-- Description :  Proposals on simplices
+-- Copyright   :  (c) Dominik Schrempf, 2020
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  unstable
+-- Portability :  portable
+--
+-- Creation date: Mon Oct 19 15:32:31 2020.
+module Mcmc.Proposal.Simplex
+  ( -- * Elements of simplices
+    Simplex (toVector),
+    simplexUniform,
+    simplexFromVector,
+
+    -- * Proposals on simplices
+    dirichlet,
+    beta,
+  )
+where
+
+import Data.Aeson
+import Data.Aeson.TH
+import qualified Data.Vector.Unboxed as V
+import Mcmc.Proposal
+import Numeric.Log
+import Statistics.Distribution
+import Statistics.Distribution.Beta
+import Statistics.Distribution.Dirichlet
+
+-- import Debug.Trace
+
+-- | An element of a simplex.
+--
+-- A vector of non-negative values summing to one.
+--
+-- The nomenclature is not very consistent, because a K-dimensional simplex is
+-- usually considered to be the set containing all @K@-dimensional vectors with
+-- non-negative elements that sum to 1.0. However, I couldn't come up with a
+-- better name. Maybe @SimplexElement@, but that was too long.
+newtype Simplex = SimplexUnsafe {toVector :: V.Vector Double}
+  deriving (Eq, Show)
+
+$(deriveJSON defaultOptions ''Simplex)
+
+-- Tolerance.
+eps :: Double
+eps = 1e-14
+
+-- Check if vector is normalized with tolerance 'eps'.
+isNormalized :: V.Vector Double -> Bool
+isNormalized v
+  | abs (V.sum v - 1.0) > eps = False
+  | otherwise = True
+
+-- Check if vector contains negative elements.
+isNegative :: V.Vector Double -> Bool
+isNegative = V.any (< 0)
+
+-- | Ensure that the value vector is an element of a simplex.
+--
+-- Return 'Left' if:
+-- - The value vector is empty.
+-- - The value vector contains negative elements.
+-- - The value vector is not normalized.
+simplexFromVector :: V.Vector Double -> Either String Simplex
+simplexFromVector v
+  | V.null v = Left "simplexFromVector: Vector is empty."
+  | isNegative v = Left "simplexFromVector: Vector contains negative elements."
+  | not (isNormalized v) = Left "simplexFromVector: Vector is not normalized."
+  | otherwise = Right $ SimplexUnsafe v
+
+-- | Create the uniform element of the K-dimensional simplex.
+--
+-- Set all values to \(1/D\).
+simplexUniform :: Int -> Simplex
+simplexUniform k = either error id $ simplexFromVector $ V.replicate k (1.0 / fromIntegral k)
+
+-- Tuning function is inverted (high alpha means small steps).
+getTuningFunction :: Double -> (Double -> Double)
+getTuningFunction t = (/ t'')
+  where
+    -- Start with small steps.
+    t' = t / 10000
+    -- Extremely small tuning parameters lead to numeric overflow. The square
+    -- root pulls the tuning parameter closer to 1.0. However, overflow may
+    -- still occur (the involved Gamma functions grow faster than the
+    -- exponential). I did not observe numeric underflow in my tests.
+    t'' = sqrt t'
+
+-- The tuning parameter is proportional to the inverted mean of the shape
+-- parameter values.
+--
+-- The values determining the proposal size have been set using an example
+-- analysis. They are good values for this analysis, but may fail for other
+-- analyses.
+dirichletSimple :: Double -> ProposalSimple Simplex
+dirichletSimple t (SimplexUnsafe xs) g = do
+  -- If @t@ is high and above 1.0, the parameter vector will be low, and the
+  -- variance will be high. If @t@ is low and below 1.0, the parameter vector
+  -- will be high, and the Dirichlet distribution will be very concentrated with
+  -- low variance.
+  let ddXs = either error id $ dirichletDistribution $ V.map tf xs
+  -- traceShowM $ V.map tf xs
+  ys <- dirichletSample ddXs g
+  -- traceShowM ys
+  -- Have to check if parameters are valid (because zeroes do occur).
+  let eitherDdYs = dirichletDistribution $ V.map tf ys
+  let r = case eitherDdYs of
+        -- Set ratio to 0; so that the proposal will not be accepted.
+        Left _ -> 0
+        Right ddYs -> dirichletDensity ddYs xs / dirichletDensity ddXs ys
+  -- I do not think a Jacobian is necessary in this case. I do know that if a
+  -- subset of states is updated a Jacobian would be necessary.
+  --
+  -- traceShowM mhRatio
+  return (SimplexUnsafe ys, r, 1.0)
+  where
+    tf = getTuningFunction t
+
+-- | Dirichlet proposal on a simplex.
+--
+-- For a given element of a K-dimensional simplex, propose a new element of the
+-- K-dimensional simplex. The new element is sampled from the multivariate
+-- Dirichlet distribution with parameter vector being the old element of the
+-- simplex.
+--
+-- The tuning parameter is used to determine the concentration of the Dirichlet
+-- distribution: the lower the tuning parameter, the higher the concentration.
+--
+-- This proposal may have low acceptance ratios. In this case, please see the
+-- coordinate wise 'beta' proposal.
+dirichlet :: PName -> PWeight -> Tune -> Proposal Simplex
+dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple
+
+-- The tuning parameter is the inverted mean of the shape values.
+--
+-- The values determining the proposal size have been set using an example
+-- analysis. They are good values for this analysis, but may fail for other
+-- analyses.
+--
+-- See also the 'dirichlet' proposal.
+betaSimple :: Int -> Double -> ProposalSimple Simplex
+betaSimple i t (SimplexUnsafe xs) g = do
+  -- Shape parameters of beta distribution. Do not assume that the sum of the
+  -- elements of 'xs' is 1.0, because then repeated proposals let the sum of the
+  -- vector diverge.
+  let aX = xI
+      bX = xsSum - xI
+      bdXI = betaDistr (tf aX) (tf bX)
+  -- New value of element i.
+  yI <- genContVar bdXI g
+  -- Shape parameters of beta distribution.
+  let aY = yI
+      bY = 1.0 - yI
+      eitherBdYI = betaDistrE (tf aY) (tf bY)
+  -- See 'dirichlet', which has the same construct.
+  let r = case eitherBdYI of
+        Nothing -> 0
+        Just bdYI -> Exp $ logDensity bdYI xI - logDensity bdXI yI
+      -- The absolute value of the determinant of the Jacobian. Derivation takes
+      -- a while...
+      ja1 = bY / bX
+      jac = Exp $ fromIntegral (V.length xs - 2) * log ja1
+  -- Construct new vector.
+  let -- Normalization function for other elements.
+      -- nf x = x * bY / bX
+      --
+      -- It turns out, that this factor is also needed to compute the determinant
+      -- of the Jacobian above.
+      nf x = x * ja1
+      ys = V.generate (V.length xs) (\j -> if i == j then yI else nf (xs V.! j))
+  return (either error id $ simplexFromVector ys, r, jac)
+  where
+    xI = xs V.! i
+    xsSum = V.sum xs
+    tf = getTuningFunction t
+
+-- | Beta proposal on a specific coordinate @i@ on a simplex.
+--
+-- For a given element of a K-dimensional simplex, propose a new element of the
+-- K-dimensional simplex. The coordinate @i@ of the new element is sampled from
+-- the beta distribution. The other coordinates are normalized such that the
+-- values sum to 1.0. The parameters of the beta distribution are chosen such
+-- that the expected value of the beta distribution is the value of the old
+-- coordinate.
+--
+-- The tuning parameter is used to determine the concentration of the beta
+-- distribution: the lower the tuning parameter, the higher the concentration.
+--
+-- See also the 'dirichlet' proposal.
+--
+-- No "out of bounds" checks are performed during compile time. Run time errors
+-- can occur if @i@ is negative, or if @i-1@ is larger than the length of the
+-- element vector of the simplex.
+beta :: Int -> PName -> PWeight -> Tune -> Proposal Simplex
+beta i = createProposal description (betaSimple i)
+  where
+    description = PDescription $ "Beta; coordinate: " ++ show i
diff --git a/src/Mcmc/Proposal/Slide.hs b/src/Mcmc/Proposal/Slide.hs
--- a/src/Mcmc/Proposal/Slide.hs
+++ b/src/Mcmc/Proposal/Slide.hs
@@ -14,7 +14,7 @@
 module Mcmc.Proposal.Slide
   ( slide,
     slideSymmetric,
-    slideUniform,
+    slideUniformSymmetric,
     slideContrarily,
   )
 where
@@ -26,7 +26,8 @@
 
 -- The actual proposal with tuning parameter.
 slideSimple :: Double -> Double -> Double -> ProposalSimple Double
-slideSimple m s t = genericContinuous (normalDistr m (s * t)) (+) (Just negate)
+slideSimple m s t =
+  genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing
 
 -- | Additive proposal with normally distributed kernel.
 slide ::
@@ -35,17 +36,20 @@
   -- | Standard deviation.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-slide m s = createProposal (slideSimple m s)
+slide m s = createProposal description (slideSimple m s)
+  where
+    description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s
 
 -- The actual proposal with tuning parameter.
 slideSymmetricSimple :: Double -> Double -> ProposalSimple Double
-slideSymmetricSimple s t = genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing
+slideSymmetricSimple s t =
+  genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing
 
 -- | Additive proposal with normally distributed kernel with mean zero. This
 -- proposal is very fast, because the Metropolis-Hastings ratio does not include
@@ -54,39 +58,44 @@
   -- | Standard deviation.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-slideSymmetric s = createProposal (slideSymmetricSimple s)
+slideSymmetric s = createProposal description (slideSymmetricSimple s)
+  where
+    description = PDescription $ "Slide symmetric; sd: " ++ show s
 
 -- The actual proposal with tuning parameter.
 slideUniformSimple :: Double -> Double -> ProposalSimple Double
 slideUniformSimple d t =
-  genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing
+  genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing
 
--- | Additive proposal with uniformly distributed kernel. This proposal is very fast,
--- because the Metropolis-Hastings ratio does not include calculation of the
--- forwards and backwards kernels.
-slideUniform ::
+-- | Additive proposal with uniformly distributed kernel with mean zero. This
+-- proposal is very fast, because the Metropolis-Hastings ratio does not include
+-- calculation of the forwards and backwards kernels.
+slideUniformSymmetric ::
   -- | Delta.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal Double
-slideUniform d = createProposal (slideUniformSimple d)
+slideUniformSymmetric d = createProposal description (slideUniformSimple d)
+  where
+    description = PDescription $ "Slide uniform symmetric; delta: " ++ show d
 
 contra :: (Double, Double) -> Double -> (Double, Double)
-contra (x, y) d = (x + d, y - d)
+contra (x, y) u = (x + u, y - u)
 
 slideContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)
-slideContrarilySimple m s t = genericContinuous (normalDistr m (s * t)) contra (Just negate)
+slideContrarilySimple m s t =
+  genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing
 
 -- | Additive proposal with normally distributed kernel.
 --
@@ -98,10 +107,12 @@
   -- | Standard deviation.
   Double ->
   -- | Name.
-  String ->
-  -- | Weight.
-  Int ->
+  PName ->
+  -- | PWeight.
+  PWeight ->
   -- | Enable tuning.
-  Bool ->
+  Tune ->
   Proposal (Double, Double)
-slideContrarily m s = createProposal (slideContrarilySimple m s)
+slideContrarily m s = createProposal description (slideContrarilySimple m s)
+  where
+    description = PDescription $ "Slide contrarily; mean: " ++ show m ++ ", sd: " ++ show s
diff --git a/src/Mcmc/Save.hs b/src/Mcmc/Save.hs
--- a/src/Mcmc/Save.hs
+++ b/src/Mcmc/Save.hs
@@ -68,7 +68,7 @@
 $(deriveJSON defaultOptions ''Save)
 
 toSave :: Status a -> Save a
-toSave (Status nm it i tr ac br at is f sv vb g _ _ _ _ c _) =
+toSave (Status nm it i tr ac br at is f sv vb g _ _ _ _ _ c _) =
   Save
     nm
     it
@@ -104,9 +104,10 @@
   (a -> Log Double) ->
   Cycle a ->
   Monitor a ->
+  Maybe (Cleaner a) ->
   Save a ->
   Status a
-fromSave p l c m (Save nm it i tr ac' br at is f sv vb g' ts) =
+fromSave pr lh cc m cl (Save nm it i tr ac' br at is f sv vb g' ts) =
   Status
     nm
     it
@@ -122,16 +123,17 @@
     g
     Nothing
     Nothing
-    p
-    l
-    c'
+    pr
+    lh
+    cl
+    cc'
     m
   where
-    ac = transformKeysA [0 ..] (ccProposals c) ac'
+    ac = transformKeysA [0 ..] (ccProposals cc) ac'
     -- TODO: Splitmix. Remove as soon as split mix is used and is available with
     -- the statistics package.
     g = unsafePerformIO $ restore $ toSeed g'
-    c' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals c) ts) c
+    cc' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals cc) ts) cc
 
 -- | Load a 'Status' from file.
 --
@@ -139,23 +141,29 @@
 -- a chain is restored:
 -- - prior function
 -- - likelihood function
+-- - cleaning function
 -- - cycle
 -- - monitor
 --
 -- To avoid incomplete continued runs, the @.mcmc@ file is removed after load.
 loadStatus ::
   FromJSON a =>
+  -- | Prior function.
   (a -> Log Double) ->
+  -- | Likelihood function.
   (a -> Log Double) ->
   Cycle a ->
   Monitor a ->
+  -- | Cleaner, if needed.
+  Maybe (Cleaner a) ->
+  -- | Path of status to load.
   FilePath ->
   IO (Status a)
-loadStatus p l c m fn = do
+loadStatus pr lh cc mn cl fn = do
   res <- eitherDecode . decompress <$> BL.readFile fn
   let s = case res of
         Left err -> error err
-        Right sv -> fromSave p l c m sv
+        Right sv -> fromSave pr lh cc mn cl sv
   -- Check if prior and likelihood matches.
   let Item x svp svl = item s
   -- Recompute and check the prior and likelihood for the last state because the
@@ -163,10 +171,10 @@
   -- function, but having the same prior and likelihood at the last state is
   -- already a good indicator.
   when
-    (p x /= svp)
+    (pr x /= svp)
     (error "loadStatus: Provided prior function does not match the saved prior.")
   when
-    (l x /= svl)
+    (lh x /= svl)
     (error "loadStatus: Provided likelihood function does not match the saved likelihood.")
   removeFile fn
   return s
diff --git a/src/Mcmc/Status.hs b/src/Mcmc/Status.hs
--- a/src/Mcmc/Status.hs
+++ b/src/Mcmc/Status.hs
@@ -21,12 +21,15 @@
 --
 -- Creation date: Tue May  5 18:01:15 2020.
 module Mcmc.Status
-  ( Status (..),
+  ( Cleaner (..),
+    Status (..),
     status,
+    cleanWith,
     saveWith,
     force,
     quiet,
     debug,
+    noData,
   )
 where
 
@@ -42,6 +45,28 @@
 import System.Random.MWC hiding (save)
 import Prelude hiding (cycle)
 
+-- | Clean the state periodically.
+--
+-- The prior and the likelihood will be updated after the cleaning process.
+--
+-- For long chains, successive numerical errors can accumulate such that the
+-- state diverges from honoring specific required constraints. In these cases, a
+-- 'Cleaner' can be used to ensure that the required constraints of the state
+-- are honored. For example, the branches of an ultrametric phylogeny may
+-- diverge slightly after successful many proposals such that the phylogeny is
+-- not anymore ultrametric.
+--
+-- Please be aware that the Markov chain will not converge to the true posterior
+-- distribution if the state is changed substantially! Only apply subtle changes
+-- that are absolutely necessary to preserve the required properties of the
+-- state such as specific numerical constraints.
+data Cleaner a = Cleaner
+  { -- | Clean every given number of iterations.
+    clEvery :: Int,
+    -- | Cleaning function. Executed before monitoring the state.
+    clFunction :: a -> a
+  }
+
 -- | The 'Status' contains all information to run an MCMC chain. It is
 -- constructed using the function 'status'.
 --
@@ -101,6 +126,8 @@
     -- | The likelihood function. The un-normalized posterior is the product of
     -- the prior and the likelihood.
     likelihoodF :: a -> Log Double,
+    -- | Clean the state periodically.
+    cleaner :: Maybe (Cleaner a),
     --
     -- Variables related to the algorithm; not saved.
 
@@ -118,8 +145,8 @@
   (a -> Log Double) ->
   -- | The likelihood function.
   (a -> Log Double) ->
-  -- | A list of 'Proposal's executed in forward order. The
-  -- chain will be logged after each cycle.
+  -- | A list of 'Proposal's executed in forward order. The chain will be logged
+  -- after each cycle.
   Cycle a ->
   -- | A 'Monitor' observing the chain.
   Monitor a ->
@@ -127,14 +154,14 @@
   a ->
   -- | Number of burn in iterations; deactivate burn in with 'Nothing'.
   Maybe Int ->
-  -- | Auto tuning period (only during burn in); deactivate
-  -- auto tuning with 'Nothing'.
+  -- | Auto tuning period (only during burn in); deactivate auto tuning with
+  -- 'Nothing'.
   Maybe Int ->
-  -- | Number of normal iterations excluding burn in. Note
-  -- that auto tuning only happens during burn in.
+  -- | Number of normal iterations excluding burn in. Note that auto tuning only
+  -- happens during burn in.
   Int ->
-  -- | A source of randomness. For reproducible runs, make
-  -- sure to use generators with the same, fixed seed.
+  -- | A source of randomness. For reproducible runs, make sure to use
+  -- generators with the same, fixed seed.
   GenIO ->
   Status a
 status n p l c m x mB mT nI g
@@ -157,11 +184,17 @@
       Nothing
       p
       l
+      Nothing
       c
       m
   where
     i = Item x (p x) (l x)
 
+-- | Clean the state every given number of generations using the given function.
+-- See 'Cleaner'.
+cleanWith :: Cleaner a -> Status a -> Status a
+cleanWith c s = s {cleaner = Just c}
+
 -- | Save the Markov chain with trace of given length.
 saveWith :: Int -> Status a -> Status a
 saveWith n s = s {save = Just n}
@@ -179,3 +212,7 @@
 -- | Be verbose.
 debug :: Status a -> Status a
 debug s = s {verbosity = Debug}
+
+-- | Set the likelihood function to 1.0. Useful for debugging and testing.
+noData :: Status a -> Status a
+noData s = s {likelihoodF = const 1.0}
diff --git a/src/Mcmc/Tools/Shuffle.hs b/src/Mcmc/Tools/Shuffle.hs
deleted file mode 100644
--- a/src/Mcmc/Tools/Shuffle.hs
+++ /dev/null
@@ -1,64 +0,0 @@
--- |
--- Module      :  Mcmc.Tools.Shuffle
--- Description :  Shuffle a list
--- Copyright   :  (c) Dominik Schrempf 2020
--- License     :  GPL-3.0-or-later
---
--- Maintainer  :  dominik.schrempf@gmail.com
--- Stability   :  unstable
--- Portability :  portable
---
--- Creation date: Wed May 20 14:37:09 2020.
---
--- From https://wiki.haskell.org/Random_shuffle.
-module Mcmc.Tools.Shuffle
-  ( shuffle,
-    shuffleN,
-    grabble,
-  )
-where
-
-import Control.Monad
-import Control.Monad.ST
-import Data.Vector (Vector)
-import qualified Data.Vector as V
-import qualified Data.Vector.Mutable as M
-import System.Random.MWC
-  ( GenIO,
-    uniformR,
-  )
-
--- | Shuffle a list.
-shuffle :: [a] -> GenIO -> IO [a]
-shuffle xs g = head <$> grabble xs 1 (length xs) g
-
--- | Shuffle a list @n@ times.
-shuffleN :: [a] -> Int -> GenIO -> IO [[a]]
-shuffleN xs n = grabble xs n (length xs)
-
--- -- Using System.Random.Shuffle. Speed is the same, so stay without additional dependency.
--- -- | Shuffle a list @n@ times.
--- shuffleN :: [a] -> Int -> GenIO -> IO [[a]]
--- shuffleN xs n g = replicateM n $ fmap (shuffle xs) (rseqM (length xs - 1) g)
---   where
---     rseqM :: Int -> GenIO -> IO [Int]
---     rseqM 0 _ = return []
---     rseqM i gen = liftM2 (:) (uniformR (0, i) gen) (rseqM (i - 1) gen)
-
--- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@
--- elements from @xs@, without replacement, and that @m@ times.
-grabble :: [a] -> Int -> Int -> GenIO -> IO [[a]]
-grabble xs m n gen = do
-  swapss <- replicateM m $
-    forM [0 .. min (l - 1) n] $ \i -> do
-      j <- uniformR (i, l) gen
-      return (i, j)
-  return $ map (V.toList . V.take n . swapElems (V.fromList xs)) swapss
-  where
-    l = length xs - 1
-
-swapElems :: Vector a -> [(Int, Int)] -> Vector a
-swapElems xs swaps = runST $ do
-  mxs <- V.unsafeThaw xs
-  mapM_ (uncurry $ M.unsafeSwap mxs) swaps
-  V.unsafeFreeze mxs
diff --git a/test/Mcmc/ProposalSpec.hs b/test/Mcmc/ProposalSpec.hs
--- a/test/Mcmc/ProposalSpec.hs
+++ b/test/Mcmc/ProposalSpec.hs
@@ -20,10 +20,10 @@
 import Test.Hspec
 
 p1 :: Proposal Double
-p1 = slideSymmetric 1.0 "test1" 1 True
+p1 = slideSymmetric 1.0 (PName "Test 1") (PWeight 1) Tune
 
 p2 :: Proposal Double
-p2 = slideSymmetric 1.0 "test2" 3 True
+p2 = slideSymmetric 1.0 (PName "Test 2") (PWeight 3) Tune
 
 c :: Cycle Double
 c = fromList [p1, p2]
diff --git a/test/Mcmc/SaveSpec.hs b/test/Mcmc/SaveSpec.hs
--- a/test/Mcmc/SaveSpec.hs
+++ b/test/Mcmc/SaveSpec.hs
@@ -38,10 +38,10 @@
 proposals :: Cycle Double
 proposals =
   fromList
-    [ slideSymmetric 0.1 "small" 5 True,
-      slideSymmetric 1.0 "medium" 2 True,
-      slideSymmetric 5.0 "large" 2 True,
-      slide 1.0 4.0 "skewed" 1 True
+    [ slideSymmetric 0.1 (PName "Small") (PWeight 5) Tune,
+      slideSymmetric 1.0 (PName "Medium") (PWeight 2) Tune,
+      slideSymmetric 5.0 (PName "Large") (PWeight 2) Tune,
+      slide 1.0 4.0 (PName "Skewed") (PWeight 1) Tune
     ]
 
 monStd :: MonitorStdOut Double
@@ -71,7 +71,7 @@
                   saveWith 100 $
                     status "SaveSpec" (const 1) lh proposals mon 0 nBurn nAutoTune nIter gen
         saveStatus "SaveSpec.json" s
-        s' <- loadStatus (const 1) lh proposals mon "SaveSpec.json"
+        s' <- loadStatus (const 1) lh proposals mon Nothing "SaveSpec.json"
         r <- mh s
         r' <- mh s'
         -- Done during 'loadStatus'.