diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -4,6 +4,24 @@
 
 ## Unreleased changes
 
+
+## 0.7.0.0
+
+-   Use random-1.2. This means, that the random number generator type has changed
+    from `GenIO` (`System.Random.MWC`) to `IOGenM StdGen`
+    (`System.Random.Stateful`).
+-   Always try to save chain, not only on `UserInterupt`.
+-   Hamiltonian proposals:
+    -   Split `HSettings` into `HParams` and `HStructure` (preparation for NUTS).
+    -   Rename `HTune` to `HTuningConf`.
+    -   No-U-Turn sampler (NUTS).
+-   Significant changes to `Proposal` and `Tuner`.
+    -   `PResult`.
+    -   `ProposalSimple` -> `PFunction`; Return type is `(PResult, Maybe AcceptanceCounts)`.
+
+
+## 0.6.2.5
+
 -   Provide unsafe loader for MHG algorithm (useful for initializing chains with
     different prior or likelihood functions from saves).
 -   Documentation and readme.
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      :  Main
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -15,11 +15,11 @@
 
 import Criterion.Main
 import Data.List
-import Mcmc.Internal.Gamma
+import Mcmc.Internal.SpecFunctions
 import Normal
 import Numeric.SpecFunctions
 import Poisson
-import System.Random.MWC
+import System.Random.Stateful
 
 gammaBenchG :: RealFloat a => (a -> a) -> [a] -> a
 gammaBenchG f = foldl' (\acc x -> acc + (f x)) 0
@@ -30,7 +30,7 @@
 
 main :: IO ()
 main = do
-  g <- create
+  g <- newIOGenM $ mkStdGen 0
   defaultMain
     [ bgroup
         "Normal"
@@ -38,8 +38,12 @@
           bench "Bactrian" $ nfIO (normalBactrianBench g),
           bench "LargeCycle" $ nfIO (normalLargeCycleBench g)
         ],
-      bench "Poisson" $ nfIO (poissonBench g),
       bgroup
+        "Poisson"
+        [ bench "Random walk" $ nfIO (poissonBench g),
+          bench "Hamiltonian" $ nfIO (poissonHamiltonianBench g)
+        ],
+      bgroup
         "MC3"
         [ bench "MC3 2" $ nfIO (normalMC3 g 2),
           bench "MC3 3" $ nfIO (normalMC3 g 3),
@@ -57,58 +61,76 @@
     ]
 
 -- benchmarking Normal/Slide
--- time                 42.31 ms   (41.88 ms .. 42.60 ms)
---                      1.000 R²   (0.999 R² .. 1.000 R²)
--- mean                 42.75 ms   (42.52 ms .. 43.29 ms)
--- std dev              661.4 μs   (347.4 μs .. 1.074 ms)
+-- time                 24.12 ms   (23.52 ms .. 24.76 ms)
+--                      0.997 R²   (0.995 R² .. 0.999 R²)
+-- mean                 23.70 ms   (23.44 ms .. 24.02 ms)
+-- std dev              658.2 μs   (448.5 μs .. 837.2 μs)
 
 -- benchmarking Normal/Bactrian
--- time                 45.51 ms   (45.30 ms .. 45.92 ms)
---                      1.000 R²   (0.999 R² .. 1.000 R²)
--- mean                 45.41 ms   (45.31 ms .. 45.61 ms)
--- std dev              276.0 μs   (141.5 μs .. 460.8 μs)
+-- time                 26.65 ms   (26.08 ms .. 27.26 ms)
+--                      0.998 R²   (0.996 R² .. 1.000 R²)
+-- mean                 26.53 ms   (26.32 ms .. 26.91 ms)
+-- std dev              626.9 μs   (322.3 μs .. 1.040 ms)
 
 -- benchmarking Normal/LargeCycle
--- time                 68.82 ms   (67.18 ms .. 70.81 ms)
---                      0.999 R²   (0.997 R² .. 1.000 R²)
--- mean                 67.68 ms   (67.26 ms .. 68.59 ms)
--- std dev              1.074 ms   (618.2 μs .. 1.602 ms)
+-- time                 52.48 ms   (51.51 ms .. 53.83 ms)
+--                      0.999 R²   (0.998 R² .. 1.000 R²)
+-- mean                 52.54 ms   (51.62 ms .. 54.90 ms)
+-- std dev              2.765 ms   (758.0 μs .. 4.788 ms)
+-- variance introduced by outliers: 14% (moderately inflated)
 
--- benchmarking Poisson
--- time                 72.94 ms   (63.07 ms .. 87.73 ms)
---                      0.953 R²   (0.920 R² .. 1.000 R²)
--- mean                 64.76 ms   (62.84 ms .. 71.62 ms)
--- std dev              5.785 ms   (783.9 μs .. 10.08 ms)
--- variance introduced by outliers: 26% (moderately inflated)
+-- benchmarking Poisson/Random walk
+-- time                 86.57 ms   (86.06 ms .. 86.97 ms)
+--                      1.000 R²   (1.000 R² .. 1.000 R²)
+-- mean                 86.57 ms   (86.20 ms .. 87.14 ms)
+-- std dev              775.9 μs   (415.4 μs .. 1.205 ms)
 
+-- benchmarking Poisson/Hamiltonian
+-- time                 2.890 s    (2.313 s .. 3.359 s)
+--                      0.993 R²   (0.992 R² .. 1.000 R²)
+-- mean                 3.082 s    (2.960 s .. 3.153 s)
+-- std dev              120.4 ms   (39.93 ms .. 164.0 ms)
+-- variance introduced by outliers: 19% (moderately inflated)
+
 -- benchmarking MC3/MC3 2
--- time                 13.08 ms   (12.73 ms .. 13.44 ms)
---                      0.993 R²   (0.986 R² .. 0.997 R²)
--- mean                 13.41 ms   (13.16 ms .. 13.72 ms)
--- std dev              682.7 μs   (520.3 μs .. 874.0 μs)
--- variance introduced by outliers: 22% (moderately inflated)
+-- time                 7.771 ms   (7.573 ms .. 8.008 ms)
+--                      0.996 R²   (0.993 R² .. 0.999 R²)
+-- mean                 7.810 ms   (7.722 ms .. 7.974 ms)
+-- std dev              328.3 μs   (223.0 μs .. 517.6 μs)
+-- variance introduced by outliers: 19% (moderately inflated)
 
 -- benchmarking MC3/MC3 3
--- time                 19.19 ms   (18.86 ms .. 19.59 ms)
---                      0.998 R²   (0.996 R² .. 1.000 R²)
--- mean                 19.28 ms   (19.11 ms .. 19.51 ms)
--- std dev              454.1 μs   (339.1 μs .. 608.4 μs)
+-- time                 10.69 ms   (10.57 ms .. 10.87 ms)
+--                      0.997 R²   (0.993 R² .. 1.000 R²)
+-- mean                 10.73 ms   (10.65 ms .. 10.90 ms)
+-- std dev              308.1 μs   (132.2 μs .. 528.9 μs)
 
 -- benchmarking MC3/MC3 4
--- time                 25.01 ms   (24.21 ms .. 25.66 ms)
---                      0.997 R²   (0.996 R² .. 0.999 R²)
--- mean                 24.21 ms   (23.99 ms .. 24.55 ms)
--- std dev              606.1 μs   (414.7 μs .. 738.5 μs)
+-- time                 13.81 ms   (13.66 ms .. 14.08 ms)
+--                      0.999 R²   (0.997 R² .. 1.000 R²)
+-- mean                 13.73 ms   (13.67 ms .. 13.85 ms)
+-- std dev              213.2 μs   (95.38 μs .. 334.4 μs)
 
 -- benchmarking MC3/MC3 5
--- time                 28.39 ms   (26.99 ms .. 29.34 ms)
---                      0.995 R²   (0.990 R² .. 0.999 R²)
--- mean                 31.13 ms   (30.09 ms .. 33.49 ms)
--- std dev              3.009 ms   (984.0 μs .. 4.319 ms)
--- variance introduced by outliers: 40% (moderately inflated)
+-- time                 16.80 ms   (16.65 ms .. 16.93 ms)
+--                      0.999 R²   (0.997 R² .. 1.000 R²)
+-- mean                 16.99 ms   (16.87 ms .. 17.23 ms)
+-- std dev              376.6 μs   (162.3 μs .. 589.3 μs)
 
 -- benchmarking MC3/MC3 10
--- time                 57.25 ms   (56.98 ms .. 57.61 ms)
+-- time                 32.59 ms   (32.10 ms .. 33.20 ms)
+--                      0.999 R²   (0.998 R² .. 1.000 R²)
+-- mean                 32.21 ms   (32.07 ms .. 32.46 ms)
+-- std dev              390.9 μs   (219.6 μs .. 640.0 μs)
+
+-- benchmarking probability/gamma function general
+-- time                 60.73 ms   (60.37 ms .. 61.26 ms)
 --                      1.000 R²   (1.000 R² .. 1.000 R²)
--- mean                 57.46 ms   (57.34 ms .. 57.56 ms)
--- std dev              192.4 μs   (140.1 μs .. 284.0 μs)
+-- mean                 60.70 ms   (60.49 ms .. 60.91 ms)
+-- std dev              385.5 μs   (287.3 μs .. 513.5 μs)
+
+-- benchmarking probability/gamma function specialized
+-- time                 59.00 ms   (58.74 ms .. 59.32 ms)
+--                      1.000 R²   (1.000 R² .. 1.000 R²)
+-- mean                 58.83 ms   (58.76 ms .. 58.97 ms)
+-- std dev              168.3 μs   (105.7 μs .. 245.7 μs)
diff --git a/bench/Normal.hs b/bench/Normal.hs
--- a/bench/Normal.hs
+++ b/bench/Normal.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Normal
 -- Description :  Benchmark Metropolis-Hastings-Green algorithm
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -19,7 +19,7 @@
 
 import Control.Monad
 import Mcmc
-import System.Random.MWC
+import System.Random.Stateful
 
 trueMean :: Double
 trueMean = 5
@@ -42,7 +42,7 @@
 mon :: Monitor Double
 mon = Monitor monStd [] []
 
-normalSlideBench :: GenIO -> IO ()
+normalSlideBench :: IOGenM StdGen -> IO ()
 normalSlideBench g = do
   let s =
         Settings
@@ -64,7 +64,7 @@
     [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (pWeight 1) Tune | i <- [0 .. 100 :: Int]]
 
 -- Should have the same run time as 'normalSlide'.
-normalLargeCycleBench :: GenIO -> IO ()
+normalLargeCycleBench :: IOGenM StdGen -> IO ()
 normalLargeCycleBench g = do
   let s =
         Settings
@@ -83,7 +83,7 @@
 ccBactrian :: Cycle Double
 ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (pWeight 1) Tune]
 
-normalBactrianBench :: GenIO -> IO ()
+normalBactrianBench :: IOGenM StdGen -> IO ()
 normalBactrianBench g = do
   let s =
         Settings
@@ -99,7 +99,7 @@
   a <- mhg s noPrior lh ccBactrian mon 0 g
   void $ mcmc s a
 
-normalMC3 :: GenIO -> Int -> IO ()
+normalMC3 :: IOGenM StdGen -> Int -> IO ()
 normalMC3 g n = do
   let mcmcS =
         Settings
diff --git a/bench/Poisson.hs b/bench/Poisson.hs
--- a/bench/Poisson.hs
+++ b/bench/Poisson.hs
@@ -1,7 +1,14 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# OPTIONS_GHC -Wno-orphans #-}
+
 -- |
 -- Module      :  Poisson
 -- Description :  Poisson regression model for airline fatalities
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -13,16 +20,26 @@
 -- See https://revbayes.github.io/tutorials/mcmc/poisson.html.
 module Poisson
   ( poissonBench,
+    poissonHamiltonianBench,
   )
 where
 
 import Control.Monad
-import Lens.Micro
+import Data.Aeson
+import Data.Typeable
+import qualified Data.Vector.Fixed as VB
+import qualified Data.Vector.Fixed.Boxed as VB
+import qualified Data.Vector.Storable as VS
 import Mcmc
-import System.Random.MWC
+import System.Random.Stateful
 
-type I = (Double, Double)
+type IG = VB.Vec2
 
+type I = IG Double
+
+instance (VB.Arity n, ToJSON a) => ToJSON (VB.Vec n a) where
+  toJSON = toJSON . VB.toList
+
 fatalities :: [Int]
 fatalities = [24, 25, 31, 31, 22, 21, 26, 20, 16, 22]
 
@@ -32,43 +49,48 @@
     ys = [1976.0 .. 1985.0]
     m = sum ys / fromIntegral (length ys)
 
-f :: Int -> Double -> I -> Likelihood
-f ft yr (a, b) = poisson l (fromIntegral ft)
+f :: (RealFloat a, Typeable a) => Int -> Double -> IG a -> LikelihoodG a
+f ft yr xs = poisson l ft
   where
-    l = exp $ a + b * yr
+    a = xs VB.! 0
+    b = xs VB.! 1
+    l = exp $ a + b * realToFrac yr
 
-lh :: LikelihoodFunction I
+lh :: (RealFloat a, Typeable a) => LikelihoodFunctionG (IG a) a
 lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]
 
 proposalAlpha :: Proposal I
-proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune
+proposalAlpha = (VB.element 0) @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune
 
 proposalBeta :: Proposal I
-proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune
+proposalBeta = (VB.element 1) @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune
 
 proposals :: Cycle I
 proposals = cycleFromList [proposalAlpha, proposalBeta]
 
 initial :: I
-initial = (0, 0)
+initial = VB.mk2 0 0
 
 monAlpha :: MonitorParameter I
-monAlpha = fst >$< monitorDouble "alpha"
+monAlpha = (VB.! 0) >$< monitorDouble "alpha"
 
 monBeta :: MonitorParameter I
-monBeta = snd >$< monitorDouble "beta"
+monBeta = (VB.! 1) >$< monitorDouble "beta"
 
 monStd :: MonitorStdOut I
 monStd = monitorStdOut [monAlpha, monBeta] 150
 
+-- monF :: MonitorFile I
+-- monF = monitorFile "params" [monAlpha, monBeta] 10
+
 mon :: Monitor I
 mon = Monitor monStd [] []
 
-poissonBench :: GenIO -> IO ()
+poissonBench :: IOGenM StdGen -> IO ()
 poissonBench g = do
   let s =
         Settings
-          (AnalysisName "Poisson")
+          (AnalysisName "PoissonR")
           (BurnInWithAutoTuning 2000 200)
           (Iterations 10000)
           TraceAuto
@@ -78,4 +100,41 @@
           LogStdOutOnly
           Quiet
   a <- mhg s noPrior lh proposals mon initial g
+  void $ mcmc s a
+
+toVec :: I -> VS.Vector Double
+toVec xs = VS.generate 2 (\i -> xs VB.! i)
+
+fromVec :: I -> VS.Vector Double -> I
+fromVec _ xs = VB.mk2 (xs VS.! 0) (xs VS.! 1)
+
+hparams :: HParams
+hparams = HParams Nothing Nothing Nothing
+
+htconf :: HTuningConf
+htconf = HTuningConf HTuneLeapfrog HTuneAllMasses
+
+hstruct :: HStructure IG
+hstruct = HStructure initial toVec fromVec
+
+target :: HTarget IG
+target = HTarget Nothing lh Nothing
+
+hmcProposal :: Cycle I
+hmcProposal = cycleFromList [hamiltonian hparams htconf hstruct target (PName "Hamiltonian") (pWeight 1)]
+
+poissonHamiltonianBench :: IOGenM StdGen -> IO ()
+poissonHamiltonianBench g = do
+  let s =
+        Settings
+          (AnalysisName "PoissonH")
+          (BurnInWithAutoTuning 2000 200)
+          (Iterations 10000)
+          TraceAuto
+          Overwrite
+          Sequential
+          NoSave
+          LogStdOutOnly
+          Quiet
+  a <- mhg s noPrior lh hmcProposal mon initial g
   void $ mcmc s a
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.6.2.5
+version:            0.7.0.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>
@@ -10,7 +10,7 @@
 bug-reports:        https://github.com/dschrempf/mcmc/issues
 author:             Dominik Schrempf
 maintainer:         dominik.schrempf@gmail.com
-copyright:          Dominik Schrempf (2021)
+copyright:          2021 Dominik Schrempf
 license:            GPL-3.0-or-later
 build-type:         Simple
 extra-source-files:
@@ -53,7 +53,8 @@
     Mcmc.Chain.Trace
     Mcmc.Cycle
     Mcmc.Environment
-    Mcmc.Internal.Gamma
+    Mcmc.Internal.SpecFunctions
+    Mcmc.Jacobian
     Mcmc.Likelihood
     Mcmc.Logger
     Mcmc.MarginalLikelihood
@@ -68,7 +69,9 @@
     Mcmc.Proposal
     Mcmc.Proposal.Bactrian
     Mcmc.Proposal.Generic
-    Mcmc.Proposal.Hamiltonian
+    Mcmc.Proposal.Hamiltonian.Common
+    Mcmc.Proposal.Hamiltonian.Hamiltonian
+    Mcmc.Proposal.Hamiltonian.Nuts
     Mcmc.Proposal.Scale
     Mcmc.Proposal.Simplex
     Mcmc.Proposal.Slide
@@ -79,11 +82,14 @@
     Mcmc.Internal.ByteString
     Mcmc.Internal.Random
     Mcmc.Internal.Shuffle
+    Mcmc.Proposal.Hamiltonian.Internal
+    Mcmc.Proposal.Hamiltonian.Masses
     Paths_mcmc
 
   autogen-modules: Paths_mcmc
   hs-source-dirs:  src
   build-depends:
+    , ad
     , aeson
     , async
     , base               >=4.7 && <5
@@ -92,7 +98,6 @@
     , containers
     , covariance         >=0.2
     , data-default
-    , deepseq
     , directory
     , dirichlet
     , double-conversion
@@ -104,6 +109,8 @@
     , parallel
     , pretty-show
     , primitive
+    , random
+    , splitmix
     , statistics
     , time
     , transformers
@@ -124,7 +131,7 @@
     , base        >=4.7 && <5
     , hspec
     , mcmc
-    , mwc-random
+    , random
     , statistics
 
 benchmark mcmc-bench
@@ -139,9 +146,11 @@
   autogen-modules: Paths_mcmc
   hs-source-dirs:  bench
   build-depends:
+    , aeson
     , base            >=4.7 && <5
     , criterion
+    , fixed-vector
     , math-functions
     , mcmc
-    , microlens
-    , mwc-random
+    , random
+    , vector
diff --git a/src/Mcmc.hs b/src/Mcmc.hs
--- a/src/Mcmc.hs
+++ b/src/Mcmc.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc
 -- Description :  Markov chain Monte Carlo samplers, batteries included
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -110,8 +110,8 @@
     -- tuning parameter @t@ is implemented as
     --
     -- @
-    -- slideSimple :: Double -> Double -> Double -> ProposalSimple Double
-    -- slideSimple m s t =
+    -- slide :: Double -> Double -> Double -> PFunction Double
+    -- slide m s t =
     --   genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing
     -- @
     --
@@ -129,7 +129,6 @@
     PWeight (..),
     pWeight,
     Proposal,
-    JacobianFunction,
     (@~),
     liftProposal,
     liftProposalWith,
@@ -143,8 +142,10 @@
     slideUniformSymmetric,
     slideContrarily,
     slideBactrian,
-    module Mcmc.Proposal.Hamiltonian,
     module Mcmc.Proposal.Simplex,
+    module Mcmc.Proposal.Hamiltonian.Common,
+    module Mcmc.Proposal.Hamiltonian.Hamiltonian,
+    module Mcmc.Proposal.Hamiltonian.Nuts,
 
     -- ** Cycles
     Cycle,
@@ -181,9 +182,10 @@
     module Mcmc.Monitor.Parameter,
     module Mcmc.Monitor.ParameterBatch,
 
-    -- * Priors, likelihoodl, and posteriors
+    -- * Priors, likelihoods, Jacobians, and posteriors
     module Mcmc.Prior,
     module Mcmc.Likelihood,
+    module Mcmc.Jacobian,
     module Mcmc.Posterior,
 
     -- * Run MCMC samplers
@@ -209,6 +211,7 @@
 import Mcmc.Algorithm.MC3
 import Mcmc.Algorithm.MHG
 import Mcmc.Cycle
+import Mcmc.Jacobian
 import Mcmc.Likelihood
 import Mcmc.MarginalLikelihood
 import Mcmc.Mcmc
@@ -219,7 +222,9 @@
 import Mcmc.Prior
 import Mcmc.Proposal
 import Mcmc.Proposal.Bactrian
-import Mcmc.Proposal.Hamiltonian
+import Mcmc.Proposal.Hamiltonian.Common
+import Mcmc.Proposal.Hamiltonian.Hamiltonian
+import Mcmc.Proposal.Hamiltonian.Nuts
 import Mcmc.Proposal.Scale
 import Mcmc.Proposal.Simplex
 import Mcmc.Proposal.Slide
diff --git a/src/Mcmc/Acceptance.hs b/src/Mcmc/Acceptance.hs
--- a/src/Mcmc/Acceptance.hs
+++ b/src/Mcmc/Acceptance.hs
@@ -1,7 +1,9 @@
+{-# LANGUAGE TemplateHaskell #-}
+
 -- |
 -- Module      :  Mcmc.Acceptance
 -- Description :  Handle acceptance rates
--- Copyright   :  (c) 2021 Dominik Schrempf
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -12,9 +14,12 @@
 module Mcmc.Acceptance
   ( -- * Acceptance rates
     AcceptanceRate,
+    AcceptanceCounts (..),
     Acceptance (fromAcceptance),
     emptyA,
-    pushA,
+    pushAccept,
+    pushReject,
+    pushAcceptanceCounts,
     resetA,
     transformKeysA,
     acceptanceRate,
@@ -22,18 +27,36 @@
   )
 where
 
-import Control.DeepSeq
 import Data.Aeson
-import Data.Bifunctor
+import Data.Aeson.TH
 import Data.Foldable
 import qualified Data.Map.Strict as M
 
 -- | Acceptance rate.
 type AcceptanceRate = Double
 
+-- | Number of accepted and rejected proposals.
+data AcceptanceCounts = AcceptanceCounts
+  { nAccepted :: !Int,
+    nRejected :: !Int
+  }
+  deriving (Show, Eq, Ord)
+
+$(deriveJSON defaultOptions ''AcceptanceCounts)
+
+addAccept :: AcceptanceCounts -> AcceptanceCounts
+addAccept (AcceptanceCounts a r) = AcceptanceCounts (a + 1) r
+
+addReject :: AcceptanceCounts -> AcceptanceCounts
+addReject (AcceptanceCounts a r) = AcceptanceCounts a (r + 1)
+
+addAcceptanceCounts :: AcceptanceCounts -> AcceptanceCounts -> AcceptanceCounts
+addAcceptanceCounts (AcceptanceCounts al rl) (AcceptanceCounts ar rr) =
+  AcceptanceCounts (al + ar) (rl + rr)
+
 -- | For each key @k@, store the number of accepted and rejected proposals.
-newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)}
-  deriving (Eq, Read, Show)
+newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k AcceptanceCounts}
+  deriving (Eq, Show)
 
 instance ToJSONKey k => ToJSON (Acceptance k) where
   toJSON (Acceptance m) = toJSON m
@@ -46,13 +69,20 @@
 --
 -- Initialize an empty storage of accepted/rejected values.
 emptyA :: Ord k => [k] -> Acceptance k
-emptyA ks = Acceptance $ M.fromList [(k, (0, 0)) | k <- ks]
+emptyA ks = Acceptance $ M.fromList [(k, AcceptanceCounts 0 0) | k <- ks]
 
--- | For key @k@, prepend an accepted (True) or rejected (False) proposal.
-pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k
-pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance
-pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance
+-- | For key @k@, add an accept.
+pushAccept :: Ord k => k -> Acceptance k -> Acceptance k
+pushAccept k = Acceptance . M.adjust addAccept k . fromAcceptance
 
+-- | For key @k@, add a reject.
+pushReject :: Ord k => k -> Acceptance k -> Acceptance k
+pushReject k = Acceptance . M.adjust addReject k . fromAcceptance
+
+-- | For key @k@, add acceptance counts.
+pushAcceptanceCounts :: Ord k => k -> AcceptanceCounts -> Acceptance k -> Acceptance k
+pushAcceptanceCounts k c = Acceptance . M.adjust (addAcceptanceCounts c) k . fromAcceptance
+
 -- | Reset acceptance storage.
 resetA :: Ord k => Acceptance k -> Acceptance k
 resetA = emptyA . M.keys . fromAcceptance
@@ -73,8 +103,8 @@
 -- zero).
 acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, AcceptanceRate)
 acceptanceRate k a = case fromAcceptance a M.!? k of
-  Just (0, 0) -> Nothing
-  Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))
+  Just (AcceptanceCounts 0 0) -> Nothing
+  Just (AcceptanceCounts as rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))
   Nothing -> error "acceptanceRate: Key not found in map."
 
 -- | Acceptance rates for all proposals.
@@ -84,7 +114,7 @@
 acceptanceRates :: Acceptance k -> M.Map k (Maybe AcceptanceRate)
 acceptanceRates =
   M.map
-    ( \(as, rs) ->
+    ( \(AcceptanceCounts as rs) ->
         if as + rs == 0
           then Nothing
           else Just $ fromIntegral as / fromIntegral (as + rs)
diff --git a/src/Mcmc/Algorithm.hs b/src/Mcmc/Algorithm.hs
--- a/src/Mcmc/Algorithm.hs
+++ b/src/Mcmc/Algorithm.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Algorithm
 -- Description :  Algortihms for Markov chain Monte Carlo samplers
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -17,6 +17,7 @@
 import qualified Data.ByteString.Lazy.Char8 as BL
 import Data.Time
 import Mcmc.Cycle
+import Mcmc.Proposal
 import Mcmc.Settings
 
 -- | Class for algorithms used by MCMC samplers.
@@ -38,10 +39,14 @@
   --
   -- NOTE: Computation in the 'IO' Monad is necessary because the trace is
   -- mutable.
-  aAutoTune :: Int -> a -> IO a
+  aAutoTune :: TuningType -> Int -> a -> IO a
 
   -- | Reset acceptance counts.
   aResetAcceptance :: a -> a
+
+  -- | Clean after burn in. In particular, this is used to reduce the length of
+  -- the trace, if required.
+  aCleanAfterBurnIn :: TraceLength -> a -> IO a
 
   -- | Summarize the cycle.
   aSummarizeCycle :: IterationMode -> a -> BL.ByteString
diff --git a/src/Mcmc/Algorithm/MC3.hs b/src/Mcmc/Algorithm/MC3.hs
--- a/src/Mcmc/Algorithm/MC3.hs
+++ b/src/Mcmc/Algorithm/MC3.hs
@@ -5,7 +5,7 @@
 -- |
 -- Module      :  Mcmc.Algorithm.MC3
 -- Description :  Metropolis-coupled Markov chain Monte Carlo algorithm
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -79,9 +79,7 @@
 import Mcmc.Proposal
 import Mcmc.Settings
 import Numeric.Log hiding (sum)
-import System.Random.MWC
-
--- import Debug.Trace hiding (trace)
+import System.Random.Stateful
 
 -- | Total number of parallel chains.
 --
@@ -130,9 +128,9 @@
     savedMC3ReciprocalTemperatures :: ReciprocalTemperatures,
     savedMC3Iteration :: Int,
     savedMC3SwapAcceptance :: Acceptance Int,
-    savedMC3Generator :: U.Vector Word32
+    savedMC3Generator :: (Word64, Word64)
   }
-  deriving (Eq, Read, Show)
+  deriving (Eq, Show)
 
 $(deriveJSON defaultOptions ''SavedMC3)
 
@@ -175,7 +173,7 @@
     mc3Iteration :: Int,
     -- | Number of accepted and rejected swaps.
     mc3SwapAcceptance :: Acceptance Int,
-    mc3Generator :: GenIO
+    mc3Generator :: IOGenM StdGen
   }
 
 instance ToJSON a => Algorithm (MC3 a) where
@@ -185,6 +183,7 @@
   aIterate = mc3Iterate
   aAutoTune = mc3AutoTune
   aResetAcceptance = mc3ResetAcceptance
+  aCleanAfterBurnIn = mc3CleanAfterBurnIn
   aSummarizeCycle = mc3SummarizeCycle
   aOpenMonitors = mc3OpenMonitors
   aExecuteMonitors = mc3ExecuteMonitors
@@ -266,8 +265,9 @@
     l = link c
     t = trace c
 
--- TODO: Splitmix. Initialization of the MC3 algorithm is an IO action because
--- the generators have to be split. And also because of the mutable trace.
+-- TODO (medium): Splitmix. Initialization of the MC3 algorithm is an IO action
+-- because the generators have to be split. And also because of the mutable
+-- trace.
 
 -- | Initialize an MC3 algorithm with a given number of chains.
 --
@@ -284,7 +284,7 @@
   Cycle a ->
   Monitor a ->
   InitialState a ->
-  GenIO ->
+  IOGenM StdGen ->
   IO (MC3 a)
 mc3 sMc3 s pr lh cc mn i0 g
   | n < 2 = error "mc3: The number of chains must be two or larger."
@@ -292,7 +292,8 @@
   | sn < 1 || sn > n - 1 = error "mc3: The number of swaps must be in [1, NChains - 1]."
   | otherwise = do
       -- Split random number generators.
-      gs <- V.fromList <$> splitGen n g
+      rs <- replicateM n $ splitGenM g
+      gs <- V.fromList <$> mapM newIOGenM rs
       cs <- V.mapM (mhg s pr lh cc mn i0) gs
       hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs
       return $ MC3 sMc3 hcs bs 0 (emptyA [0 .. n - 2]) g
@@ -407,10 +408,10 @@
       -- traceIO $ "Log priors (left, right, after swap): " <> show (ln prL') <> " " <> show (ln prR')
       -- traceIO $ "Log likelihoods (left, right, before swap): " <> show (ln lhL) <> " " <> show (ln lhR)
       -- traceIO $ "Log likelihood (left, right, after swap): " <> show (ln lhL') <> " " <> show (ln lhR')
-      let !ac' = pushA i True (mc3SwapAcceptance a)
+      let !ac' = pushAccept i (mc3SwapAcceptance a)
       return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'}
     else do
-      let !ac' = pushA i False (mc3SwapAcceptance a)
+      let !ac' = pushReject i (mc3SwapAcceptance a)
       return $ a {mc3SwapAcceptance = ac'}
   where
     g = mc3Generator a
@@ -420,8 +421,9 @@
   where
     mhgs = mc3MHGChains a
 
--- TODO: Splitmix. 'mc3Iterate' is actually not parallel, but concurrent because
--- of the IO constraint. Use pure parallel code when we have a pure generator.
+-- TODO (medium): Splitmix. 'mc3Iterate' is actually not parallel, but
+-- concurrent because of the IO constraint. Use pure parallel code when we have
+-- a pure generator.
 --
 -- However, we have to take care of the mutable traces.
 mc3Iterate ::
@@ -475,10 +477,10 @@
     rNew = (brOld / blOld) ** xi
     brNew = blNew * rNew
 
-mc3AutoTune :: ToJSON a => Int -> MC3 a -> IO (MC3 a)
-mc3AutoTune l a = do
+mc3AutoTune :: ToJSON a => TuningType -> Int -> MC3 a -> IO (MC3 a)
+mc3AutoTune b l a = do
   -- 1. Auto tune all chains.
-  mhgs' <- V.mapM (aAutoTune l) $ mc3MHGChains a
+  mhgs' <- V.mapM (aAutoTune b l) $ mc3MHGChains a
   -- 2. Auto tune temperatures.
   let optimalRate = getOptimalRate PDimensionUnknown
       mCurrentRates = acceptanceRates $ mc3SwapAcceptance a
@@ -519,6 +521,13 @@
     ac' = resetA $ mc3SwapAcceptance a
     --
     a' = a {mc3MHGChains = mhgs', mc3SwapAcceptance = ac'}
+
+mc3CleanAfterBurnIn :: ToJSON a => TraceLength -> MC3 a -> IO (MC3 a)
+mc3CleanAfterBurnIn tl a = do
+  cs' <- V.mapM (aCleanAfterBurnIn tl) cs
+  pure $ a {mc3MHGChains = cs'}
+  where
+    cs = mc3MHGChains a
 
 -- Information in cycle summary:
 --
diff --git a/src/Mcmc/Algorithm/MHG.hs b/src/Mcmc/Algorithm/MHG.hs
--- a/src/Mcmc/Algorithm/MHG.hs
+++ b/src/Mcmc/Algorithm/MHG.hs
@@ -4,7 +4,7 @@
 -- |
 -- Module      :  Mcmc.Algorithm.MHG
 -- Description :  Metropolis-Hastings-Green algorithm
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -34,6 +34,7 @@
 import Control.Parallel.Strategies
 import Data.Aeson
 import qualified Data.ByteString.Lazy.Char8 as BL
+import Data.Maybe
 import Data.Time
 import qualified Data.Vector as VB
 import Mcmc.Acceptance
@@ -50,7 +51,7 @@
 import Mcmc.Proposal
 import Mcmc.Settings
 import Numeric.Log
-import System.Random.MWC
+import System.Random.Stateful
 import Text.Printf
 import Prelude hiding (cycle)
 
@@ -64,6 +65,7 @@
   aIterate = mhgIterate
   aAutoTune = mhgAutoTune
   aResetAcceptance = mhgResetAcceptance
+  aCleanAfterBurnIn = mhgCleanAfterBurnIn
   aSummarizeCycle = mhgSummarizeCycle
   aOpenMonitors = mhgOpenMonitors
   aExecuteMonitors = mhgExecuteMonitors
@@ -71,6 +73,23 @@
   aCloseMonitors = mhgCloseMonitors
   aSave = mhgSave
 
+getTraceLength ::
+  Maybe BurnInSettings ->
+  TraceLength ->
+  Monitor a ->
+  Cycle a ->
+  Int
+getTraceLength burnIn tl mn cc = maximum $ minimumTraceLength : bi : batchMonitorSizes
+  where
+    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn
+    minimumTraceLength = case tl of
+      TraceAuto -> 1
+      TraceMinimum n -> n
+    bi = case (ccRequireTrace cc, burnIn) of
+      (True, Just (BurnInWithAutoTuning _ n)) -> n
+      (True, Just (BurnInWithCustomAutoTuning ns ms)) -> max (maximum $ 0 : ns) (maximum $ 0 : ms)
+      _ -> 0
+
 -- | Initialize an MHG algorithm.
 --
 -- NOTE: Computation in the 'IO' Monad is necessary because the trace is
@@ -82,25 +101,17 @@
   Cycle a ->
   Monitor a ->
   InitialState a ->
-  GenIO ->
+  IOGenM StdGen ->
   IO (MHG a)
 mhg s pr lh cc mn i0 g = do
   -- The trace is a mutable vector and the mutable state needs to be handled by
   -- a monad.
-  tr <- replicateT traceLength l0
+  tr <- replicateT tl l0
   return $ MHG $ Chain Nothing l0 0 tr ac g 0 pr lh cc mn
   where
     l0 = Link i0 (pr i0) (lh i0)
     ac = emptyA $ ccProposals cc
-    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn
-    minimumTraceLength = case sTraceLength s of
-      TraceAuto -> 1
-      TraceMinimum n -> n
-    bi = case sBurnIn s of
-      BurnInWithAutoTuning _ n -> n
-      BurnInWithCustomAutoTuning ns ms -> max (maximum $ 0 : ns) (maximum $ 0 : ms)
-      _ -> 0
-    traceLength = maximum $ minimumTraceLength : bi : batchMonitorSizes
+    tl = getTraceLength (Just $ sBurnIn s) (sTraceLength s) mn cc
 
 mhgFn :: AnalysisName -> FilePath
 mhgFn (AnalysisName nm) = nm ++ ".mcmc.mhg"
@@ -162,11 +173,14 @@
 
 -- The MHG ratio. This implementation has the following properties:
 --
+-- - The kernel ratio and the Jacobian are checked carefully and should be
+-- - strictly positive, finite numbers.
+--
 -- - The ratio is 'Infinity' if fX is zero. In this case, the proposal is always
 --   accepted.
 --
--- - The ratio 'NaN' if (fY or q or j) and fX are zero. In this case, the
---   proposal is always rejected.
+-- - The ratio is 'NaN' if fY and fX are zero. In this case, the proposal is
+--   always rejected.
 --
 -- This means that a chain in a state with posterior probability zero (fX=0) can
 -- only move if a state with non-zero posterior probability is proposed.
@@ -182,49 +196,71 @@
 -- Since I trust the author of Chapter 1 (Charles Geyer) I choose to follow
 -- option (b). However, Option (a) is more user-friendly.
 --
--- [1] Handbook of markov chain monte carlo (2011), CRC press.
+-- [1] Handbook of Markov chain Monte Carlo (2011), CRC press.
 --
--- [2] Dellaportas, P., & Roberts, G. O., An introduction to mcmc, Lecture Notes
+-- [2] Dellaportas, P., & Roberts, G. O., An introduction to MCMC, Lecture Notes
 -- in Statistics, (), 1–41 (2003).
 -- http://dx.doi.org/10.1007/978-0-387-21811-3_1.
 mhgRatio :: Posterior -> Posterior -> KernelRatio -> Jacobian -> MHGRatio
--- q = qYX / qXY * jXY; see 'ProposalSimple'.
--- j = Jacobian.
-mhgRatio fX fY q j = fY / fX * q * j
+mhgRatio fX fY q j
+  | q == 0.0 = error "mhgRatio: Kernel ratio is negative infinity. Use 'ForceReject'."
+  | q == 1.0 / 0.0 = error "mhgRatio: Kernel ratio is infinity. Use 'ForceAccept'."
+  | q == 0.0 / 0.0 = error "mhgRatio: Kernel ratio is NaN."
+  | j == 0.0 = error "mhgRatio: Jacobian is negative infinity. Use 'ForceReject'."
+  | j == 1.0 / 0.0 = error "mhgRatio: Jacobian is infinity. Use 'ForceAccept'."
+  | j == 0.0 / 0.0 = error "mhgRatio: Jacobian is NaN."
+  | otherwise = fY / fX * q * j
 {-# INLINE mhgRatio #-}
 
 -- | Accept or reject a proposal with given MHG ratio?
-mhgAccept :: MHGRatio -> GenIO -> IO Bool
+mhgAccept :: MHGRatio -> IOGenM StdGen -> IO Bool
 mhgAccept r g
   | ln r >= 0.0 = return True
   | otherwise = do
-      b <- uniform g
+      b <- uniformRM (0, 1) g
       return $ b < exp (ln r)
 
 mhgPropose :: MHG a -> Proposal a -> IO (MHG a)
 mhgPropose (MHG c) p = do
   -- 1. Sample new state.
-  (!y, !q, !j) <- liftIO $ s x g
-  -- 2. Calculate Metropolis-Hastings-Green ratio.
+  !(pres, mcs) <- liftIO $ s x g
+  -- 2. Define new prior and likelihood calculation functions. Avoid actual
+  -- calculation of the values.
   --
   -- Most often, parallelization is not helpful, because the prior and
   -- likelihood functions are too fast; see
   -- https://stackoverflow.com/a/46603680/3536806.
-  let (pY, lY) = (pF y, lF y) `using` parTuple2 rdeepseq rdeepseq
-  -- let !pY = pF y
-  --     !lY = lF y
-  let !r = mhgRatio (pX * lX) (pY * lY) q j
+  let calcPrLh y = (pF y, lF y) `using` parTuple2 rdeepseq rdeepseq
+      accept y pr lh =
+        let !ac' = case mcs of
+              Nothing -> pushAccept p ac
+              Just cs -> pushAcceptanceCounts p cs ac
+         in pure $ MHG $ c {link = Link y pr lh, acceptance = ac'}
+      reject =
+        let !ac' = case mcs of
+              Nothing -> pushReject p ac
+              Just cs -> pushAcceptanceCounts p cs ac
+         in pure $ MHG $ c {acceptance = ac'}
   -- 3. Accept or reject.
-  accept <- mhgAccept r g
-  if accept
-    then do
-      let !ac' = pushA p True ac
-      return $ MHG $ c {link = Link y pY lY, acceptance = ac'}
-    else do
-      let !ac' = pushA p False ac
-      return $ MHG $ c {acceptance = pushA p False ac'}
+  --
+  -- 3a. When rejection is inevitable, avoid calculation of the prior, the
+  -- likelihood and the MHG ratio.
+  case pres of
+    ForceReject -> reject
+    ForceAccept y -> let (pY, lY) = calcPrLh y in accept y pY lY
+    (Propose y q j) ->
+      if q <= 0.0 || j <= 0.0
+        then reject
+        else do
+          -- 3b. Calculate Metropolis-Hastings-Green ratio.
+          let (pY, lY) = calcPrLh y
+              !r = mhgRatio (pX * lX) (pY * lY) q j
+          isAccept <- mhgAccept r g
+          if isAccept
+            then accept y pY lY
+            else reject
   where
-    s = prSimple p
+    s = prFunction p
     (Link x pX lX) = link c
     pF = priorFunction c
     lF = likelihoodFunction c
@@ -245,12 +281,13 @@
 -- At the moment this just checks whether the prior, likelihood, or posterior
 -- are NaN or infinite.
 mhgIsInValidState :: MHG a -> Bool
-mhgIsInValidState a = check p || check l || check (p * l)
+mhgIsInValidState a = checkSoft p || check l || check (p * l)
   where
     x = link $ fromMHG a
     p = prior x
     l = likelihood x
     check v = let v' = ln v in isNaN v' || isInfinite v' || v' == 0
+    checkSoft v = let v' = ln v in isNaN v' || isInfinite v'
 
 -- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG
 -- algorithm is just inherently sequential. Parallelization can be achieved by
@@ -266,18 +303,34 @@
     cc = cycle c
     g = generator c
 
-mhgAutoTune :: Int -> MHG a -> IO (MHG a)
-mhgAutoTune n (MHG c) = do
-  tr <- VB.map state <$> takeT n (trace c)
-  return $ MHG $ c {cycle = autoTuneCycle ac tr cc}
+mhgAutoTune :: TuningType -> Int -> MHG a -> IO (MHG a)
+mhgAutoTune b n (MHG c) = do
+  mxs <-
+    if ccRequireTrace cc
+      then Just . VB.map state <$> takeT n tr
+      else pure Nothing
+  return $ MHG $ c {cycle = autoTuneCycle b ac mxs cc}
   where
     ac = acceptance c
     cc = cycle c
+    tr = trace c
 
 mhgResetAcceptance :: MHG a -> MHG a
 mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}
   where
     ac = acceptance c
+
+mhgCleanAfterBurnIn :: TraceLength -> MHG a -> IO (MHG a)
+mhgCleanAfterBurnIn tl (MHG c) = do
+  xs <- takeT l tr
+  tr' <- fromVectorT xs
+  let c' = c {trace = tr'}
+  pure $ MHG c'
+  where
+    mn = monitor c
+    cc = cycle c
+    tr = trace c
+    l = getTraceLength Nothing tl mn cc
 
 mhgSummarizeCycle :: IterationMode -> MHG a -> BL.ByteString
 mhgSummarizeCycle m (MHG c) = summarizeCycle m ac cc
diff --git a/src/Mcmc/Chain/Chain.hs b/src/Mcmc/Chain/Chain.hs
--- a/src/Mcmc/Chain/Chain.hs
+++ b/src/Mcmc/Chain/Chain.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Chain.Chain
 -- Description :  Simple representation of a Markov chain
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -59,7 +59,7 @@
 import Mcmc.Monitor
 import Mcmc.Prior
 import Mcmc.Proposal
-import System.Random.MWC hiding (save)
+import System.Random.Stateful
 import Prelude hiding (cycle)
 
 -- | Type synonym to indicate the initial state.
@@ -98,7 +98,7 @@
     -- order.
     acceptance :: Acceptance (Proposal a),
     -- | The random number generator.
-    generator :: GenIO,
+    generator :: IOGenM StdGen,
     --
     -- Variables and functions; not saved.
 
diff --git a/src/Mcmc/Chain/Link.hs b/src/Mcmc/Chain/Link.hs
--- a/src/Mcmc/Chain/Link.hs
+++ b/src/Mcmc/Chain/Link.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Chain.Link
 -- Description :  The state combined with auxiliary variables
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Chain/Save.hs b/src/Mcmc/Chain/Save.hs
--- a/src/Mcmc/Chain/Save.hs
+++ b/src/Mcmc/Chain/Save.hs
@@ -14,7 +14,7 @@
 -- |
 -- Module      :  Mcmc.Chain.Save
 -- Description :  Save and load a Markov chain
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 module Mcmc.Chain.Save
   ( SavedChain (..),
@@ -31,7 +31,6 @@
 import Data.Maybe
 import qualified Data.Stack.Circular as C
 import qualified Data.Vector as VB
-import qualified Data.Vector.Unboxed as VU
 import Data.Word
 import Mcmc.Acceptance
 import Mcmc.Chain.Chain
@@ -54,10 +53,10 @@
     savedIteration :: Int,
     savedTrace :: C.Stack VB.Vector (Link a),
     savedAcceptance :: Acceptance Int,
-    savedSeed :: VU.Vector Word32,
+    savedSeed :: (Word64, Word64),
     savedTuningParameters :: [Maybe (TuningParameter, AuxiliaryTuningParameters)]
   }
-  deriving (Eq, Read, Show)
+  deriving (Eq, Show)
 
 $(deriveJSON defaultOptions ''SavedChain)
 
diff --git a/src/Mcmc/Chain/Trace.hs b/src/Mcmc/Chain/Trace.hs
--- a/src/Mcmc/Chain/Trace.hs
+++ b/src/Mcmc/Chain/Trace.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Chain.Trace
 -- Description :  History of a Markov chain
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -12,6 +12,7 @@
 module Mcmc.Chain.Trace
   ( Trace,
     replicateT,
+    fromVectorT,
     lengthT,
     pushT,
     headT,
@@ -37,6 +38,12 @@
 -- Call 'error' if the maximum size is zero or negative.
 replicateT :: Int -> Link a -> IO (Trace a)
 replicateT n l = Trace <$> C.replicate n l
+
+-- | Create a trace from a vector. The length is determined by the vector.
+--
+-- Call 'error' if the vector is empty.
+fromVectorT :: VB.Vector (Link a) -> IO (Trace a)
+fromVectorT xs = Trace <$> C.fromVector xs
 
 -- | Get the length of the trace.
 lengthT :: Trace a -> Int
diff --git a/src/Mcmc/Cycle.hs b/src/Mcmc/Cycle.hs
--- a/src/Mcmc/Cycle.hs
+++ b/src/Mcmc/Cycle.hs
@@ -4,7 +4,7 @@
 -- |
 -- Module      :  Mcmc.Cycle
 -- Description :  A cycle is a list of proposals
--- Copyright   :  (c) 2021 Dominik Schrempf
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -15,7 +15,7 @@
 module Mcmc.Cycle
   ( -- * Cycles
     Order (..),
-    Cycle (ccProposals),
+    Cycle (ccProposals, ccRequireTrace),
     cycleFromList,
     setOrder,
     IterationMode (..),
@@ -31,14 +31,13 @@
 import qualified Data.ByteString.Builder as BB
 import qualified Data.ByteString.Lazy.Char8 as BL
 import Data.Default
-import Data.Either
 import Data.List
 import qualified Data.Map.Strict as M
 import qualified Data.Vector as VB
 import Mcmc.Acceptance
 import Mcmc.Internal.Shuffle
 import Mcmc.Proposal
-import System.Random.MWC
+import System.Random.Stateful
 
 -- | Define the order in which 'Proposal's are executed in a 'Cycle'. The total
 -- number of 'Proposal's per 'Cycle' may differ between 'Order's (e.g., compare
@@ -94,7 +93,9 @@
 -- that they can be uniquely identified.
 data Cycle a = Cycle
   { ccProposals :: [Proposal a],
-    ccOrder :: Order
+    ccOrder :: Order,
+    -- | Does the cycle require the trace when auto tuning? See 'tRequireTrace'.
+    ccRequireTrace :: Bool
   }
 
 -- | Create a 'Cycle' from a list of 'Proposal's.
@@ -103,7 +104,7 @@
   error "cycleFromList: Received an empty list but cannot create an empty Cycle."
 cycleFromList xs =
   if length uniqueXs == length xs
-    then Cycle xs def
+    then Cycle xs def (any needsTrace xs)
     else error $ "\n" ++ msg ++ "cycleFromList: Proposals are not unique."
   where
     uniqueXs = nub xs
@@ -112,6 +113,9 @@
     removedDescriptions = map (show . prDescription) removedXs
     removedMsgs = zipWith (\n d -> n ++ " " ++ d) removedNames removedDescriptions
     msg = unlines removedMsgs
+    needsTrace p = case prTuner p of
+      Nothing -> False
+      Just t -> tRequireTrace t
 
 -- | Set the order of 'Proposal's in a 'Cycle'.
 setOrder :: Order -> Cycle a -> Cycle a
@@ -122,14 +126,17 @@
   deriving (Eq)
 
 -- | Replicate 'Proposal's according to their weights and possibly shuffle them.
-prepareProposals :: IterationMode -> Cycle a -> GenIO -> IO [Proposal a]
-prepareProposals m (Cycle xs o) g = case o of
-  RandomO -> shuffle ps g
-  SequentialO -> return ps
-  RandomReversibleO -> do
-    psR <- shuffle ps g
-    return $ psR ++ reverse psR
-  SequentialReversibleO -> return $ ps ++ reverse ps
+prepareProposals :: StatefulGen g m => IterationMode -> Cycle a -> g -> m [Proposal a]
+prepareProposals m (Cycle xs o _) g =
+  if null ps
+    then error "prepareProposals: No proposals found."
+    else case o of
+      RandomO -> shuffle ps g
+      SequentialO -> return ps
+      RandomReversibleO -> do
+        psR <- shuffle ps g
+        return $ psR ++ reverse psR
+      SequentialReversibleO -> return $ ps ++ reverse ps
   where
     !ps =
       concat
@@ -143,7 +150,7 @@
 
 -- The number of proposals depends on the order.
 getNProposalsPerCycle :: IterationMode -> Cycle a -> Int
-getNProposalsPerCycle m (Cycle xs o) = case o of
+getNProposalsPerCycle m (Cycle xs o _) = case o of
   RandomO -> once
   SequentialO -> once
   RandomReversibleO -> 2 * once
@@ -155,28 +162,34 @@
     once = sum $ map (fromPWeight . prWeight) xs'
 
 -- See 'tuneWithTuningParameters' and 'Tuner'.
-tuneWithChainParameters :: AcceptanceRate -> VB.Vector a -> Proposal a -> Either String (Proposal a)
-tuneWithChainParameters ar xs p = case prTuner p of
-  Nothing -> Left "tuneWithChainParameters: Proposal is not tunable."
-  Just (Tuner t fT ts fTs _) ->
-    let t' = fT ar t
-        ts' = fTs xs ts
-     in tuneWithTuningParameters t' ts' p
+tuneWithChainParameters :: TuningType -> AcceptanceRate -> Maybe (VB.Vector a) -> Proposal a -> Either String (Proposal a)
+tuneWithChainParameters b ar mxs p = case prTuner p of
+  Nothing -> Right p
+  Just (Tuner t ts rt fT _) -> case (rt, mxs) of
+    (True, Nothing) -> error "tuneWithChainParameters: trace required"
+    _ ->
+      let (t', ts') = fT b d ar mxs (t, ts)
+       in tuneWithTuningParameters t' ts' p
+      where
+        d = prDimension p
 
 -- | Calculate acceptance rates and auto tunes the 'Proposal's in the 'Cycle'.
 --
 -- Do not change 'Proposal's that are not tuneable.
-autoTuneCycle :: Acceptance (Proposal a) -> VB.Vector a -> Cycle a -> Cycle a
-autoTuneCycle a xs c =
-  if sort (M.keys ar) == sort ps
-    then c {ccProposals = map tuneF ps}
-    else error "autoTuneCycle: Proposals in map and cycle do not match."
-  where
-    ar = acceptanceRates a
-    ps = ccProposals c
-    tuneF p = case ar M.!? p of
-      Just (Just x) -> fromRight p (tuneWithChainParameters x xs p)
-      _ -> p
+autoTuneCycle :: TuningType -> Acceptance (Proposal a) -> Maybe (VB.Vector a) -> Cycle a -> Cycle a
+autoTuneCycle b a mxs c = case (ccRequireTrace c, mxs) of
+  (False, Just _) -> error "autoTuneCycle: trace not required"
+  (True, Nothing) -> error "autoTuneCycle: trace required"
+  _ ->
+    if sort (M.keys ar) == sort ps
+      then c {ccProposals = map tuneF ps}
+      else error "autoTuneCycle: Proposals in map and cycle do not match."
+    where
+      ar = acceptanceRates a
+      ps = ccProposals c
+      tuneF p = case ar M.!? p of
+        Just (Just x) -> either error id $ tuneWithChainParameters b x mxs p
+        _ -> p
 
 -- | Horizontal line of proposal summaries.
 proposalHLine :: BL.ByteString
diff --git a/src/Mcmc/Environment.hs b/src/Mcmc/Environment.hs
--- a/src/Mcmc/Environment.hs
+++ b/src/Mcmc/Environment.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Environment
 -- Description :  Runtime environment
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Internal/ByteString.hs b/src/Mcmc/Internal/ByteString.hs
--- a/src/Mcmc/Internal/ByteString.hs
+++ b/src/Mcmc/Internal/ByteString.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Internal.ByteString
 -- Description :  ByteString tools
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Internal/Gamma.hs b/src/Mcmc/Internal/Gamma.hs
deleted file mode 100644
--- a/src/Mcmc/Internal/Gamma.hs
+++ /dev/null
@@ -1,205 +0,0 @@
--- |
--- Module      :  Mcmc.Internal.Gamma
--- Description :  Generalized gamma function for automatic differentiation
--- Copyright   :  (c) 2021 Dominik Schrempf
--- License     :  GPL-3.0-or-later
---
--- Maintainer  :  dominik.schrempf@gmail.com
--- Stability   :  experimental
--- Portability :  portable
---
--- Creation date: Tue Jul 13 12:53:09 2021.
---
--- The code is taken from "Numeric.SpecFunctions".
-module Mcmc.Internal.Gamma
-  ( logGammaG,
-  )
-where
-
-import Data.Typeable
-import qualified Data.Vector as VB
-import Numeric.Polynomial
-import Numeric.SpecFunctions
-import Unsafe.Coerce
-
-mSqrtEpsG :: RealFloat a => a
-mSqrtEpsG = 1.4901161193847656e-8
-
-mEulerMascheroniG :: RealFloat a => a
-mEulerMascheroniG = 0.5772156649015328606065121
-
--- | Generalized version of the log gamma distribution. See
--- 'Numeric.SpecFunctions.logGamma'.
-logGammaG :: (Typeable a, RealFloat a) => a -> a
-logGammaG z
-  | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z
-  | otherwise = logGammaNonDouble z
-{-# SPECIALIZE logGammaG :: Double -> Double #-}
-
--- | See 'Numeric.SpecFunctions.logGamma'.
-logGammaNonDouble :: RealFloat a => a -> a
-logGammaNonDouble z
-  | z <= 0 = 1 / 0
-  | z < mSqrtEpsG = log (1 / z - mEulerMascheroniG)
-  | z < 0.5 = lgamma1_15G z (z - 1) - log z
-  | z < 1 = lgamma15_2G z (z - 1) - log z
-  | z <= 1.5 = lgamma1_15G (z - 1) (z - 2)
-  | z < 2 = lgamma15_2G (z - 1) (z - 2)
-  | z < 15 = lgammaSmallG z
-  | otherwise = lanczosApproxG z
-
-lgamma1_15G :: RealFloat a => a -> a -> a
-lgamma1_15G zm1 zm2 =
-  r * y + r
-    * ( evaluatePolynomial zm1 tableLogGamma_1_15PG
-          / evaluatePolynomial zm1 tableLogGamma_1_15QG
-      )
-  where
-    r = zm1 * zm2
-    y = 0.52815341949462890625
-
-tableLogGamma_1_15PG :: RealFloat a => VB.Vector a
-tableLogGamma_1_15PG =
-  VB.fromList
-    [ 0.490622454069039543534e-1,
-      -0.969117530159521214579e-1,
-      -0.414983358359495381969e0,
-      -0.406567124211938417342e0,
-      -0.158413586390692192217e0,
-      -0.240149820648571559892e-1,
-      -0.100346687696279557415e-2
-    ]
-{-# NOINLINE tableLogGamma_1_15PG #-}
-
-tableLogGamma_1_15QG :: RealFloat a => VB.Vector a
-tableLogGamma_1_15QG =
-  VB.fromList
-    [ 1,
-      0.302349829846463038743e1,
-      0.348739585360723852576e1,
-      0.191415588274426679201e1,
-      0.507137738614363510846e0,
-      0.577039722690451849648e-1,
-      0.195768102601107189171e-2
-    ]
-{-# NOINLINE tableLogGamma_1_15QG #-}
-
-lgamma15_2G :: RealFloat a => a -> a -> a
-lgamma15_2G zm1 zm2 =
-  r * y + r
-    * ( evaluatePolynomial (-zm2) tableLogGamma_15_2PG
-          / evaluatePolynomial (-zm2) tableLogGamma_15_2QG
-      )
-  where
-    r = zm1 * zm2
-    y = 0.452017307281494140625
-
-tableLogGamma_15_2PG :: RealFloat a => VB.Vector a
-tableLogGamma_15_2PG =
-  VB.fromList
-    [ -0.292329721830270012337e-1,
-      0.144216267757192309184e0,
-      -0.142440390738631274135e0,
-      0.542809694055053558157e-1,
-      -0.850535976868336437746e-2,
-      0.431171342679297331241e-3
-    ]
-{-# NOINLINE tableLogGamma_15_2PG #-}
-
-tableLogGamma_15_2QG :: RealFloat a => VB.Vector a
-tableLogGamma_15_2QG =
-  VB.fromList
-    [ 1,
-      -0.150169356054485044494e1,
-      0.846973248876495016101e0,
-      -0.220095151814995745555e0,
-      0.25582797155975869989e-1,
-      -0.100666795539143372762e-2,
-      -0.827193521891290553639e-6
-    ]
-{-# NOINLINE tableLogGamma_15_2QG #-}
-
-lgammaSmallG :: RealFloat a => a -> a
-lgammaSmallG = go 0
-  where
-    go acc z
-      | z < 3 = acc + lgamma2_3G z
-      | otherwise = go (acc + log zm1) zm1
-      where
-        zm1 = z - 1
-
-lgamma2_3G :: RealFloat a => a -> a
-lgamma2_3G z =
-  r * y + r
-    * ( evaluatePolynomial zm2 tableLogGamma_2_3PG
-          / evaluatePolynomial zm2 tableLogGamma_2_3QG
-      )
-  where
-    r = zm2 * (z + 1)
-    zm2 = z - 2
-    y = 0.158963680267333984375e0
-
-tableLogGamma_2_3PG :: RealFloat a => VB.Vector a
-tableLogGamma_2_3PG =
-  VB.fromList
-    [ -0.180355685678449379109e-1,
-      0.25126649619989678683e-1,
-      0.494103151567532234274e-1,
-      0.172491608709613993966e-1,
-      -0.259453563205438108893e-3,
-      -0.541009869215204396339e-3,
-      -0.324588649825948492091e-4
-    ]
-{-# NOINLINE tableLogGamma_2_3PG #-}
-
-tableLogGamma_2_3QG :: RealFloat a => VB.Vector a
-tableLogGamma_2_3QG =
-  VB.fromList
-    [ 1,
-      0.196202987197795200688e1,
-      0.148019669424231326694e1,
-      0.541391432071720958364e0,
-      0.988504251128010129477e-1,
-      0.82130967464889339326e-2,
-      0.224936291922115757597e-3,
-      -0.223352763208617092964e-6
-    ]
-{-# NOINLINE tableLogGamma_2_3QG #-}
-
-lanczosApproxG :: RealFloat a => a -> a
-lanczosApproxG z =
-  (log (z + g - 0.5) - 1) * (z - 0.5)
-    + log (evalRatioG tableLanczosG z)
-  where
-    g = 6.024680040776729583740234375
-
-tableLanczosG :: RealFloat a => VB.Vector (a, a)
-tableLanczosG =
-  VB.fromList
-    [ (56906521.91347156388090791033559122686859, 0),
-      (103794043.1163445451906271053616070238554, 39916800),
-      (86363131.28813859145546927288977868422342, 120543840),
-      (43338889.32467613834773723740590533316085, 150917976),
-      (14605578.08768506808414169982791359218571, 105258076),
-      (3481712.15498064590882071018964774556468, 45995730),
-      (601859.6171681098786670226533699352302507, 13339535),
-      (75999.29304014542649875303443598909137092, 2637558),
-      (6955.999602515376140356310115515198987526, 357423),
-      (449.9445569063168119446858607650988409623, 32670),
-      (19.51992788247617482847860966235652136208, 1925),
-      (0.5098416655656676188125178644804694509993, 66),
-      (0.006061842346248906525783753964555936883222, 1)
-    ]
-{-# NOINLINE tableLanczosG #-}
-
-data LG a = LG !a !a
-
-evalRatioG :: RealFloat a => VB.Vector (a, a) -> a -> a
-evalRatioG coef x
-  | x > 1 = fini $ VB.foldl' stepL (LG 0 0) coef
-  | otherwise = fini $ VB.foldr' stepR (LG 0 0) coef
-  where
-    fini (LG num den) = num / den
-    stepR (a, b) (LG num den) = LG (num * x + a) (den * x + b)
-    stepL (LG num den) (a, b) = LG (num * rx + a) (den * rx + b)
-    rx = recip x
diff --git a/src/Mcmc/Internal/Random.hs b/src/Mcmc/Internal/Random.hs
--- a/src/Mcmc/Internal/Random.hs
+++ b/src/Mcmc/Internal/Random.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Internal.Random
 -- Description :  Tools for random calculations
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -12,37 +12,23 @@
 --
 -- Creation date: Wed Nov 25 07:14:52 2020.
 module Mcmc.Internal.Random
-  ( splitGen,
-    saveGen,
+  ( saveGen,
     loadGen,
   )
 where
 
-import Control.Monad
-import Control.Monad.Primitive
-import qualified Data.Vector.Unboxed as V
+import Data.IORef
 import Data.Word
-import System.Random.MWC
-
--- | Split a generator.
---
--- Splitting an MWC pseudo number generator is not good practice. However, I
--- have to go with this solution for now, and wait for proper support of
--- spittable pseudo random number generators such as @splitmix@.
-splitGen :: PrimMonad m => Int -> Gen (PrimState m) -> m [Gen (PrimState m)]
-splitGen n gen
-  | n <= 0 = return []
-  | otherwise = do
-      seeds :: [V.Vector Word32] <- replicateM n $ uniformVector gen 256
-      mapM initialize seeds
-
--- TODO: Splitmix. Remove or amend these functions as soon as split mix is used
--- and is available with the statistics package.
+import System.Random.Internal
+import System.Random.SplitMix
+import System.Random.Stateful
 
 -- | Save a generator to a seed.
-saveGen :: GenIO -> IO (V.Vector Word32)
-saveGen = fmap fromSeed . save
+saveGen :: IOGenM StdGen -> IO (Word64, Word64)
+saveGen (IOGenM r) = do
+  (StdGen g) <- readIORef r
+  pure $ unseedSMGen g
 
 -- | Load a generator from a seed.
-loadGen :: V.Vector Word32 -> IO GenIO
-loadGen = restore . toSeed
+loadGen :: (Word64, Word64) -> IO (IOGenM StdGen)
+loadGen s = newIOGenM $ StdGen $ seedSMGen' s
diff --git a/src/Mcmc/Internal/Shuffle.hs b/src/Mcmc/Internal/Shuffle.hs
--- a/src/Mcmc/Internal/Shuffle.hs
+++ b/src/Mcmc/Internal/Shuffle.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Internal.Shuffle
 -- Description :  Shuffle a list
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -21,23 +21,26 @@
 import Control.Monad.ST
 import qualified Data.Vector as V
 import qualified Data.Vector.Mutable as M
-import System.Random.MWC
+import System.Random.Stateful
 
 -- Fisher-Yates shuffle. See also
 -- 'System.Random.MWC.Distributions.uniformPermutation' which is a little
 -- cleaner, in my opinion. However, I would like to move away from MWC so I
 -- leave the custom implementation for now.
+--
+-- NOTE: Moving away from MWC is not anymore true because I can use SplitMix
+-- with it.
 
 -- | Shuffle a vector.
-shuffle :: [a] -> GenIO -> IO [a]
+shuffle :: StatefulGen g m => [a] -> g -> m [a]
 shuffle xs = grabble xs (length xs)
 
 -- | @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 -> GenIO -> IO [a]
+grabble :: StatefulGen g m => [a] -> Int -> g -> m [a]
 grabble xs m gen = do
   swaps <- forM [0 .. min (l - 1) m] $ \i -> do
-    j <- uniformR (i, l) gen
+    j <- uniformRM (i, l) gen
     return (i, j)
   return $ (V.toList . V.take m . swapElems (V.fromList xs)) swaps
   where
diff --git a/src/Mcmc/Internal/SpecFunctions.hs b/src/Mcmc/Internal/SpecFunctions.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Internal/SpecFunctions.hs
@@ -0,0 +1,410 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- |
+-- Module      :  Mcmc.Internal.Gamma
+-- Description :  Generalized gamma function for automatic differentiation
+-- Copyright   :  2021 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- Creation date: Tue Jul 13 12:53:09 2021.
+--
+-- The code is taken from "Numeric.SpecFunctions".
+module Mcmc.Internal.SpecFunctions
+  ( logGammaG,
+    logFactorialG,
+  )
+where
+
+import Data.Typeable
+import qualified Data.Vector as VB
+import Numeric.Polynomial
+import Numeric.SpecFunctions
+import Unsafe.Coerce
+
+mSqrtEpsG :: RealFloat a => a
+mSqrtEpsG = 1.4901161193847656e-8
+
+mEulerMascheroniG :: RealFloat a => a
+mEulerMascheroniG = 0.5772156649015328606065121
+
+-- | Generalized version of the log gamma distribution. See
+-- 'Numeric.SpecFunctions.logGamma'.
+logGammaG :: (Typeable a, RealFloat a) => a -> a
+logGammaG z
+  | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z
+  | otherwise = logGammaNonDouble z
+{-# SPECIALIZE logGammaG :: Double -> Double #-}
+
+-- | See 'Numeric.SpecFunctions.logGamma'.
+logGammaNonDouble :: RealFloat a => a -> a
+logGammaNonDouble z
+  | z <= 0 = 1 / 0
+  | z < mSqrtEpsG = log (1 / z - mEulerMascheroniG)
+  | z < 0.5 = lgamma1_15G z (z - 1) - log z
+  | z < 1 = lgamma15_2G z (z - 1) - log z
+  | z <= 1.5 = lgamma1_15G (z - 1) (z - 2)
+  | z < 2 = lgamma15_2G (z - 1) (z - 2)
+  | z < 15 = lgammaSmallG z
+  | otherwise = lanczosApproxG z
+
+lgamma1_15G :: RealFloat a => a -> a -> a
+lgamma1_15G zm1 zm2 =
+  r * y + r
+    * ( evaluatePolynomial zm1 tableLogGamma_1_15PG
+          / evaluatePolynomial zm1 tableLogGamma_1_15QG
+      )
+  where
+    r = zm1 * zm2
+    y = 0.52815341949462890625
+
+tableLogGamma_1_15PG :: RealFloat a => VB.Vector a
+tableLogGamma_1_15PG =
+  VB.fromList
+    [ 0.490622454069039543534e-1,
+      -0.969117530159521214579e-1,
+      -0.414983358359495381969e0,
+      -0.406567124211938417342e0,
+      -0.158413586390692192217e0,
+      -0.240149820648571559892e-1,
+      -0.100346687696279557415e-2
+    ]
+{-# NOINLINE tableLogGamma_1_15PG #-}
+
+tableLogGamma_1_15QG :: RealFloat a => VB.Vector a
+tableLogGamma_1_15QG =
+  VB.fromList
+    [ 1,
+      0.302349829846463038743e1,
+      0.348739585360723852576e1,
+      0.191415588274426679201e1,
+      0.507137738614363510846e0,
+      0.577039722690451849648e-1,
+      0.195768102601107189171e-2
+    ]
+{-# NOINLINE tableLogGamma_1_15QG #-}
+
+lgamma15_2G :: RealFloat a => a -> a -> a
+lgamma15_2G zm1 zm2 =
+  r * y + r
+    * ( evaluatePolynomial (-zm2) tableLogGamma_15_2PG
+          / evaluatePolynomial (-zm2) tableLogGamma_15_2QG
+      )
+  where
+    r = zm1 * zm2
+    y = 0.452017307281494140625
+
+tableLogGamma_15_2PG :: RealFloat a => VB.Vector a
+tableLogGamma_15_2PG =
+  VB.fromList
+    [ -0.292329721830270012337e-1,
+      0.144216267757192309184e0,
+      -0.142440390738631274135e0,
+      0.542809694055053558157e-1,
+      -0.850535976868336437746e-2,
+      0.431171342679297331241e-3
+    ]
+{-# NOINLINE tableLogGamma_15_2PG #-}
+
+tableLogGamma_15_2QG :: RealFloat a => VB.Vector a
+tableLogGamma_15_2QG =
+  VB.fromList
+    [ 1,
+      -0.150169356054485044494e1,
+      0.846973248876495016101e0,
+      -0.220095151814995745555e0,
+      0.25582797155975869989e-1,
+      -0.100666795539143372762e-2,
+      -0.827193521891290553639e-6
+    ]
+{-# NOINLINE tableLogGamma_15_2QG #-}
+
+lgammaSmallG :: RealFloat a => a -> a
+lgammaSmallG = go 0
+  where
+    go acc z
+      | z < 3 = acc + lgamma2_3G z
+      | otherwise = go (acc + log zm1) zm1
+      where
+        zm1 = z - 1
+
+lgamma2_3G :: RealFloat a => a -> a
+lgamma2_3G z =
+  r * y + r
+    * ( evaluatePolynomial zm2 tableLogGamma_2_3PG
+          / evaluatePolynomial zm2 tableLogGamma_2_3QG
+      )
+  where
+    r = zm2 * (z + 1)
+    zm2 = z - 2
+    y = 0.158963680267333984375e0
+
+tableLogGamma_2_3PG :: RealFloat a => VB.Vector a
+tableLogGamma_2_3PG =
+  VB.fromList
+    [ -0.180355685678449379109e-1,
+      0.25126649619989678683e-1,
+      0.494103151567532234274e-1,
+      0.172491608709613993966e-1,
+      -0.259453563205438108893e-3,
+      -0.541009869215204396339e-3,
+      -0.324588649825948492091e-4
+    ]
+{-# NOINLINE tableLogGamma_2_3PG #-}
+
+tableLogGamma_2_3QG :: RealFloat a => VB.Vector a
+tableLogGamma_2_3QG =
+  VB.fromList
+    [ 1,
+      0.196202987197795200688e1,
+      0.148019669424231326694e1,
+      0.541391432071720958364e0,
+      0.988504251128010129477e-1,
+      0.82130967464889339326e-2,
+      0.224936291922115757597e-3,
+      -0.223352763208617092964e-6
+    ]
+{-# NOINLINE tableLogGamma_2_3QG #-}
+
+lanczosApproxG :: RealFloat a => a -> a
+lanczosApproxG z =
+  (log (z + g - 0.5) - 1) * (z - 0.5)
+    + log (evalRatioG tableLanczosG z)
+  where
+    g = 6.024680040776729583740234375
+
+tableLanczosG :: RealFloat a => VB.Vector (a, a)
+tableLanczosG =
+  VB.fromList
+    [ (56906521.91347156388090791033559122686859, 0),
+      (103794043.1163445451906271053616070238554, 39916800),
+      (86363131.28813859145546927288977868422342, 120543840),
+      (43338889.32467613834773723740590533316085, 150917976),
+      (14605578.08768506808414169982791359218571, 105258076),
+      (3481712.15498064590882071018964774556468, 45995730),
+      (601859.6171681098786670226533699352302507, 13339535),
+      (75999.29304014542649875303443598909137092, 2637558),
+      (6955.999602515376140356310115515198987526, 357423),
+      (449.9445569063168119446858607650988409623, 32670),
+      (19.51992788247617482847860966235652136208, 1925),
+      (0.5098416655656676188125178644804694509993, 66),
+      (0.006061842346248906525783753964555936883222, 1)
+    ]
+{-# NOINLINE tableLanczosG #-}
+
+data LG a = LG !a !a
+
+evalRatioG :: RealFloat a => VB.Vector (a, a) -> a -> a
+evalRatioG coef x
+  | x > 1 = fini $ VB.foldl' stepL (LG 0 0) coef
+  | otherwise = fini $ VB.foldr' stepR (LG 0 0) coef
+  where
+    fini (LG num den) = num / den
+    stepR (a, b) (LG num den) = LG (num * x + a) (den * x + b)
+    stepL (LG num den) (a, b) = LG (num * rx + a) (den * rx + b)
+    rx = recip x
+
+-- | Generalized version of the log factorial function. See
+-- 'Numeric.SpecFunctions.logFactorial'.
+logFactorialG :: forall a b. (Integral a, RealFloat b, Typeable b) => a -> b
+logFactorialG n
+  | typeOf (undefined :: b) == typeOf (0 :: Double) = unsafeCoerce $ logFactorial n
+  | otherwise = logFactorialNonDouble n
+{-# SPECIALIZE logFactorialG :: Int -> Double #-}
+
+logFactorialNonDouble :: (Integral a, RealFloat b) => a -> b
+logFactorialNonDouble n
+  | n < 0 = error "logFactorialNonDouble: Negative input."
+  | n <= 170 = log $ VB.unsafeIndex factorialTable (fromIntegral n)
+  | n < 1500 = stirling + rx * ((1 / 12) - (1 / 360) * rx * rx)
+  | otherwise = stirling + (1 / 12) * rx
+  where
+    stirling = (x - 0.5) * log x - x + mLnSqrt2Pi
+    x = fromIntegral n + 1
+    rx = recip x
+{-# SPECIALIZE logFactorialNonDouble :: RealFloat a => Int -> a #-}
+
+mLnSqrt2Pi :: RealFloat a => a
+mLnSqrt2Pi = 0.9189385332046727417803297364056176398613974736377834128171
+{-# INLINE mLnSqrt2Pi #-}
+
+factorialTable :: RealFloat a => VB.Vector a
+{-# NOINLINE factorialTable #-}
+factorialTable =
+  VB.fromListN
+    171
+    [ 1.0,
+      1.0,
+      2.0,
+      6.0,
+      24.0,
+      120.0,
+      720.0,
+      5040.0,
+      40320.0,
+      362880.0,
+      3628800.0,
+      3.99168e7,
+      4.790016e8,
+      6.2270208e9,
+      8.71782912e10,
+      1.307674368e12,
+      2.0922789888e13,
+      3.55687428096e14,
+      6.402373705728e15,
+      1.21645100408832e17,
+      2.43290200817664e18,
+      5.109094217170944e19,
+      1.1240007277776077e21,
+      2.5852016738884974e22,
+      6.204484017332394e23,
+      1.5511210043330984e25,
+      4.032914611266056e26,
+      1.0888869450418352e28,
+      3.0488834461171384e29,
+      8.841761993739702e30,
+      2.6525285981219103e32,
+      8.222838654177922e33,
+      2.631308369336935e35,
+      8.683317618811886e36,
+      2.9523279903960412e38,
+      1.0333147966386144e40,
+      3.719933267899012e41,
+      1.3763753091226343e43,
+      5.23022617466601e44,
+      2.0397882081197442e46,
+      8.159152832478977e47,
+      3.3452526613163803e49,
+      1.4050061177528798e51,
+      6.041526306337383e52,
+      2.6582715747884485e54,
+      1.1962222086548019e56,
+      5.5026221598120885e57,
+      2.5862324151116818e59,
+      1.2413915592536073e61,
+      6.082818640342675e62,
+      3.0414093201713376e64,
+      1.5511187532873822e66,
+      8.065817517094388e67,
+      4.2748832840600255e69,
+      2.308436973392414e71,
+      1.2696403353658275e73,
+      7.109985878048634e74,
+      4.0526919504877214e76,
+      2.3505613312828785e78,
+      1.386831185456898e80,
+      8.32098711274139e81,
+      5.075802138772247e83,
+      3.146997326038793e85,
+      1.9826083154044399e87,
+      1.2688693218588415e89,
+      8.24765059208247e90,
+      5.44344939077443e92,
+      3.647111091818868e94,
+      2.4800355424368305e96,
+      1.711224524281413e98,
+      1.197857166996989e100,
+      8.504785885678623e101,
+      6.1234458376886085e103,
+      4.470115461512684e105,
+      3.307885441519386e107,
+      2.4809140811395396e109,
+      1.88549470166605e111,
+      1.4518309202828586e113,
+      1.1324281178206297e115,
+      8.946182130782974e116,
+      7.15694570462638e118,
+      5.797126020747368e120,
+      4.753643337012841e122,
+      3.9455239697206583e124,
+      3.314240134565353e126,
+      2.81710411438055e128,
+      2.422709538367273e130,
+      2.1077572983795275e132,
+      1.8548264225739844e134,
+      1.650795516090846e136,
+      1.4857159644817613e138,
+      1.352001527678403e140,
+      1.2438414054641305e142,
+      1.1567725070816416e144,
+      1.087366156656743e146,
+      1.0329978488239058e148,
+      9.916779348709496e149,
+      9.619275968248211e151,
+      9.426890448883246e153,
+      9.332621544394413e155,
+      9.332621544394415e157,
+      9.425947759838358e159,
+      9.614466715035125e161,
+      9.902900716486179e163,
+      1.0299016745145626e166,
+      1.0813967582402908e168,
+      1.1462805637347082e170,
+      1.2265202031961378e172,
+      1.3246418194518288e174,
+      1.4438595832024934e176,
+      1.5882455415227428e178,
+      1.7629525510902446e180,
+      1.974506857221074e182,
+      2.2311927486598134e184,
+      2.543559733472187e186,
+      2.9250936934930154e188,
+      3.393108684451898e190,
+      3.9699371608087206e192,
+      4.68452584975429e194,
+      5.574585761207606e196,
+      6.689502913449126e198,
+      8.094298525273443e200,
+      9.875044200833601e202,
+      1.214630436702533e205,
+      1.5061417415111406e207,
+      1.8826771768889257e209,
+      2.372173242880047e211,
+      3.0126600184576594e213,
+      3.856204823625804e215,
+      4.974504222477286e217,
+      6.466855489220473e219,
+      8.471580690878819e221,
+      1.1182486511960041e224,
+      1.4872707060906857e226,
+      1.9929427461615188e228,
+      2.6904727073180504e230,
+      3.6590428819525483e232,
+      5.012888748274991e234,
+      6.917786472619488e236,
+      9.615723196941088e238,
+      1.3462012475717523e241,
+      1.898143759076171e243,
+      2.6953641378881624e245,
+      3.8543707171800725e247,
+      5.5502938327393044e249,
+      8.047926057471992e251,
+      1.1749972043909107e254,
+      1.7272458904546386e256,
+      2.5563239178728654e258,
+      3.808922637630569e260,
+      5.713383956445854e262,
+      8.62720977423324e264,
+      1.3113358856834524e267,
+      2.0063439050956823e269,
+      3.0897696138473508e271,
+      4.789142901463393e273,
+      7.471062926282894e275,
+      1.1729568794264143e278,
+      1.8532718694937346e280,
+      2.946702272495038e282,
+      4.714723635992061e284,
+      7.590705053947218e286,
+      1.2296942187394494e289,
+      2.0044015765453023e291,
+      3.287218585534296e293,
+      5.423910666131589e295,
+      9.003691705778436e297,
+      1.5036165148649988e300,
+      2.526075744973198e302,
+      4.269068009004705e304,
+      7.257415615307998e306
+    ]
diff --git a/src/Mcmc/Jacobian.hs b/src/Mcmc/Jacobian.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Jacobian.hs
@@ -0,0 +1,32 @@
+-- |
+-- Module      :  Mcmc.Jacobian
+-- Description :  Types and convenience functions for Jacobians
+-- Copyright   :  2022 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- Creation date: Tue May 31 09:37:54 2022.
+module Mcmc.Jacobian
+  ( Jacobian,
+    JacobianG,
+    JacobianFunction,
+    JacobianFunctionG,
+  )
+where
+
+import Numeric.Log
+
+-- | Absolute value of the determinant of the Jacobian matrix.
+type Jacobian = Log Double
+
+-- | Generalized Jacobian.
+type JacobianG a = Log a
+
+-- | Function calculating the 'Jacobian'.
+type JacobianFunction a = JacobianFunctionG a Double
+
+-- | Function calculating the 'Jacobian'.
+type JacobianFunctionG a b = a -> JacobianG b
diff --git a/src/Mcmc/Likelihood.hs b/src/Mcmc/Likelihood.hs
--- a/src/Mcmc/Likelihood.hs
+++ b/src/Mcmc/Likelihood.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Likelihood
 -- Description :  Types and convenience functions for computing likelihoods
--- Copyright   :  (c) 2021 Dominik Schrempf
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Logger.hs b/src/Mcmc/Logger.hs
--- a/src/Mcmc/Logger.hs
+++ b/src/Mcmc/Logger.hs
@@ -4,7 +4,7 @@
 -- |
 -- Module      :  Mcmc.Logger
 -- Description :  Minimal monad logger
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/MarginalLikelihood.hs b/src/Mcmc/MarginalLikelihood.hs
--- a/src/Mcmc/MarginalLikelihood.hs
+++ b/src/Mcmc/MarginalLikelihood.hs
@@ -4,7 +4,7 @@
 -- |
 -- Module      :  Mcmc.MarginalLikelihood
 -- Description :  Calculate the marginal likelihood
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -38,7 +38,6 @@
 import Mcmc.Chain.Trace
 import Mcmc.Cycle
 import Mcmc.Environment
-import Mcmc.Internal.Random
 import Mcmc.Likelihood
 import Mcmc.Logger
 import Mcmc.Mcmc
@@ -47,7 +46,7 @@
 import Mcmc.Settings
 import Numeric.Log hiding (sum)
 import System.Directory
-import System.Random.MWC
+import System.Random.Stateful
 import Text.Printf
 import Text.Show.Pretty
 import Prelude hiding (cycle)
@@ -224,7 +223,7 @@
   Cycle a ->
   Monitor a ->
   a ->
-  GenIO ->
+  IOGenM StdGen ->
   -- For each point a vector of likelihoods stored in log domain.
   ML [VU.Vector Likelihood]
 mlRun k xs em vb prf lhf cc mn i0 g = do
@@ -271,11 +270,12 @@
   Cycle a ->
   Monitor a ->
   a ->
-  GenIO ->
+  IOGenM StdGen ->
   ML MarginalLikelihood
-tiWrapper s prf lhf cc mn i0 g = do
+tiWrapper s prf lhf cc mn i0 g0 = do
   logInfoB "Path integral (thermodynamic integration)."
-  [g0, g1] <- splitGen 2 g
+  r1 <- splitGenM g0
+  g1 <- newIOGenM r1
 
   -- Parallel execution of both path integrals.
   r <- ask
@@ -354,7 +354,7 @@
   Cycle a ->
   Monitor a ->
   a ->
-  GenIO ->
+  IOGenM StdGen ->
   ML MarginalLikelihood
 sssWrapper s prf lhf cc mn i0 g = do
   logInfoB "Stepping stone sampling."
@@ -378,9 +378,7 @@
   Cycle a ->
   Monitor a ->
   InitialState a ->
-  -- | A source of randomness. For reproducible runs, make sure to use
-  -- generators with the same seed.
-  GenIO ->
+  IOGenM StdGen ->
   IO MarginalLikelihood
 marginalLikelihood s prf lhf cc mn i0 g = do
   -- Initialize.
@@ -401,7 +399,7 @@
           ThermodynamicIntegration -> tiWrapper s prf lhf cc mn i0 g
           SteppingStoneSampling -> sssWrapper s prf lhf cc mn i0 g
         logInfoS $ "Marginal log likelihood: " ++ show (ln val)
-        -- TODO: Simulation variance.
+        -- TODO (low): Simulation variance.
         logInfoS "The simulation variance is not yet available."
         return val
     )
diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs
--- a/src/Mcmc/Mcmc.hs
+++ b/src/Mcmc/Mcmc.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Mcmc
 -- Description :  Framework for running Markov chain Monte Carlo samplers
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -29,6 +29,7 @@
 import Mcmc.Cycle
 import Mcmc.Environment
 import Mcmc.Logger
+import Mcmc.Proposal (TuningType (LastTuningStep, NormalTuningStep))
 import Mcmc.Settings
 import System.Exit
 import System.IO
@@ -55,25 +56,25 @@
   return $ aResetAcceptance a
 
 mcmcExceptionHandler :: Algorithm a => Environment Settings -> a -> AsyncException -> IO b
-mcmcExceptionHandler e a UserInterrupt = do
+mcmcExceptionHandler e a err = do
   putStrLn ""
-  putStrLn "USER INTERRUPT!"
+  putStrLn "INTERRUPT!"
   putStrLn "Try to terminate gracefully and save chain for continuation."
-  putStrLn "Press CTRL-C again to terminate now."
-  putStrLn "Close output files."
+  putStrLn "Press CTRL-C (again) to terminate now."
+  putStrLn "Closing output files."
   _ <- aCloseMonitors a
   closeEnvironment e
-  putStrLn "Try to save settings."
+  putStrLn "Saving settings."
   let s = settings e
   settingsSave s
-  putStrLn "Try to save compressed MCMC analysis."
+  putStrLn "Saving compressed MCMC analysis."
   putStrLn "For long traces, or complex objects, this may take a while."
   let nm = sAnalysisName s
   aSave nm a
   putStrLn "Markov chain saved. Analysis can be continued."
-  putStrLn "Terminate gracefully."
-  exitWith $ ExitFailure 1
-mcmcExceptionHandler _ _ e = throw e
+  putStrLn "Graceful termination successful."
+  putStrLn "Rethrowing error."
+  throw err
 
 -- XXX: Exception handling. Is it enough to mask execution of monitors and catch
 -- UserInterrupt during iterations?
@@ -110,10 +111,13 @@
   mcmcExecuteMonitors a
   when (aIsInValidState a) (logWarnB "The initial state is invalid!")
   a' <- mcmcBurnIn a
+  logInfoS $ "Clean chain after burn in."
+  let tl = sTraceLength s
+  a'' <- liftIO $ aCleanAfterBurnIn tl a'
   let i = fromIterations $ sIterations s
   logInfoS $ "Run chain for " ++ show i ++ " iterations."
-  logInfoB $ aStdMonitorHeader a'
-  mcmcIterate AllProposals i a'
+  logInfoB $ aStdMonitorHeader a''
+  mcmcIterate AllProposals i a''
 
 mcmcContinueRun :: Algorithm a => a -> MCMC a
 mcmcContinueRun a = do
@@ -183,23 +187,26 @@
       return a''
 
 -- Auto tune the proposals.
-mcmcAutotune :: Algorithm a => Int -> a -> MCMC a
-mcmcAutotune n a = do
+mcmcAutotune :: Algorithm a => TuningType -> Int -> a -> MCMC a
+mcmcAutotune NormalTuningStep n a = do
   logDebugB "Auto tune."
-  liftIO $ aAutoTune n a
+  liftIO $ aAutoTune NormalTuningStep n a
+mcmcAutotune LastTuningStep n a = do
+  logDebugB "Last auto tune."
+  liftIO $ aAutoTune LastTuningStep n a
 
 mcmcBurnInWithAutoTuning :: Algorithm a => IterationMode -> [Int] -> a -> MCMC a
 mcmcBurnInWithAutoTuning _ [] _ = error "mcmcBurnInWithAutoTuning: Empty list."
 mcmcBurnInWithAutoTuning m [x] a = do
   -- Last round.
   a' <- mcmcIterate m x a
-  a'' <- mcmcAutotune x a'
+  a'' <- mcmcAutotune LastTuningStep x a'
   logInfoB $ aSummarizeCycle m a''
   logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations."
   mcmcResetAcceptance a''
 mcmcBurnInWithAutoTuning m (x : xs) a = do
   a' <- mcmcIterate m x a
-  a'' <- mcmcAutotune x a'
+  a'' <- mcmcAutotune NormalTuningStep x a'
   logDebugB $ aSummarizeCycle m a''
   logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations."
   logDebugB $ aStdMonitorHeader a''
diff --git a/src/Mcmc/Monitor.hs b/src/Mcmc/Monitor.hs
--- a/src/Mcmc/Monitor.hs
+++ b/src/Mcmc/Monitor.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Monitor
 -- Description :  Monitor a Markov chain
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Monitor/Log.hs b/src/Mcmc/Monitor/Log.hs
--- a/src/Mcmc/Monitor/Log.hs
+++ b/src/Mcmc/Monitor/Log.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Monitor.Log
 -- Description :  Monitor logarithmic values
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
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
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Monitor.Parameter
 -- Description :  Monitor parameters
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Monitor/ParameterBatch.hs b/src/Mcmc/Monitor/ParameterBatch.hs
--- a/src/Mcmc/Monitor/ParameterBatch.hs
+++ b/src/Mcmc/Monitor/ParameterBatch.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Monitor.ParameterBatch
 -- Description :  Batch monitor parameters
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Monitor/Time.hs b/src/Mcmc/Monitor/Time.hs
--- a/src/Mcmc/Monitor/Time.hs
+++ b/src/Mcmc/Monitor/Time.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Monitor.Time
 -- Description :  Print time related values
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Posterior.hs b/src/Mcmc/Posterior.hs
--- a/src/Mcmc/Posterior.hs
+++ b/src/Mcmc/Posterior.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Posterior
 -- Description :  Types for posterior values and functions
--- Copyright   :  (c) 2021 Dominik Schrempf
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Prior.hs b/src/Mcmc/Prior.hs
--- a/src/Mcmc/Prior.hs
+++ b/src/Mcmc/Prior.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Prior
 -- Description :  Types and convenience functions for computing priors
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -46,12 +46,19 @@
 import Control.Monad
 import Data.Maybe (fromMaybe)
 import Data.Typeable
-import Mcmc.Internal.Gamma
+import Mcmc.Internal.SpecFunctions
 import Mcmc.Statistics.Types
 import Numeric.Log
-import qualified Statistics.Distribution as S
-import qualified Statistics.Distribution.Poisson as S
 
+-- TODO (high): Think about using a "structure" variable.
+--
+-- For example,
+--
+-- type PriorFunctionG s a = s a -> PriorG a
+--
+-- Many of the prior functions would need to use the Identity Functor. This may
+-- be slow.
+
 -- | Prior values are stored in log domain.
 type Prior = PriorG Double
 
@@ -90,7 +97,7 @@
 
 -- | Improper uniform prior; strictly less than zero.
 negative :: RealFloat a => PriorFunctionG a a
-negative = lessThan 0
+negative = lessThan 0.0
 {-# SPECIALIZE negative :: PriorFunction Double #-}
 
 -- | Exponential distributed prior.
@@ -98,8 +105,8 @@
 -- Call 'error' if the rate is zero or negative.
 exponential :: RealFloat a => Rate a -> PriorFunctionG a a
 exponential l x
-  | l <= 0 = error "exponential: Rate is zero or negative."
-  | x < 0 = error "exponential: Negative value."
+  | l <= 0.0 = error "exponential: Rate is zero or negative."
+  | x <= 0.0 = 0.0
   | otherwise = ll * Exp (negate l * x)
   where
     ll = Exp $ log l
@@ -110,11 +117,10 @@
 -- Call 'error' if the shape or scale are zero or negative.
 gamma :: (Typeable a, RealFloat a) => Shape a -> Scale a -> PriorFunctionG a a
 gamma k t x
-  | k <= 0 = error "gamma: Shape is zero or negative."
-  | t <= 0 = error "gamma: Scale is zero or negative."
-  | x < 0 = error "gamma: Negative value."
-  | x == 0 = 0.0
-  | otherwise = Exp $ log x * (k - 1) - (x / t) - logGammaG k - log t * k
+  | k <= 0.0 = error "gamma: Shape is zero or negative."
+  | t <= 0.0 = error "gamma: Scale is zero or negative."
+  | x <= 0.0 = 0.0
+  | otherwise = Exp $ log x * (k - 1.0) - (x / t) - logGammaG k - log t * k
 {-# SPECIALIZE gamma :: Double -> Double -> PriorFunction Double #-}
 
 -- | See 'gamma' but parametrized using mean and variance.
@@ -167,12 +173,12 @@
 -- Call 'error' if the standard deviation is zero or negative.
 logNormal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a
 logNormal m s x
-  | s <= 0 = error "logNormal: Standard deviation is zero or negative."
-  | x <= 0 = 0
+  | s <= 0.0 = error "logNormal: Standard deviation is zero or negative."
+  | x <= 0.0 = 0.0
   | otherwise = Exp $ t + e
   where
     t = negate $ mLnSqrt2Pi + log (x * s)
-    a = recip $ 2 * s * s
+    a = recip $ 2.0 * s * s
     b = log x - m
     e = negate $ a * b * b
 
@@ -202,13 +208,11 @@
 -- | Poisson distributed prior.
 --
 -- Call 'error' if the rate is zero or negative.
-poisson :: Rate Double -> PriorFunction Int
+poisson :: (RealFloat a, Typeable a) => Rate a -> PriorFunctionG Int a
 poisson l n
-  | l <= 0 = error "poisson: Rate is zero or negative."
-  | n < 0 = error "poisson: Negative value."
-  | otherwise = Exp $ S.logProbability d n
-  where
-    d = S.poisson l
+  | l < 0.0 = error "poisson: Rate is zero or negative."
+  | n < 0 = 0.0
+  | otherwise = Exp $ log l * fromIntegral n - logFactorialG n - l
 
 -- | Intelligent product that stops when encountering a zero.
 --
@@ -220,5 +224,5 @@
 
 -- The type could be generalized to any MonadPlus Integer
 prodM :: RealFloat a => [Log a] -> Maybe (Log a)
-prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0)) 1
+prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0.0)) 1.0
 {-# SPECIALIZE prodM :: [Log Double] -> Maybe (Log Double) #-}
diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs
--- a/src/Mcmc/Proposal.hs
+++ b/src/Mcmc/Proposal.hs
@@ -5,7 +5,7 @@
 -- |
 -- Module      :  Mcmc.Proposal
 -- Description :  Proposals are instruction to move around the state space
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -23,24 +23,25 @@
     PSpeed (..),
     Proposal (..),
     KernelRatio,
+    PResult (..),
     Jacobian,
     JacobianFunction,
     (@~),
     liftProposal,
     liftProposalWith,
-    ProposalSimple,
+    PFunction,
     createProposal,
 
     -- * Tuners
     Tuner (..),
     Tune (..),
     TuningParameter,
+    TuningType (..),
     TuningFunction,
     AuxiliaryTuningParameters,
-    AuxiliaryTuningFunction,
-    defaultTuningFunctionWith,
-    noTuningFunction,
-    noAuxiliaryTuningFunction,
+    tuningFunction,
+    tuningFunctionWithAux,
+    tuningFunctionOnlyAux,
     tuningParameterMin,
     tuningParameterMax,
     tuneWithTuningParameters,
@@ -52,6 +53,7 @@
   )
 where
 
+import Data.Bifunctor
 import qualified Data.ByteString.Builder as BB
 import qualified Data.ByteString.Lazy.Char8 as BL
 import qualified Data.Double.Conversion.ByteString as BC
@@ -61,8 +63,9 @@
 import Lens.Micro
 import Mcmc.Acceptance
 import Mcmc.Internal.ByteString
+import Mcmc.Jacobian
 import Numeric.Log hiding (sum)
-import System.Random.MWC
+import System.Random.Stateful
 
 -- | Proposal name.
 newtype PName = PName {fromPName :: String}
@@ -164,8 +167,8 @@
     -- | The weight determines how often a 'Proposal' is executed per iteration of
     -- the Markov chain.
     prWeight :: PWeight,
-    -- | Simple proposal without name, weight, and tuning information.
-    prSimple :: ProposalSimple a,
+    -- | Simple proposal function without name, weight, and tuning information.
+    prFunction :: PFunction a,
     -- | Tuning is disabled if set to 'Nothing'.
     prTuner :: Maybe (Tuner a)
   }
@@ -178,24 +181,29 @@
 
 -- | Ratio of the proposal kernels.
 --
--- Part of the MHG acceptance ratio.
---
--- See also 'Jacobian'.
+-- For unbiased, volume preserving proposals, the values is 1.0.
 --
--- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,
--- it is more declarative to have them separate. Like so, we are constantly
--- reminded: Is the Jacobian modifier different from 1.0?
+-- For biased proposals, the kernel ratio is qYX / qXY, where qAB is the
+-- probability density to move from A to B.
 type KernelRatio = Log Double
 
--- | Absolute value of the determinant of the Jacobian matrix.
---
--- Part of the MHG acceptance ratio.
---
--- See also 'Jacobian'.
-type Jacobian = Log Double
-
--- | Function calculating the 'Jacobian' of a proposal.
-type JacobianFunction a = a -> Jacobian
+-- | Proposal result.
+data PResult a
+  = -- | Accept the new value regardless of the prior, likelihood or Jacobian.
+    ForceAccept !a
+  | -- | Reject the proposal regardless of the prior, likelihood or Jacobian.
+    ForceReject
+  | -- | Propose a new value.
+    --
+    -- In order to calculate the Metropolis-Hastings-Green ratio, we need to know
+    -- the ratio of the backward to forward kernels (the 'KernelRatio' or the
+    -- probability masses or probability densities) and the 'Jacobian'.
+    --
+    -- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,
+    -- it is more declarative to have them separate. Like so, we are constantly
+    -- reminded: Is the Jacobian modifier different from 1.0?
+    Propose !a !KernelRatio !Jacobian
+  deriving (Show, Eq)
 
 -- | Lift a proposal from one data type to another.
 --
@@ -219,55 +227,52 @@
 
 -- | Lift a proposal from one data type to another.
 --
--- A function to calculate the Jacobian has to be provided (but see '(@~)').
+-- A function to calculate the Jacobian has to be provided. If the Jabobian is
+-- 1.0, use 'liftProposal'.
 --
 -- For further reference, please see the [example
 -- @Pair@](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs).
 liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b
 liftProposalWith jf l (Proposal n r d p w s t) =
-  Proposal n r d p w (liftProposalSimpleWith jf l s) (liftTunerWith jf l <$> t)
+  Proposal n r d p w (liftPFunctionWith jf l s) (liftTunerWith jf l <$> t)
 
--- | Simple proposal without tuning information.
+-- TODO @Dominik (high, feature): Proposals should be aware of: Is this Burn in, or not?
+
+-- | Simple proposal function without tuning information.
 --
 -- Instruction about randomly moving from the current state to a new state,
 -- given some source of randomness.
 --
--- In order to calculate the Metropolis-Hastings-Green ratio, we need to know
--- the ratio of the backward to forward kernels (the 'KernelRatio' or the
--- probability masses or probability densities) and the 'Jacobian'.
---
--- 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, KernelRatio, Jacobian)
+-- Maybe report acceptance counts internal to the proposal (e.g., used by
+-- proposals based on Hamiltonian dynamics).
+type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceCounts)
 
--- Lift a simple proposal from one data type to another.
-liftProposalSimpleWith :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b
-liftProposalSimpleWith jf l s = s'
+-- Lift a proposal function from one data type to another.
+liftPFunctionWith :: JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b
+liftPFunctionWith jf l s = s'
   where
     s' y g = do
-      (x', r, j) <- s (y ^. l) g
-      let y' = set l x' y
-          jxy = jf y
-          jyx = jf y'
-          j' = j * jyx / jxy
-      return (y', r, j')
+      (pr, ac) <- s (y ^. l) g
+      let pr' = case pr of
+            ForceAccept x' -> ForceAccept $ set l x' y
+            ForceReject -> ForceReject
+            Propose x' r j ->
+              let y' = set l x' y
+                  jxy = jf y
+                  jyx = jf y'
+                  j' = j * jyx / jxy
+               in Propose y' r j'
+      pure (pr', ac)
 
 -- | Create a proposal with a single tuning parameter.
 --
--- Proposals with arbitrary tuning parameters have to be created manually. See
+-- Proposals with auxiliary tuning parameters have to be created manually. See
 -- 'Tuner' for more information, and 'Mcmc.Proposal.Hamiltonian' for an example.
 createProposal ::
   -- | Description of the proposal type and parameters.
   PDescription ->
-  -- | Function creating a simple proposal for a given tuning parameter.
-  (TuningParameter -> ProposalSimple a) ->
+  -- | Function creating a simple proposal function for a given tuning parameter.
+  (TuningParameter -> PFunction a) ->
   -- | Speed.
   PSpeed ->
   -- | Dimension.
@@ -282,36 +287,36 @@
 createProposal r f s d n w Tune =
   Proposal n r s d w (f 1.0) (Just tuner)
   where
-    fT = defaultTuningFunctionWith d
-    fTs = noAuxiliaryTuningFunction
+    fT = tuningFunction
     g t _ = Right $ f t
-    tuner = Tuner 1.0 fT VU.empty fTs g
+    tuner = Tuner 1.0 VU.empty False fT g
 createProposal r f s d n w NoTune =
   Proposal n r s d w (f 1.0) Nothing
 
 -- | Required information to tune 'Proposal's.
 data Tuner a = Tuner
   { tTuningParameter :: TuningParameter,
-    tTuningFunction :: TuningFunction,
     tAuxiliaryTuningParameters :: AuxiliaryTuningParameters,
-    tAuxiliaryTuningFunction :: AuxiliaryTuningFunction a,
+    -- | Does the tuner require the trace over the last tuning period?
+    tRequireTrace :: Bool,
+    tTuningFunction :: TuningFunction a,
     -- | Given the tuning parameter, and the auxiliary tuning parameters, get
-    -- the tuned simple proposal.
+    -- the tuned propose function.
     --
     -- Should return 'Left' if the vector of auxiliary tuning parameters is
     -- invalid.
-    tGetSimpleProposal ::
+    tGetPFunction ::
       TuningParameter ->
       AuxiliaryTuningParameters ->
-      Either String (ProposalSimple a)
+      Either String (PFunction a)
   }
 
 -- Lift tuner from one data type to another.
 liftTunerWith :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b
-liftTunerWith jf l (Tuner p fP ps fPs g) = Tuner p fP ps fPs' g'
+liftTunerWith jf l (Tuner p ps nt fP g) = Tuner p ps nt fP' g'
   where
-    fPs' = fPs . VB.map (^. l)
-    g' x xs = liftProposalSimpleWith jf l <$> g x xs
+    fP' b d r = fP b d r . fmap (VB.map (^. l))
+    g' x xs = liftPFunctionWith jf l <$> g x xs
 
 -- | Tune proposal?
 data Tune = Tune | NoTune
@@ -323,38 +328,53 @@
 -- expected acceptance rate; and vice versa.
 type TuningParameter = Double
 
--- | Compute new tuning parameter from a given acceptance rate and the old
--- tuning parameter.
-type TuningFunction = AcceptanceRate -> TuningParameter -> TuningParameter
+-- | The last tuning step may be special.
+data TuningType = NormalTuningStep | LastTuningStep
 
--- | Auxiliary tuning parameters; vector may be empty.
+-- | Compute new tuning parameters.
+type TuningFunction a =
+  TuningType ->
+  PDimension ->
+  -- | Acceptance rate of last tuning period.
+  AcceptanceRate ->
+  -- | Trace of last tuning period. Only available when requested by proposal.
+  Maybe (VB.Vector a) ->
+  (TuningParameter, AuxiliaryTuningParameters) ->
+  (TuningParameter, AuxiliaryTuningParameters)
+
+-- | Auxiliary tuning parameters.
 --
 -- Auxiliary tuning parameters are not shown in proposal summaries.
+--
+-- Vector may be empty.
 type AuxiliaryTuningParameters = VU.Vector TuningParameter
 
--- | Compute new auxiliary tuning parameters from a given trace and the old
--- auxiliary tuning parameters.
-type AuxiliaryTuningFunction a = VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters
+tuningFunctionSimple :: PDimension -> AcceptanceRate -> TuningParameter -> TuningParameter
+tuningFunctionSimple d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t
 
 -- | Default tuning function.
 --
--- Subject to change.
-defaultTuningFunctionWith ::
-  -- Optimal acceptance rate.
-  PDimension ->
-  TuningFunction
-defaultTuningFunctionWith d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t
+-- The default tuning function only uses the acceptance rate. In particular, it
+-- does not handle auxiliary tuning parameters and ignores the actual samples
+-- attained during the last tuning period.
+tuningFunction :: TuningFunction a
+tuningFunction _ d r _ = bimap (tuningFunctionSimple d r) id
 
--- | Do not tune.
---
--- Useful if auxiliary tuning parameters are tuned, but not the main tuning
--- parameter.
-noTuningFunction :: TuningFunction
-noTuningFunction _ = id
+-- | Also tune auxiliary tuning parameters.
+tuningFunctionWithAux ::
+  -- | Auxiliary tuning function.
+  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->
+  TuningFunction a
+tuningFunctionWithAux _ _ _ _ Nothing = error "tuningFunctionWithAux: empty trace"
+tuningFunctionWithAux f b d r (Just xs) = bimap (tuningFunctionSimple d r) (f b xs)
 
--- | Do not tune auxiliary parameters.
-noAuxiliaryTuningFunction :: AuxiliaryTuningFunction a
-noAuxiliaryTuningFunction _ = id
+-- | Only tune auxiliary tuning parameters.
+tuningFunctionOnlyAux ::
+  -- | Auxiliary tuning function.
+  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->
+  TuningFunction a
+tuningFunctionOnlyAux _ _ _ _ Nothing = error "tuningFunctionOnlyAux: empty trace"
+tuningFunctionOnlyAux f b _ _ (Just xs) = bimap id (f b xs)
 
 -- IDEA: Per proposal type tuning parameter boundaries. For example, a sliding
 -- proposal with a large tuning parameter is not a problem. But then, if the
@@ -391,14 +411,14 @@
   Either String (Proposal a)
 tuneWithTuningParameters t ts p = case prTuner p of
   Nothing -> Left "tuneWithTuningParameters: Proposal is not tunable."
-  Just (Tuner _ fT _ fTs g) ->
+  Just (Tuner _ _ nt fT g) ->
     -- Ensure that the tuning parameter is strictly positive and well bounded.
     let t' = max tuningParameterMin t
         t'' = min tuningParameterMax t'
         psE = g t'' ts
      in case psE of
           Left err -> Left $ "tune: " <> err
-          Right ps -> Right $ p {prSimple = ps, prTuner = Just $ Tuner t'' fT ts fTs g}
+          Right ps -> Right $ p {prFunction = ps, prTuner = Just $ Tuner t'' ts nt fT g}
 
 -- | See 'PDimension'.
 getOptimalRate :: PDimension -> Double
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
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Proposal.Bactrian
 -- Description :  Bactrian proposals
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -23,8 +23,8 @@
 import Numeric.Log
 import Statistics.Distribution
 import Statistics.Distribution.Normal
-import System.Random.MWC
 import System.Random.MWC.Distributions
+import System.Random.Stateful
 
 -- | Type synonym indicating the spike parameter.
 type SpikeParameter = Double
@@ -32,7 +32,7 @@
 genBactrian ::
   SpikeParameter ->
   StandardDeviation Double ->
-  GenIO ->
+  IOGenM StdGen ->
   IO Double
 genBactrian m s g = do
   let mn = m * s
@@ -55,21 +55,20 @@
 bactrianAdditive ::
   SpikeParameter ->
   StandardDeviation Double ->
-  ProposalSimple Double
+  PFunction Double
 bactrianAdditive m s x g = do
   dx <- genBactrian m s g
-  return (x + dx, 1.0, 1.0)
+  pure (Propose (x + dx) 1.0 1.0, Nothing)
 
--- bactrianSimple lens spike stdDev tune forwardOp backwardOp
-bactrianAdditiveSimple ::
+bactrianAdditivePFunction ::
   SpikeParameter ->
   StandardDeviation Double ->
   TuningParameter ->
-  ProposalSimple Double
-bactrianAdditiveSimple m s t
-  | m < 0 = error "bactrianAdditiveSimple: Spike parameter negative."
-  | m >= 1 = error "bactrianAdditiveSimple: Spike parameter 1.0 or larger."
-  | s <= 0 = error "bactrianAdditiveSimple: Standard deviation 0.0 or smaller."
+  PFunction Double
+bactrianAdditivePFunction m s t
+  | m < 0 = error "bactrianAdditivePFunction: Spike parameter negative."
+  | m >= 1 = error "bactrianAdditivePFunction: Spike parameter 1.0 or larger."
+  | s <= 0 = error "bactrianAdditivePFunction: Standard deviation 0.0 or smaller."
   | otherwise = bactrianAdditive m (t * s)
 
 -- | Additive symmetric proposal with kernel similar to the silhouette of a
@@ -90,7 +89,7 @@
   PWeight ->
   Tune ->
   Proposal Double
-slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) PFast (PDimension 1)
+slideBactrian m s = createProposal description (bactrianAdditivePFunction m s) PFast (PDimension 1)
   where
     description = PDescription $ "Slide Bactrian; spike: " ++ show m ++ ", sd: " ++ show s
 
@@ -106,24 +105,24 @@
 bactrianMult ::
   SpikeParameter ->
   StandardDeviation Double ->
-  ProposalSimple Double
+  PFunction Double
 bactrianMult m s x g = do
   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)
+  pure (Propose (x * u) (qYX / qXY) jac, Nothing)
 
-bactrianMultSimple ::
+bactrianMultPFunction ::
   SpikeParameter ->
   StandardDeviation Double ->
   TuningParameter ->
-  ProposalSimple Double
-bactrianMultSimple m s t
-  | m < 0 = error "bactrianMultSimple: Spike parameter negative."
-  | m >= 1 = error "bactrianMultSimple: Spike parameter 1.0 or larger."
-  | s <= 0 = error "bactrianMultSimple: Standard deviation 0.0 or smaller."
+  PFunction Double
+bactrianMultPFunction m s t
+  | m < 0 = error "bactrianMultPFunction: Spike parameter negative."
+  | m >= 1 = error "bactrianMultPFunction: Spike parameter 1.0 or larger."
+  | s <= 0 = error "bactrianMultPFunction: Standard deviation 0.0 or smaller."
   | otherwise = bactrianMult m (t * s)
 
 -- | Multiplicative proposal with kernel similar to the silhouette of a Bactrian
@@ -137,6 +136,6 @@
   PWeight ->
   Tune ->
   Proposal Double
-scaleBactrian m s = createProposal description (bactrianMultSimple m s) PFast (PDimension 1)
+scaleBactrian m s = createProposal description (bactrianMultPFunction m s) PFast (PDimension 1)
   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
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Proposal.Generic
 -- Description :  Generic interface for creating proposals
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -28,7 +28,7 @@
 -- 1. Let \(D\) be a continuous probability distribution on \(\mathbb{D}\);
 --    sample an auxiliary variable \(epsilon \sim D\).
 --
--- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). Propose a
+-- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). PFunction a
 --    new state \(x' = x \odot \epsilon\).
 --
 -- 3. If the proposal is unbiased, the Metropolis-Hastings-Green ratio can
@@ -71,7 +71,7 @@
   -- Conversion to log domain is necessary, because some determinants of
   -- Jacobians are very small (or large).
   Maybe (a -> Double -> Jacobian) ->
-  ProposalSimple a
+  PFunction a
 genericContinuous d f mInv mJac x g = do
   u <- genContVar d g
   let r = case mInv of
@@ -83,7 +83,7 @@
       j = case mJac of
         Nothing -> 1.0
         Just fJac -> fJac x u
-  return (x `f` u, r, j)
+  pure (Propose (x `f` u) r j, Nothing)
 {-# INLINEABLE genericContinuous #-}
 
 -- | Generic function to create proposals for discrete parameters (e.g., 'Int').
@@ -103,7 +103,7 @@
   --
   -- Only required for biased proposals.
   Maybe (Int -> Int) ->
-  ProposalSimple a
+  PFunction a
 genericDiscrete d f mfInv x g = do
   u <- genDiscreteVar d g
   let r = case mfInv of
@@ -112,5 +112,5 @@
           let qXY = Exp $ logProbability d u
               qYX = Exp $ logProbability d (fInv u)
            in qYX / qXY
-  return (x `f` u, r, 1.0)
+  pure (Propose (x `f` u) r 1.0, Nothing)
 {-# INLINEABLE genericDiscrete #-}
diff --git a/src/Mcmc/Proposal/Hamiltonian.hs b/src/Mcmc/Proposal/Hamiltonian.hs
deleted file mode 100644
--- a/src/Mcmc/Proposal/Hamiltonian.hs
+++ /dev/null
@@ -1,601 +0,0 @@
-{-# LANGUAGE TupleSections #-}
-{-# LANGUAGE UndecidableInstances #-}
-
--- |
--- Module      :  Mcmc.Proposal.Hamiltonian
--- Description :  Hamiltonian Monte Carlo proposal
--- Copyright   :  (c) 2021 Dominik Schrempf
--- License     :  GPL-3.0-or-later
---
--- Maintainer  :  dominik.schrempf@gmail.com
--- Stability   :  experimental
--- Portability :  portable
---
--- Creation date: Mon Jul  5 12:59:42 2021.
---
--- The Hamiltonian Monte Carlo (HMC) proposal.
---
--- For references, see:
---
--- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using
---   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),
---   Handbook of Markov Chain Monte Carlo (2011), CRC press.
---
--- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data
---   analysis (2014), CRC Press.
---
--- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual
---   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).
---
--- NOTE on implementation:
---
--- - The implementation assumes the existence of the 'Gradient'. Like so, the
---   user can use automatic or manual differentiation, depending on the problem
---   at hand.
---
--- - The Hamiltonian proposal acts on a vector of storable 'Values'. Functions
---   converting the state to and from this vector have to be provided. See
---   'HSettings'.
---
--- - The desired acceptance rate is 0.65, although the dimension of the proposal
---   is high.
---
--- - The speed of this proposal can change drastically when tuned because the
---   leapfrog trajectory length is changed.
---
--- - The Hamiltonian proposal is agnostic of the actual prior and likelihood
---   functions, and so, points with zero posterior probability cannot be
---   detected. This affects models with constrained parameters. See Gelman p.
---   303. This problem can be ameliorated by providing a 'Validate' function so
---   that the proposal can gracefully fail as soon as the state becomes invalid.
-module Mcmc.Proposal.Hamiltonian
-  ( Values,
-    Gradient,
-    Validate,
-    Masses,
-    LeapfrogTrajectoryLength,
-    LeapfrogScalingFactor,
-    HTuneLeapfrog (..),
-    HTuneMasses (..),
-    HTune (..),
-    HSettings (..),
-    hamiltonian,
-  )
-where
-
-import qualified Data.Vector as VB
-import qualified Data.Vector.Storable as VS
-import qualified Data.Vector.Unboxed as VU
-import Mcmc.Proposal
-import qualified Numeric.LinearAlgebra as L
-import Numeric.Log
-import Numeric.MathFunctions.Constants
-import qualified Statistics.Covariance as S
-import qualified Statistics.Function as S
-import qualified Statistics.Sample as S
-import System.Random.MWC
-
--- TODO: No-U-turn sampler (NUTS). Ameliorates necessity to determine the
--- leapfrog trajectory length L. (I think this is a necessary extension.)
-
--- TODO: Riemannian adaptation: State-dependent mass matrix. (Seems a little bit
--- of an overkill.)
-
--- | The Hamiltonian proposal acts on a vector of floating point values.
-type Values = L.Vector Double
-
--- | Gradient of the log posterior function.
---
--- The gradient has to be provided for the complete state. The reason is that
--- the gradient may change if parameters untouched by the Hamiltonian proposal
--- are altered by other proposals.
-type Gradient a = a -> a
-
--- | Function validating the state.
---
--- Useful if parameters are constrained.
---
--- Also the validity of the state may depend on parameters untouched by the
--- Hamiltonian proposal.
-type Validate a = a -> Bool
-
--- | Parameter mass matrix.
---
--- The masses roughly describe how reluctant the particles move through the
--- state space. If a parameter has higher mass, the momentum in this direction
--- will be changed less by the provided gradient, than when the same parameter
--- has lower mass. Off-diagonal entries describe the covariance structure. If
--- two parameters are negatively correlated, their generated initial momenta are
--- likely to have opposite signs.
---
--- The proposal is more efficient if masses are assigned according to the
--- inverse (co)-variance structure of the posterior function. That is,
--- parameters changing on larger scales should have lower masses than parameters
--- changing on lower scales. In particular, the optimal entries of the diagonal
--- of the mass matrix are the inverted variances of the parameters distributed
--- according to the posterior function.
---
--- Of course, the scales of the parameters of the posterior function are usually
--- unknown. Often, it is sufficient to
---
--- - set the diagonal entries of the mass matrix to identical values roughly
---   scaled with the inverted estimated average variance of the posterior
---   function; or even to
---
--- - set all diagonal entries of the mass matrix to 1.0, and all other entries
---   to 0.0, and trust the tuning algorithm (see 'HTune') to find the correct
---   values.
-type Masses = L.Herm Double
-
--- | Mean leapfrog trajectory length \(L\).
---
--- Number of leapfrog steps per proposal.
---
--- To avoid problems with ergodicity, the actual number of leapfrog steps is
--- sampled per proposal from a discrete uniform distribution over the interval
--- \([\text{floor}(0.8L),\text{ceiling}(1.2L)]\).
---
--- For a discussion of ergodicity and reasons why randomization is important,
--- see [1] p. 15; also mentioned in [2] p. 304.
---
--- Usually set to 10, but larger values may be desirable.
---
--- NOTE: To avoid errors, the left bound has an additional hard minimum of 1,
--- and the right bound is required to be larger equal than the left bound.
---
--- NOTE: Call 'error' if value is less than 1.
-type LeapfrogTrajectoryLength = Int
-
--- | Mean of leapfrog scaling factor \(\epsilon\).
---
--- Determines the size of each leapfrog step.
---
--- To avoid problems with ergodicity, the actual leapfrog scaling factor is
--- sampled per proposal from a continuous uniform distribution over the interval
--- \((0.8\epsilon,1.2\epsilon]\).
---
--- For a discussion of ergodicity and reasons why randomization is important,
--- see [1] p. 15; also mentioned in [2] p. 304.
---
--- Usually set such that \( L \epsilon = 1.0 \), but smaller values may be
--- required if acceptance rates are low.
---
--- NOTE: Call 'error' if value is zero or negative.
-type LeapfrogScalingFactor = Double
-
--- Internal. Values; target state containing parameters.
-type Positions = Values
-
--- Internal. Momenta of the parameters.
-type Momenta = L.Vector Double
-
--- | Tune leapfrog parameters?
-data HTuneLeapfrog
-  = HNoTuneLeapfrog
-  | -- | We expect that the larger the leapfrog scaling factor the lower the
-    -- acceptance ratio. Consequently, if the acceptance rate is too low, the
-    -- leapfrog scaling factor is decreased and vice versa. Further, the leapfrog
-    -- trajectory length is scaled such that the product of the leapfrog scaling
-    -- factor and leapfrog trajectory length stays roughly constant.
-    HTuneLeapfrog
-  deriving (Eq, Show)
-
--- | Tune masses?
---
--- The masses are tuned according to the (co)variances of the parameters
--- obtained from the posterior distribution over the last auto tuning interval.
-data HTuneMasses
-  = HNoTuneMasses
-  | -- | Diagonal only: The variances of the parameters are calculated and the
-    -- masses are amended using the old masses and the inverted variances. If, for
-    -- a specific coordinate, the sample size is 60 or lower, or if the calculated
-    -- variance is out of predefined bounds [1e-6, 1e6], the mass of the affected
-    -- position is not changed.
-    HTuneDiagonalMassesOnly
-  | -- | All masses: The covariance matrix of the parameters is estimated and the
-    -- inverted matrix (sometimes called precision matrix) is used as mass matrix.
-    -- This procedure is error prone, but models with high correlations between
-    -- parameters it is necessary to tune off-diagonal entries. The full mass
-    -- matrix is only tuned if more than 200 samples are available. For these
-    -- reasons, when tuning all masses it is recommended to use tuning settings
-    -- such as
-    --
-    -- @
-    -- BurnInWithCustomAutoTuning ([10, 20 .. 200] ++ replicate 5 500)
-    -- @
-    HTuneAllMasses
-  deriving (Eq, Show)
-
--- | Tuning settings.
-data HTune = HTune HTuneLeapfrog HTuneMasses
-  deriving (Eq, Show)
-
--- | Specifications of the Hamilton Monte Carlo proposal.
-data HSettings a = HSettings
-  { -- | The sample state is used for error checks and to calculate the dimension
-    -- of the proposal.
-    hSample :: a,
-    -- | Extract values to be manipulated by the Hamiltonian proposal from the
-    -- state.
-    hToVector :: a -> Values,
-    -- | Put those values back into the state.
-    hFromVectorWith :: a -> Values -> a,
-    hGradient :: Gradient a,
-    hMaybeValidate :: Maybe (Validate a),
-    hMasses :: Masses,
-    hLeapfrogTrajectoryLength :: LeapfrogTrajectoryLength,
-    hLeapfrogScalingFactor :: LeapfrogScalingFactor,
-    hTune :: HTune
-  }
-
-checkHSettings :: Eq a => HSettings a -> Maybe String
-checkHSettings (HSettings x toVec fromVec _ _ masses l eps _)
-  | any (<= 0) diagonalMasses = Just "checkHSettings: Some diagonal entries of the mass matrix are zero or negative."
-  | nrows /= ncols = Just "checkHSettings: Mass matrix is not square."
-  | fromVec x xVec /= x = Just "checkHSettings: 'fromVectorWith x (toVector x) /= x' for sample state."
-  | L.size xVec /= nrows = Just "checkHSettings: Mass matrix has different size than 'toVector x', where x is sample state."
-  | l < 1 = Just "checkHSettings: Leapfrog trajectory length is zero or negative."
-  | eps <= 0 = Just "checkHSettings: Leapfrog scaling factor is zero or negative."
-  | otherwise = Nothing
-  where
-    ms = L.unSym masses
-    diagonalMasses = L.toList $ L.takeDiag ms
-    nrows = L.rows ms
-    ncols = L.cols ms
-    xVec = toVec x
-
--- Internal. Mean vector containing zeroes.
-type HMu = L.Vector Double
-
--- Internal. Symmetric, inverted mass matrix.
-type HMassesInv = L.Herm Double
-
--- Internal. Symmetric, inverted mass matrix scaled with the leapfrog step size
--- epsilon.
-type HMassesInvEps = L.Herm Double
-
--- Internal. Logarithm of the determinant of the mass matrix.
-type HLogDetMasses = Double
-
--- Internal data type containing memoized values.
-data HData = HData
-  { _hMu :: HMu,
-    _hMassesInv :: HMassesInv,
-    _hMassesInvEps :: HMassesInvEps,
-    _hLogDetMasses :: HLogDetMasses
-  }
-
--- Call 'error' if the determinant of the covariance matrix is negative.
-getHData :: HSettings a -> HData
-getHData s =
-  -- The multivariate normal distribution requires a positive definite matrix
-  -- with positive determinant.
-  if sign == 1.0
-    then HData mu massesInvH massesInvEpsH logDetMasses
-    else
-      let msg =
-            "hamiltonianSimple: Determinant of covariance matrix is negative."
-              <> " The logarithm of the absolute value of the determinant is: "
-              <> show logDetMasses
-              <> "."
-       in error msg
-  where
-    ms = hMasses s
-    nrows = L.rows $ L.unSym ms
-    mu = L.fromList $ replicate nrows 0.0
-    (massesInv, (logDetMasses, sign)) = L.invlndet $ L.unSym ms
-    -- In theory we can trust that the matrix is symmetric here, because the
-    -- inverse of a symmetric matrix is symmetric. However, one may want to
-    -- implement a check anyways.
-    massesInvH = L.trustSym massesInv
-    eps = hLeapfrogScalingFactor s
-    massesInvEpsH = L.scale eps massesInvH
-
-generateMomenta ::
-  -- Provided so that it does not have to be recreated.
-  HMu ->
-  Masses ->
-  GenIO ->
-  IO Momenta
-generateMomenta mu masses gen = do
-  seed <- uniformM gen :: IO Int
-  let momenta = L.gaussianSample seed 1 mu masses
-  return $ L.flatten momenta
-
--- Prior distribution of momenta.
---
--- Log of density of multivariate normal distribution with given parameters.
--- https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function.
-logDensityMultivariateNormal ::
-  -- Mean vector.
-  L.Vector Double ->
-  -- Inverted covariance matrix.
-  L.Herm Double ->
-  -- Logarithm of the determinant of the covariance matrix.
-  Double ->
-  -- Value vector.
-  L.Vector Double ->
-  Log Double
-logDensityMultivariateNormal mu sigmaInvH logDetSigma xs =
-  Exp $ c + (-0.5) * (logDetSigma + ((dxs L.<# sigmaInv) L.<.> dxs))
-  where
-    dxs = xs - mu
-    k = fromIntegral $ L.size mu
-    c = negate $ m_ln_sqrt_2_pi * k
-    sigmaInv = L.unSym sigmaInvH
-
-leapfrog ::
-  Gradient Positions ->
-  Maybe (Validate Positions) ->
-  HMassesInvEps ->
-  LeapfrogTrajectoryLength ->
-  LeapfrogScalingFactor ->
-  Positions ->
-  Momenta ->
-  -- Maybe (Positions', Momenta'); fail if state is not valid.
-  Maybe (Positions, Momenta)
-leapfrog grad mVal hMassesInvEps l eps theta phi = do
-  let -- The first half step of the momenta.
-      phiHalf = leapfrogStepMomenta (0.5 * eps) grad theta phi
-  -- L-1 full steps. This gives the positions theta_{L-1}, and the momenta
-  -- phi_{L-1/2}.
-  (thetaLM1, phiLM1Half) <- go (l - 1) (Just (theta, phiHalf))
-  -- The last full step of the positions.
-  thetaL <- valF $ leapfrogStepPositions hMassesInvEps thetaLM1 phiLM1Half
-  let -- The last half step of the momenta.
-      phiL = leapfrogStepMomenta (0.5 * eps) grad thetaL phiLM1Half
-  return (thetaL, phiL)
-  where
-    valF x = case mVal of
-      Nothing -> Just x
-      Just f -> if f x then Just x else Nothing
-    go _ Nothing = Nothing
-    go 0 (Just (t, p)) = Just (t, p)
-    go n (Just (t, p)) =
-      let t' = leapfrogStepPositions hMassesInvEps t p
-          p' = leapfrogStepMomenta eps grad t' p
-          r = (,p') <$> valF t'
-       in go (n - 1) r
-
-leapfrogStepMomenta ::
-  -- Step size (usually a multiple of the leapfrog scaling factor).
-  Double ->
-  Gradient Positions ->
-  -- Current positions.
-  Positions ->
-  -- Current momenta.
-  Momenta ->
-  -- New momenta.
-  Momenta
-leapfrogStepMomenta eps grad theta phi = phi + L.scale eps (grad theta)
-
-leapfrogStepPositions ::
-  HMassesInvEps ->
-  -- Current positions.
-  Positions ->
-  -- Current momenta.
-  Momenta ->
-  -- New positions.
-  Positions
-leapfrogStepPositions hMassesInvEps theta phi = theta + (L.unSym hMassesInvEps L.#> phi)
-
-massesToTuningParameters :: Masses -> AuxiliaryTuningParameters
-massesToTuningParameters = VB.convert . L.flatten . L.unSym
-
-tuningParametersToMasses ::
-  -- Dimension of the mass matrix.
-  Int ->
-  AuxiliaryTuningParameters ->
-  Masses
-tuningParametersToMasses d = L.trustSym . L.reshape d . VB.convert
-
-hTuningParametersToSettings ::
-  HSettings a ->
-  TuningParameter ->
-  AuxiliaryTuningParameters ->
-  Either String (HSettings a)
-hTuningParametersToSettings s t ts
-  | nTsNotOK =
-      Left "hTuningParametersToSettings: Auxiliary variables do not have correct dimension."
-  | otherwise =
-      Right $
-        s
-          { hMasses = msTuned,
-            hLeapfrogTrajectoryLength = lTuned,
-            hLeapfrogScalingFactor = eTuned
-          }
-  where
-    ms = hMasses s
-    d = L.rows $ L.unSym ms
-    l = hLeapfrogTrajectoryLength s
-    e = hLeapfrogScalingFactor s
-    (HTune tlf tms) = hTune s
-    nTsNotOK =
-      let nTs = VU.length ts
-       in case tms of
-            HNoTuneMasses -> nTs /= 0
-            _ -> nTs /= d * d
-    msTuned = case tms of
-      HNoTuneMasses -> ms
-      _ -> tuningParametersToMasses d ts
-    -- The larger epsilon, the larger the proposal step size and the lower the
-    -- expected acceptance ratio.
-    --
-    -- Further, we roughly keep \( L * \epsilon = 1.0 \). The equation is not
-    -- correct, because we pull L closer to the original value to keep the
-    -- runtime somewhat acceptable.
-    (lTuned, eTuned) = case tlf of
-      HNoTuneLeapfrog -> (l, e)
-      HTuneLeapfrog -> (ceiling $ fromIntegral l / (t ** 0.9) :: Int, t * e)
-
-hamiltonianSimpleWithTuningParameters ::
-  HSettings a ->
-  TuningParameter ->
-  AuxiliaryTuningParameters ->
-  Either String (ProposalSimple a)
-hamiltonianSimpleWithTuningParameters s t ts =
-  hamiltonianSimple <$> hTuningParametersToSettings s t ts
-
--- The inverted covariance matrix and the log determinant of the covariance
--- matrix are calculated by 'hamiltonianSimple'.
-hamiltonianSimpleWithMemoizedCovariance ::
-  HSettings a ->
-  HData ->
-  ProposalSimple a
-hamiltonianSimpleWithMemoizedCovariance st dt x g = do
-  phi <- generateMomenta mu masses g
-  lRan <- uniformR (lL, lR) g
-  eRan <- uniformR (eL, eR) g
-  case leapfrog gradientVec mValVec massesInvEps lRan eRan theta phi of
-    Nothing -> return (x, 0.0, 1.0)
-    Just (theta', phi') ->
-      let -- Prior of momenta.
-          prPhi = logDensityMultivariateNormal mu massesInv logDetMasses phi
-          -- NOTE: Neal page 12: In order for the proposal to be in detailed
-          -- balance, the momenta have to be negated before proposing the new
-          -- value. This is not required here since the prior involves a
-          -- multivariate normal distribution with means 0.
-          prPhi' = logDensityMultivariateNormal mu massesInv logDetMasses phi'
-          kernelR = prPhi' / prPhi
-       in return (fromVec x theta', kernelR, 1.0)
-  where
-    (HSettings _ toVec fromVec gradient mVal masses l e _) = st
-    theta = toVec x
-    lL = maximum [1 :: Int, floor $ (0.8 :: Double) * fromIntegral l]
-    lR = maximum [lL, ceiling $ (1.2 :: Double) * fromIntegral l]
-    eL = 0.8 * e
-    eR = 1.2 * e
-    (HData mu massesInv massesInvEps logDetMasses) = dt
-    -- Vectorize the gradient and validation functions.
-    gradientVec = toVec . gradient . fromVec x
-    mValVec = mVal >>= (\f -> return $ f . fromVec x)
-
-hamiltonianSimple ::
-  HSettings a ->
-  ProposalSimple a
-hamiltonianSimple s = hamiltonianSimpleWithMemoizedCovariance s hd
-  where
-    hd = getHData s
-
--- If changed, also change help text of 'HTuneMasses'.
-massMin :: Double
-massMin = 1e-6
-
--- If changed, also change help text of 'HTuneMasses'.
-massMax :: Double
-massMax = 1e6
-
--- Minimal number of unique samples required for tuning the diagonal entries of
--- the mass matrix.
---
--- If changed, also change help text of 'HTuneMasses'.
-samplesMinDiagonal :: Int
-samplesMinDiagonal = 61
-
--- Minimal number of samples required for tuning all entries of the mass matrix.
---
--- If changed, also change help text of 'HTuneMasses'.
-samplesMinAll :: Int
-samplesMinAll = 201
-
-getSampleSize :: VS.Vector Double -> Int
-getSampleSize = VS.length . VS.uniq . S.gsort
-
--- Diagonal elements are variances which are strictly positive.
-getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double
-getNewMassDiagonalWithRescue sampleSize massOld massEstimate
-  | sampleSize < samplesMinDiagonal = massOld
-  -- NaN and negative masses could be errors.
-  | isNaN massEstimate = massOld
-  | massEstimate <= 0 = massOld
-  | massMin > massNew = massMin
-  | massNew > massMax = massMax
-  | otherwise = massNew
-  where
-    massNewSqrt = recip 3 * (sqrt massOld + 2 * sqrt massEstimate)
-    massNew = massNewSqrt ** 2
-
--- NOTE: Here, we lose time because we convert the states to vectors again,
--- something that has already been done. But then, auto tuning is not a runtime
--- determining factor.
-tuneDiagonalMassesOnly ::
-  Int ->
-  (a -> Positions) ->
-  AuxiliaryTuningFunction a
-tuneDiagonalMassesOnly dim toVec xs ts
-  -- If not enough data is available, do not tune.
-  | VB.length xs < samplesMinDiagonal = ts
-  | otherwise =
-      -- Replace the diagonal.
-      massesToTuningParameters $
-        L.trustSym $ massesOld - L.diag massesDiagonalOld + L.diag massesDiagonalNew
-  where
-    -- xs: Each vector entry contains all parameter values of one iteration.
-    -- xs': Each row contains all parameter values of one iteration.
-    xs' = L.fromRows $ VB.toList $ VB.map toVec xs
-    sampleSizes = VS.fromList $ map getSampleSize $ L.toColumns xs'
-    massesOld = L.unSym $ tuningParametersToMasses dim ts
-    massesDiagonalOld = L.takeDiag massesOld
-    massesDiagonalEstimate = VS.fromList $ map (recip . S.variance) $ L.toColumns xs'
-    massesDiagonalNew =
-      VS.zipWith3
-        getNewMassDiagonalWithRescue
-        sampleSizes
-        massesDiagonalOld
-        massesDiagonalEstimate
-
--- NOTE: Here, we lose time because we convert the states to vectors again,
--- something that has already been done. But then, auto tuning is not a runtime
--- determining factor.
-tuneAllMasses ::
-  Int ->
-  (a -> Positions) ->
-  AuxiliaryTuningFunction a
-tuneAllMasses dim toVec xs ts
-  -- If not enough data is available, do not tune.
-  | VB.length xs < samplesMinDiagonal = ts
-  -- If not enough data is available, only the diagonal masses are tuned.
-  | VB.length xs < samplesMinAll = fallbackDiagonal
-  | L.rank xs' /= dim = fallbackDiagonal
-  | otherwise = massesToTuningParameters $ L.trustSym massesNew
-  where
-    fallbackDiagonal = tuneDiagonalMassesOnly dim toVec xs ts
-    -- xs: Each vector entry contains all parameter values of one iteration.
-    -- xs': Each row contains all parameter values of one iteration.
-    xs' = L.fromRows $ VB.toList $ VB.map toVec xs
-    (_, ss, xsNormalized) = S.scale xs'
-    -- sigmaNormalized = L.unSym $ either error id $ S.oracleApproximatingShrinkage xsNormalized
-    sigmaNormalized = L.unSym $ either error fst $ S.graphicalLasso 0.5 xsNormalized
-    sigma = S.rescaleSWith ss sigmaNormalized
-    massesNew = L.inv sigma
-
--- | Hamiltonian Monte Carlo proposal.
-hamiltonian ::
-  Eq a =>
-  HSettings a ->
-  PName ->
-  PWeight ->
-  Proposal a
-hamiltonian s n w = case checkHSettings s of
-  Just err -> error err
-  Nothing ->
-    let desc = PDescription "Hamiltonian Monte Carlo (HMC)"
-        toVec = hToVector s
-        dim = (L.size $ toVec $ hSample s)
-        pDim = PSpecial dim 0.65
-        ts = massesToTuningParameters (hMasses s)
-        ps = hamiltonianSimple s
-        hamiltonianWith = Proposal n desc PSlow pDim w ps
-        tSet@(HTune tlf tms) = hTune s
-        tFun = case tlf of
-          HNoTuneLeapfrog -> noTuningFunction
-          HTuneLeapfrog -> defaultTuningFunctionWith pDim
-        tFunAux = case tms of
-          HNoTuneMasses -> noAuxiliaryTuningFunction
-          HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly dim toVec
-          HTuneAllMasses -> tuneAllMasses dim toVec
-     in case tSet of
-          (HTune HNoTuneLeapfrog HNoTuneMasses) -> hamiltonianWith Nothing
-          _ ->
-            let tuner = Tuner 1.0 tFun ts tFunAux (hamiltonianSimpleWithTuningParameters s)
-             in hamiltonianWith $ Just tuner
diff --git a/src/Mcmc/Proposal/Hamiltonian/Common.hs b/src/Mcmc/Proposal/Hamiltonian/Common.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Proposal/Hamiltonian/Common.hs
@@ -0,0 +1,215 @@
+{-# LANGUAGE RankNTypes #-}
+
+-- |
+-- Module      :  Mcmc.Proposal.Hamiltonian.Common
+-- Description :  Code shared by proposals based on Hamiltonian dynamics
+-- Copyright   :  2022 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- Creation date: Fri Jun  3 10:27:12 2022.
+--
+-- See "Mcmc.Proposal.Hamiltonian.Hamiltonian".
+--
+-- References:
+--
+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using
+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),
+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.
+--
+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data
+--   analysis (2014), CRC Press.
+--
+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual
+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).
+--
+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:
+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of
+--   Machine Learning Research.
+module Mcmc.Proposal.Hamiltonian.Common
+  ( -- * Basic types
+    Positions,
+    Momenta,
+    Masses,
+    LeapfrogTrajectoryLength,
+    LeapfrogScalingFactor,
+    LeapfrogSimulationLength,
+
+    -- * Tuning
+    HTuneLeapfrog (..),
+    HTuneMasses (..),
+    HTuningConf (..),
+
+    -- * Structure of state
+    HStructure (..),
+
+    -- * Target function
+    HTarget (..),
+  )
+where
+
+import Data.Typeable
+import qualified Data.Vector.Storable as VS
+import Mcmc.Jacobian
+import Mcmc.Likelihood
+import Mcmc.Prior
+import qualified Numeric.LinearAlgebra as L
+
+-- NOTE: Implementing the Riemannian adaptation (state-dependent mass matrix).
+-- seems a little bit of an overkill.
+
+-- | The Hamiltonian proposal acts on a vector of floating point values referred
+-- to as positions.
+--
+-- The positions can represent the complete state or a subset of the state of
+-- the Markov chain; see 'HStructure'.
+--
+-- The length of the position vector determines the size of the squared mass
+-- matrix; see 'Masses'.
+type Positions = VS.Vector Double
+
+-- | Momenta of the 'Positions'.
+type Momenta = VS.Vector Double
+
+-- | Mass matrix.
+--
+-- The masses roughly describe how reluctant the particles move through the
+-- state space. If a parameter has higher mass the velocity in this direction
+-- will be changed less by the provided gradient than when the same parameter
+-- has lower mass. Off-diagonal entries describe the covariance structure. If
+-- two parameters are negatively correlated, their generated initial momenta are
+-- likely to have opposite signs.
+--
+-- The matrix is square with the number of rows and columns being equal to the
+-- length of 'Positions'.
+--
+-- The proposal is more efficient if masses are assigned according to the
+-- inverse (co)-variance structure of the posterior function. That is,
+-- parameters changing on larger scales should have lower masses than parameters
+-- changing on lower scales. In particular, the optimal entries of the diagonal
+-- of the mass matrix are the inverted variances of the parameters distributed
+-- according to the posterior function.
+--
+-- Of course, the scales of the parameters of the posterior function are usually
+-- unknown. Often, it is sufficient to
+--
+-- - set the diagonal entries of the mass matrix to identical values roughly
+--   scaled with the inverted estimated average variance of the posterior
+--   function; or even to
+--
+-- - set all diagonal entries of the mass matrix to 1.0, and all other entries
+--   to 0.0, and trust the tuning algorithm (see 'HTuningConf') to find the
+--   correct values.
+type Masses = L.Herm Double
+
+-- | Mean leapfrog trajectory length \(L\).
+--
+-- Number of leapfrog steps per proposal.
+--
+-- To avoid problems with ergodicity, the actual number of leapfrog steps is
+-- sampled per proposal from a discrete uniform distribution over the interval
+-- \([\text{floor}(0.9L),\text{ceiling}(1.1L)]\).
+--
+-- For a discussion of ergodicity and reasons why randomization is important,
+-- see [1] p. 15; also mentioned in [2] p. 304.
+--
+-- To avoid errors, the left bound of the interval has an additional hard
+-- minimum of 1, and the right bound is required to be larger equal than the
+-- left bound.
+type LeapfrogTrajectoryLength = Int
+
+-- | Mean of leapfrog scaling factor \(\epsilon\).
+--
+-- Determines the size of each leapfrog step.
+--
+-- To avoid problems with ergodicity, the actual leapfrog scaling factor is
+-- sampled per proposal from a continuous uniform distribution over the interval
+-- \((0.9\epsilon,1.1\epsilon]\).
+--
+-- For a discussion of ergodicity and reasons why randomization is important,
+-- see [1] p. 15; also mentioned in [2] p. 304.
+--
+-- Call 'error' if value is zero or negative.
+type LeapfrogScalingFactor = Double
+
+-- | Product of 'LeapfrogTrajectoryLength' and 'LeapfrogScalingFactor'.
+--
+-- A good value is hard to find and varies between applications. For example,
+-- see Figure 6 in [4].
+--
+-- Call 'error' if value is zero or negative.
+type LeapfrogSimulationLength = Double
+
+-- | Tune leapfrog parameters?
+data HTuneLeapfrog
+  = HNoTuneLeapfrog
+  | -- | We expect that the larger the leapfrog scaling factor the lower the
+    -- acceptance ratio. Consequently, if the acceptance rate is too low, the
+    -- leapfrog scaling factor is decreased and vice versa. Further, the leapfrog
+    -- trajectory length is scaled such that the product of the leapfrog scaling
+    -- factor and leapfrog trajectory length stays constant.
+    HTuneLeapfrog
+  deriving (Eq, Show)
+
+-- | Tune masses?
+--
+-- The masses are tuned according to the (co)variances of the parameters
+-- obtained from the posterior distribution over the last auto tuning interval.
+data HTuneMasses
+  = HNoTuneMasses
+  | -- | Diagonal only: The variances of the parameters are calculated and the
+    -- masses are amended using the old masses and the inverted variances. If,
+    -- for a specific coordinate, the sample size is 60 or lower, or if the
+    -- calculated variance is out of predefined bounds [1e-8, 1e8], the mass of
+    -- the affected position is not changed.
+    HTuneDiagonalMassesOnly
+  | -- | All masses: The covariance matrix of the parameters is estimated and
+    -- the inverted matrix (sometimes called precision matrix) is used as mass
+    -- matrix. This procedure is error prone, but models with high correlations
+    -- between parameters strongly profit from tuning off-diagonal entries. The
+    -- full mass matrix is only tuned if equal or more than "(n_masses +
+    -- max(n_masses, 61)" samples are available. For these reasons, when tuning
+    -- all masses it is recommended to use tuning settings such as
+    --
+    -- @
+    -- BurnInWithCustomAutoTuning [10, 20 .. 200] [200, 220 .. 500]
+    -- @
+    HTuneAllMasses
+  deriving (Eq, Show)
+
+-- | Tuning configuration of the Hamilton Monte Carlo proposal.
+data HTuningConf = HTuningConf HTuneLeapfrog HTuneMasses
+  deriving (Eq, Show)
+
+-- | The Hamilton Monte Carlo proposal requires information about the structure
+-- of the state, which is denoted as @s@.
+--
+-- Please also refer to the top level module documentation of
+-- "Mcmc.Proposal.Hamiltonian.Hamiltonian".
+data HStructure s = HStructure
+  { -- | The sample state is used for error checks and to calculate the dimension
+    -- of the proposal.
+    hSample :: s Double,
+    -- | Extract a a subset of values to be manipulated by the Hamiltonian
+    -- proposal from the complete state.
+    hToVector :: s Double -> Positions,
+    -- | Put those values back into the complete state.
+    hFromVectorWith :: s Double -> Positions -> s Double
+  }
+
+-- | The target is composed of the prior, likelihood, and jacobian functions.
+--
+-- The structure of the state is denoted as @s@.
+--
+-- Please also refer to the top level module documentation.
+data HTarget s = HTarget
+  { -- | Function computing the log prior.
+    hPrior :: forall a. (RealFloat a, Typeable a) => Maybe (PriorFunctionG (s a) a),
+    -- | Function computing the log likelihood.
+    hLikelihood :: forall a. (RealFloat a, Typeable a) => LikelihoodFunctionG (s a) a,
+    -- | Function computing the log of the Jacobian.
+    hJacobian :: forall a. (RealFloat a, Typeable a) => Maybe (JacobianFunctionG (s a) a)
+  }
diff --git a/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs b/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs
@@ -0,0 +1,210 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- |
+-- Module      :  Mcmc.Proposal.Hamiltonian.Hamiltonian
+-- Description :  Hamiltonian Monte Carlo proposal
+-- Copyright   :  2021 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- Creation date: Mon Jul  5 12:59:42 2021.
+--
+-- The Hamiltonian Monte Carlo (HMC) proposal.
+--
+-- The HMC proposal acts on 'Positions', a vector of floating point values. The
+-- manipulated values can represent the complete state, or a subset of the
+-- complete state. Functions converting the state to and from this vector have
+-- to be provided; see 'HStructure'.
+--
+-- Even though the proposal may only act on a subset of the complete state, the
+-- prior, likelihood, and Jacobian functions of the complete state have to be
+-- provided; see 'HTarget'. This is because parameters not manipulated by the
+-- HMC proposal still influence the prior, likelihood and Jacobian functions.
+--
+-- The points given above have implications on how the HMC proposal is handled:
+-- Do not use 'liftProposalWith', 'liftProposal', or '(@~)' with the HMC
+-- proposal; instead provide proper conversion functions with 'HStructure'.
+--
+-- The gradient of the log target function is calculated using automatic
+-- differentiation; see the excellent
+-- [ad](https://hackage.haskell.org/package/ad) package.
+--
+-- The desired acceptance rate is 0.65, although the dimension of the proposal
+-- is high.
+--
+-- The speed of this proposal changes drastically with the leapfrog trajectory
+-- length and the leapfrog scaling factor. Hence, the speed will change during
+-- burn in.
+--
+-- References:
+--
+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using
+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),
+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.
+--
+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data
+--   analysis (2014), CRC Press.
+--
+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual
+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).
+--
+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:
+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of
+--   Machine Learning Research.
+module Mcmc.Proposal.Hamiltonian.Hamiltonian
+  ( -- * Hamiltonian Monte Carlo proposal
+    HParams (..),
+    defaultHParams,
+    hamiltonian,
+  )
+where
+
+import Data.Bifunctor
+import Mcmc.Acceptance
+import Mcmc.Algorithm.MHG
+import Mcmc.Proposal
+import Mcmc.Proposal.Hamiltonian.Common
+import Mcmc.Proposal.Hamiltonian.Internal
+import Mcmc.Proposal.Hamiltonian.Masses
+import Numeric.AD.Double
+import qualified Numeric.LinearAlgebra as L
+import Numeric.Log
+import System.Random.Stateful
+
+-- | Parameters of the Hamilton Monte Carlo proposal.
+--
+-- If a parameter is 'Nothing', a default value is used (see 'defaultHParams').
+data HParams = HParams
+  { hLeapfrogScalingFactor :: Maybe LeapfrogScalingFactor,
+    hLeapfrogSimulationLength :: Maybe LeapfrogSimulationLength,
+    hMasses :: Maybe Masses
+  }
+  deriving (Show)
+
+-- | Default parameters.
+--
+-- - Estimate a reasonable leapfrog scaling factor using Algorithm 4 [4]. If all
+--   fails, use 0.1.
+--
+-- - Leapfrog simulation length is set to 0.5.
+--
+-- - The mass matrix is set to the identity matrix.
+defaultHParams :: HParams
+defaultHParams = HParams Nothing Nothing Nothing
+
+hamiltonianPFunctionWithTuningParameters ::
+  Traversable s =>
+  Dimension ->
+  HStructure s ->
+  (s Double -> Target) ->
+  TuningParameter ->
+  AuxiliaryTuningParameters ->
+  Either String (PFunction (s Double))
+hamiltonianPFunctionWithTuningParameters d hstruct targetWith _ ts = do
+  hParamsI <- fromAuxiliaryTuningParameters d ts
+  pure $ hamiltonianPFunction hParamsI hstruct targetWith
+
+-- TODO @Dominik (high, issue): Acceptance counts. How to combine with values
+-- reported here and from the NUTS sampler.
+
+-- TODO @Dominik (high, feature): The expected acceptance counts should not be
+-- calculated after burn in. Rather, the actual acceptance counts should be
+-- reported. For this to work, the proposal needs to know if it is in "burn in
+-- phase" or not.
+
+-- The inverted covariance matrix and the log determinant of the covariance
+-- matrix are calculated by 'hamiltonianPFunction'.
+hamiltonianPFunction ::
+  HParamsI ->
+  HStructure s ->
+  (s Double -> Target) ->
+  PFunction (s Double)
+hamiltonianPFunction hparamsi hstruct targetWith x g = do
+  p <- generateMomenta mus ms g
+  eRan <- uniformRM (eL, eR) g
+  -- NOTE: The NUTS paper does not sample l since l varies naturally because
+  -- of epsilon. I still think it should vary because otherwise, there may be
+  -- dragons due to periodicity.
+  let lM = la / eRan
+      lL = maximum [1 :: Int, floor $ 0.9 * lM]
+      lR = maximum [lL, ceiling $ 1.1 * lM]
+  lRan <- uniformRM (lL, lR) g
+  case leapfrog (targetWith x) msI lRan eRan q p of
+    Nothing -> pure (ForceReject, Just $ AcceptanceCounts 0 100)
+    -- Check if next state is accepted here, because the Jacobian is included in
+    -- the target function. If not: pure (x, 0.0, 1.0).
+    Just (q', p', prQ, prQ') -> do
+      let -- Prior of momenta.
+          prP = exponentialKineticEnergy msI p
+          prP' = exponentialKineticEnergy msI p'
+          r = prQ' * prP' / (prQ * prP)
+      accept <- mhgAccept r g
+      -- NOTE: For example, Neal page 12: In order for the Hamiltonian proposal
+      -- to be in detailed balance, the momenta have to be negated before
+      -- proposing the new value. That is, the negated momenta would guide the
+      -- chain back to the previous state. However, we are only interested in
+      -- the positions, and are not even storing the momenta.
+      let pr = if accept then ForceAccept (fromVec x q') else ForceReject
+          ar = exp $ ln r
+          getCounts s = max 0 $ min 100 $ round $ s * 100
+          ac =
+            if ar >= 0
+              then let cs = getCounts ar in AcceptanceCounts cs (100 - cs)
+              else error $ "hamiltonianPFunction: Acceptance rate negative."
+      pure (pr, Just ac)
+  where
+    (HParamsI e la ms _ _ msI mus) = hparamsi
+    (HStructure _ toVec fromVec) = hstruct
+    q = toVec x
+    eL = 0.9 * e
+    eR = 1.1 * e
+
+-- | Hamiltonian Monte Carlo proposal.
+--
+-- The structure of the state is denoted as @s@.
+--
+-- May call 'error' during initialization.
+hamiltonian ::
+  Traversable s =>
+  HParams ->
+  HTuningConf ->
+  HStructure s ->
+  HTarget s ->
+  PName ->
+  PWeight ->
+  Proposal (s Double)
+hamiltonian hparams htconf hstruct htarget n w =
+  let -- Misc.
+      desc = PDescription "Hamiltonian Monte Carlo (HMC)"
+      (HStructure sample toVec fromVec) = hstruct
+      dim = L.size $ toVec sample
+      -- See bottom of page 1615 in [4].
+      pDim = PSpecial dim 0.65
+      -- Vectorize and derive the target function.
+      (HTarget mPrF lhF mJcF) = htarget
+      tF y = case (mPrF, mJcF) of
+        (Nothing, Nothing) -> lhF y
+        (Just prF, Nothing) -> prF y * lhF y
+        (Nothing, Just jcF) -> lhF y * jcF y
+        (Just prF, Just jcF) -> prF y * lhF y * jcF y
+      tFnG = grad' (ln . tF)
+      targetWith x = bimap Exp toVec . tFnG . fromVec x
+      (HParams mEps mLa mMs) = hparams
+      hParamsI =
+        either error id $
+          hParamsIWith (targetWith sample) (toVec sample) mEps mLa mMs
+      ps = hamiltonianPFunction hParamsI hstruct targetWith
+      hamiltonianWith = Proposal n desc PSlow pDim w ps
+      -- Tuning.
+      ts = toAuxiliaryTuningParameters hParamsI
+      tuner = do
+        tfun <- hTuningFunctionWith dim toVec htconf
+        let pfun = hamiltonianPFunctionWithTuningParameters dim hstruct targetWith
+        pure $ Tuner 1.0 ts True tfun pfun
+   in case checkHStructureWith (hpsMasses hParamsI) hstruct of
+        Just err -> error err
+        Nothing -> hamiltonianWith tuner
diff --git a/src/Mcmc/Proposal/Hamiltonian/Internal.hs b/src/Mcmc/Proposal/Hamiltonian/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Proposal/Hamiltonian/Internal.hs
@@ -0,0 +1,413 @@
+-- Creation date: Thu Jun  9 15:12:39 2022.
+--
+-- See "Mcmc.Proposal.Hamiltonian.Hamiltonian".
+--
+-- References:
+--
+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using
+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),
+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.
+--
+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data
+--   analysis (2014), CRC Press.
+--
+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual
+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).
+--
+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:
+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of
+--   Machine Learning Research.
+
+-- |
+-- Module      :  Mcmc.Proposal.Hamiltonian.Internal
+-- Description :  Internal definitions related to Hamiltonian dynamics
+-- Copyright   :  2022 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+module Mcmc.Proposal.Hamiltonian.Internal
+  ( -- * Parameters
+    HParamsI (..),
+    hParamsIWith,
+
+    -- * Tuning
+    toAuxiliaryTuningParameters,
+    fromAuxiliaryTuningParameters,
+    findReasonableEpsilon,
+    hTuningFunctionWith,
+
+    -- * Structure of state
+    checkHStructureWith,
+
+    -- * Hamiltonian dynamics
+    generateMomenta,
+    exponentialKineticEnergy,
+
+    -- * Leapfrog integrator
+    Target,
+    leapfrog,
+  )
+where
+
+import Control.Monad
+import Control.Monad.ST
+import Data.Foldable
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Unboxed as VU
+import Mcmc.Proposal
+import Mcmc.Proposal.Hamiltonian.Common
+import Mcmc.Proposal.Hamiltonian.Masses
+import qualified Numeric.LinearAlgebra as L
+import Numeric.Log
+import System.Random.Stateful
+
+-- Variable tuning parameters.
+--
+-- See Algorithm 5 or 6 in [4].
+data TParamsVar = TParamsVar
+  { -- \bar{eps} of Algorithm 5 or 6.
+    tpvLeapfrogScalingFactorMean :: LeapfrogScalingFactor,
+    -- H_i of Algorithm 5 or 6.
+    tpvHStatistics :: Double,
+    -- m of Algorithm 5 or 6.
+    tpvCurrentTuningStep :: Double
+  }
+  deriving (Show)
+
+tParamsVar :: TParamsVar
+tParamsVar = TParamsVar 1.0 0.0 1.0
+
+-- Fixed tuning parameters.
+--
+-- See Algorithm 5 and 6 in [4].
+data TParamsFixed = TParamsFixed
+  { tpfEps0 :: Double,
+    tpfMu :: Double,
+    tpfGa :: Double,
+    tpfT0 :: Double,
+    tpfKa :: Double
+  }
+  deriving (Show)
+
+-- The default tuning parameters in [4] are:
+--
+--   mu = log $ 10 * eps
+--   ga = 0.05
+--   t0 = 10
+--   ka = 0.75
+--
+-- However, these default tuning parameters will be off, because the authors
+-- suggesting these values tune the proposal after every single iteration.
+--
+-- The following values are tweaked for our case, where tuning does not happen
+-- after each iteration. Of course, we could tune the leapfrog parameters after
+-- each generation. Even the mass parameters could be tuned each iteration when
+-- the masses are estimated from more past iterations spanning many tuning
+-- intervals.
+--
+-- NOTE: In theory, these we could expose these internal tuning parameters to
+-- the user.
+tParamsFixedWith :: LeapfrogScalingFactor -> TParamsFixed
+tParamsFixedWith eps = TParamsFixed eps mu ga t0 ka
+  where
+    mu = log $ 10 * eps
+    ga = 0.1
+    t0 = 3
+    ka = 0.5
+
+-- All internal parameters.
+data HParamsI = HParamsI
+  { hpsLeapfrogScalingFactor :: LeapfrogScalingFactor,
+    hpsLeapfrogSimulationLength :: LeapfrogSimulationLength,
+    hpsMasses :: Masses,
+    hpsTParamsVar :: TParamsVar,
+    hpsTParamsFixed :: TParamsFixed,
+    hpsMassesI :: MassesI,
+    hpsMu :: Mu
+  }
+  deriving (Show)
+
+-- NOTE: If changed, amend help text of 'defaultHParams', and 'defaultNParams'.
+defaultLeapfrogScalingFactor :: LeapfrogScalingFactor
+defaultLeapfrogScalingFactor = 0.1
+
+-- NOTE: If changed, amend help text of 'defaultHParams'.
+defaultLeapfrogSimulationLength :: LeapfrogSimulationLength
+defaultLeapfrogSimulationLength = 0.5
+
+-- NOTE: If changed, amend help text of 'defaultHParams'.
+defaultMassesWith :: Int -> Masses
+defaultMassesWith d = L.trustSym $ L.ident d
+
+-- Instantiate all internal parameters.
+hParamsIWith ::
+  Target ->
+  Positions ->
+  Maybe LeapfrogScalingFactor ->
+  Maybe LeapfrogSimulationLength ->
+  Maybe Masses ->
+  Either String HParamsI
+hParamsIWith htarget p mEps mLa mMs = do
+  d <- case VS.length p of
+    0 -> eWith "Empty position vector."
+    d -> Right d
+  ms <- case mMs of
+    Nothing -> Right $ defaultMassesWith d
+    Just ms -> do
+      let ms' = cleanMatrix $ L.unSym ms
+          diagonalMs = L.toList $ L.takeDiag ms'
+      when (any (<= 0) diagonalMs) $ eWith "Some diagonal masses are zero or negative."
+      let nrows = L.rows ms'
+          ncols = L.cols ms'
+      when (nrows /= ncols) $ eWith "Mass matrix is not square."
+      Right ms
+  let msI = getMassesI ms
+      mus = getMus ms
+  la <- case mLa of
+    Nothing -> Right defaultLeapfrogSimulationLength
+    Just l
+      | l <= 0 -> eWith "Leapfrog simulation length is zero or negative."
+      | otherwise -> Right l
+  eps <- case mEps of
+    Nothing -> Right $ runST $ do
+      -- NOTE: This is not random. However, I do not want to provide a generator
+      -- when creating the proposal.
+      g <- newSTGenM $ mkStdGen 42
+      findReasonableEpsilon htarget ms p g
+    Just e
+      | e <= 0 -> eWith "Leapfrog scaling factor is zero or negative."
+      | otherwise -> Right e
+  let tParamsFixed = tParamsFixedWith eps
+  pure $ HParamsI eps la ms tParamsVar tParamsFixed msI mus
+  where
+    eWith m = Left $ "hParamsIWith: " <> m
+
+-- Save internal parameters.
+toAuxiliaryTuningParameters :: HParamsI -> AuxiliaryTuningParameters
+toAuxiliaryTuningParameters (HParamsI eps la ms tpv tpf _ _) =
+  -- Put masses to the end. Like so, conversion is easier.
+  VU.fromList $ eps : la : epsMean : h : m : eps0 : mu : ga : t0 : ka : msL
+  where
+    (TParamsVar epsMean h m) = tpv
+    (TParamsFixed eps0 mu ga t0 ka) = tpf
+    msL = VU.toList $ massesToVector ms
+
+-- Load internal parameters.
+fromAuxiliaryTuningParameters :: Dimension -> AuxiliaryTuningParameters -> Either String HParamsI
+fromAuxiliaryTuningParameters d xs
+  | (d * d) + 10 /= len = Left "fromAuxiliaryTuningParameters: Dimension mismatch."
+  | fromIntegral (d * d) /= lenMs = Left "fromAuxiliaryTuningParameters: Masses dimension mismatch."
+  | otherwise = case VU.toList $ VU.take 10 xs of
+      [eps, la, epsMean, h, m, eps0, mu, ga, t0, ka] ->
+        let tpv = TParamsVar epsMean h m
+            tpf = TParamsFixed eps0 mu ga t0 ka
+         in Right $ HParamsI eps la ms tpv tpf msI mus
+      -- To please the exhaustive pattern match checker.
+      _ -> Left "fromAuxiliaryTuningParameters: Impossible dimension mismatch."
+  where
+    len = VU.length xs
+    msV = VU.drop 10 xs
+    lenMs = VU.length msV
+    ms = vectorToMasses d msV
+    msI = getMassesI ms
+    mus = getMus ms
+
+-- See Algorithm 4 in [4].
+findReasonableEpsilon ::
+  StatefulGen g m =>
+  Target ->
+  Masses ->
+  Positions ->
+  g ->
+  m LeapfrogScalingFactor
+findReasonableEpsilon t ms q g = do
+  p <- generateMomenta mu ms g
+  case leapfrog t msI 1 eI q p of
+    Nothing -> pure defaultLeapfrogScalingFactor
+    Just (_, p', prQ, prQ') -> do
+      let expEKin = exponentialKineticEnergy msI p
+          expEKin' = exponentialKineticEnergy msI p'
+          rI :: Double
+          rI = exp $ ln $ prQ' * expEKin' / (prQ * expEKin)
+          a :: Double
+          a = if rI > 0.5 then 1 else (-1)
+          go e r =
+            if r ** a > 2 ** (negate a)
+              then case leapfrog t msI 1 e q p of
+                Nothing -> e
+                Just (_, p'', _, prQ'') ->
+                  let expEKin'' = exponentialKineticEnergy msI p''
+                      r' :: Double
+                      r' = exp $ ln $ prQ'' * expEKin'' / (prQ * expEKin)
+                      e' = (2 ** a) * e
+                   in go e' r'
+              else e
+      pure $ go eI rI
+  where
+    eI = 1.0
+    msI = getMassesI ms
+    mu = getMus ms
+
+hTuningFunctionWith ::
+  Dimension ->
+  -- Conversion from value to vector.
+  (a -> Positions) ->
+  HTuningConf ->
+  Maybe (TuningFunction a)
+hTuningFunctionWith n toVec (HTuningConf lc mc) = case (lc, mc) of
+  (HNoTuneLeapfrog, HNoTuneMasses) -> Nothing
+  (_, _) -> Just $
+    \tt pdim ar mxs (_, ts) ->
+      case mxs of
+        Nothing -> error "hTuningFunctionWith: empty trace"
+        Just xs ->
+          let (HParamsI eps la ms tpv tpf msI mus) =
+                -- NOTE: Use error here, because a dimension mismatch is a serious bug.
+                either error id $ fromAuxiliaryTuningParameters n ts
+              (TParamsVar epsMean h m) = tpv
+              (TParamsFixed eps0 mu ga t0 ka) = tpf
+              (ms', msI') = case mc of
+                HNoTuneMasses -> (ms, msI)
+                HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly toVec xs (ms, msI)
+                HTuneAllMasses -> tuneAllMasses toVec xs (ms, msI)
+              (eps'', epsMean'', h'') = case lc of
+                HNoTuneLeapfrog -> (eps, epsMean, h)
+                HTuneLeapfrog ->
+                  let delta = getOptimalRate pdim
+                      c = recip $ m + t0
+                      h' = (1.0 - c) * h + c * (delta - ar)
+                      logEps' = mu - (sqrt m / ga) * h'
+                      eps' = exp logEps'
+                      mMKa = m ** (negate ka)
+                      epsMean' = exp $ mMKa * logEps' + (1 - mMKa) * log epsMean
+                   in (eps', epsMean', h')
+              eps''' = case tt of
+                NormalTuningStep -> eps''
+                LastTuningStep -> epsMean''
+              tpv' = TParamsVar epsMean'' h'' (m + 1.0)
+           in (eps''' / eps0, toAuxiliaryTuningParameters $ HParamsI eps''' la ms' tpv' tpf msI' mus)
+
+checkHStructureWith :: Foldable s => Masses -> HStructure s -> Maybe String
+checkHStructureWith ms (HStructure x toVec fromVec)
+  | toList (fromVec x xVec) /= toList x = eWith "'fromVectorWith x (toVector x) /= x' for sample state."
+  | L.size xVec /= nrows = eWith "Mass matrix and 'toVector x' have different sizes for sample state."
+  | otherwise = Nothing
+  where
+    eWith m = Just $ "checkHStructureWith: " <> m
+    nrows = L.rows $ L.unSym ms
+    xVec = toVec x
+
+-- Generate momenta for a new iteration.
+generateMomenta ::
+  StatefulGen g m =>
+  Mu ->
+  Masses ->
+  g ->
+  m Momenta
+generateMomenta mu masses gen = do
+  seed <- uniformM gen
+  let momenta = L.gaussianSample seed 1 mu masses
+  return $ L.flatten momenta
+
+-- Compute exponent of kinetic energy.
+--
+-- Use a general matrix which has special representations for diagonal and
+-- sparse matrices, both of which are really useful here.
+exponentialKineticEnergy ::
+  MassesI ->
+  Momenta ->
+  Log Double
+exponentialKineticEnergy msI xs =
+  -- NOTE: Because of numerical errors, the following formulas exhibit different
+  -- traces (although the posterior appears to be the same):
+  -- - This one we cannot use with general matrices:
+  --   Exp $ (-0.5) * ((xs L.#> msI) L.<.> xs)
+  Exp $ (-0.5) * (xs L.<.> (msI L.!#> xs))
+
+-- Function calculating target value and gradient.
+--
+-- The function acts on the subset of the state manipulated by the proposal but
+-- the value and gradient have to be calculated for the complete state. The
+-- reason is that parameters untouched by the Hamiltonian proposal may affect
+-- the result or the gradient.
+--
+-- Make sure that the value is calculated lazily because many times, only the
+-- gradient is required.
+type Target = Positions -> (Log Double, Positions)
+
+-- Leapfrog integrator.
+leapfrog ::
+  Target ->
+  MassesI ->
+  --
+  LeapfrogTrajectoryLength ->
+  LeapfrogScalingFactor ->
+  --
+  Positions ->
+  Momenta ->
+  -- | (New positions, new momenta, old target, new target).
+  --
+  -- Fail if state is not valid.
+  Maybe (Positions, Momenta, Log Double, Log Double)
+leapfrog tF msI l eps q p = do
+  -- The first half step of the momenta.
+  (x, pHalf) <-
+    let (x, pHalf) = leapfrogStepMomenta (0.5 * eps) tF q p
+     in if x > 0.0
+          then Just (x, pHalf)
+          else Nothing
+  -- L-1 full steps for positions and momenta. This gives the positions q_{L-1},
+  -- and the momenta p_{L-1/2}.
+  (qLM1, pLM1Half) <- go (l - 1) $ Just $ (q, pHalf)
+  -- The last full step of the positions.
+  let qL = leapfrogStepPositions msI eps qLM1 pLM1Half
+  -- The last half step of the momenta.
+  (x', pL) <-
+    let (x', pL) = leapfrogStepMomenta (0.5 * eps) tF qL pLM1Half
+     in if x' > 0.0
+          then Just (x', pL)
+          else Nothing
+  return (qL, pL, x, x')
+  where
+    go _ Nothing = Nothing
+    go n (Just (qs, ps))
+      | n <= 0 = Just (qs, ps)
+      | otherwise =
+          let qs' = leapfrogStepPositions msI eps qs ps
+              (x, ps') = leapfrogStepMomenta eps tF qs' p
+           in if x > 0.0
+                then go (n - 1) $ Just $ (qs', ps')
+                else Nothing
+
+leapfrogStepMomenta ::
+  LeapfrogScalingFactor ->
+  Target ->
+  -- Current positions.
+  Positions ->
+  -- Current momenta.
+  Momenta ->
+  -- New momenta; also return value target function to be collected at the end
+  -- of the leapfrog integration.
+  (Log Double, Momenta)
+leapfrogStepMomenta eps tf q p = (x, p + L.scale eps g)
+  where
+    (x, g) = tf q
+
+leapfrogStepPositions ::
+  MassesI ->
+  LeapfrogScalingFactor ->
+  -- Current positions.
+  Positions ->
+  -- Current momenta.
+  Momenta ->
+  -- New positions.
+  Positions
+-- NOTE: Because of numerical errors, the following formulas exhibit different
+-- traces (although the posterior appears to be the same):
+-- 1. This one we cannot use with general matrices:
+--    leapfrogStepPositions msI eps q p = q + (L.scale eps msI L.!#> p)
+-- 2. This one seems to be more numerically unstable:
+--    leapfrogStepPositions msI eps q p = q + L.scale eps (msI L.!#> p)
+leapfrogStepPositions msI eps q p = q + (msI L.!#> L.scale eps p)
diff --git a/src/Mcmc/Proposal/Hamiltonian/Masses.hs b/src/Mcmc/Proposal/Hamiltonian/Masses.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Proposal/Hamiltonian/Masses.hs
@@ -0,0 +1,329 @@
+-- |
+-- Module      :  Mcmc.Proposal.Hamiltonian.Masses
+-- Description :  Mass matrices
+-- Copyright   :  2022 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- Creation date: Tue Jun 14 10:09:24 2022.
+module Mcmc.Proposal.Hamiltonian.Masses
+  ( Mu,
+    MassesI,
+    toGMatrix,
+    cleanMatrix,
+    getMassesI,
+    getMus,
+    Dimension,
+    vectorToMasses,
+    massesToVector,
+
+    -- * Tuning
+    tuneDiagonalMassesOnly,
+    tuneAllMasses,
+  )
+where
+
+import Data.Maybe
+import qualified Data.Vector as VB
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Unboxed as VU
+import Mcmc.Proposal.Hamiltonian.Common
+import qualified Numeric.LinearAlgebra as L
+import qualified Statistics.Covariance as S
+import qualified Statistics.Function as S
+import qualified Statistics.Sample as S
+
+-- Mean vector containing zeroes. We save this vector because it is required
+-- when sampling from the multivariate normal distribution.
+type Mu = L.Vector Double
+
+-- General, symmetric, inverted mass matrix.
+type MassesI = L.GMatrix
+
+-- Purge masses and inverted masses (excluding the diagonal) strictly smaller
+-- than the precision.
+--
+-- If changed, also change help text of 'HTuneMasses', which is indirectly
+-- affected via 'massMin', and 'massMax'.
+precision :: Double
+precision = 1e-8
+
+isDiag :: L.Matrix Double -> Bool
+isDiag xs = abs (sumDiag - sumFull) < precision
+  where
+    xsAbs = L.cmap abs xs
+    sumDiag = L.sumElements (L.takeDiag xsAbs)
+    sumFull = L.sumElements xsAbs
+
+-- Consider a matrix sparse if less than (5 * number of rows) elements are
+-- non-zero.
+isSparse :: L.Matrix Double -> Bool
+isSparse xs = nNonZero < fromIntegral nMax
+  where
+    n = L.rows xs
+    m = min 5 n
+    nMax = n * m
+    f x = if abs x >= precision then 1 else 0 :: Double
+    xsInd = L.cmap f xs
+    nNonZero = L.sumElements xsInd
+
+toAssocMatrix :: L.Matrix Double -> L.AssocMatrix
+toAssocMatrix xs
+  | n /= m = error "toAssocMatrix: Matrix not square."
+  | otherwise =
+      [ ((i, j), e)
+        | i <- [0 .. (n - 1)],
+          j <- [0 .. (n - 1)],
+          let e = xs `L.atIndex` (i, j),
+          abs e >= precision
+      ]
+  where
+    n = L.rows xs
+    m = L.cols xs
+
+toGMatrix :: L.Matrix Double -> L.GMatrix
+toGMatrix xs
+  | n == 0 || m == 0 = error "toGMatrix: Matrix empty."
+  | n /= m = error "toGMatrix: Matrix not square."
+  | isDiag xs = L.mkDiagR n m $ L.takeDiag xs
+  | isSparse xs = L.mkSparse $ toAssocMatrix xs
+  | otherwise = L.mkDense xs
+  where
+    n = L.rows xs
+    m = L.cols xs
+
+-- Diagonal:
+-- - NaN values are set to 1.
+-- - Negative values are set to 1.
+-- - Small elements are set to 'precision'.
+--
+-- Off-diagonal:
+-- - NaN values are set to 0.
+-- - Elements with absolute values strictly smaller than 'precision' are purged.
+--
+-- We are permissive with negative and NaN values because adequate masses are
+-- crucial. The Hamiltonian algorithms also work when the masses are off.
+cleanMatrix :: L.Matrix Double -> L.Matrix Double
+cleanMatrix xs =
+  (L.diag $ L.cmap cleanDiag xsDiag) + L.cmap cleanOffDiag xsOffDiag
+  where
+    xsDiag = L.takeDiag xs
+    cleanDiag x
+      | isNaN x = 1
+      | x < 0 = 1
+      -- The strict comparison is important.
+      | x < precision = precision
+      | otherwise = x
+    xsOffDiag = xs - L.diag xsDiag
+    cleanOffDiag x
+      | isNaN x = 0
+      -- The strict comparison is important.
+      | abs x < precision = 0
+      | otherwise = x
+
+getMassesI :: L.Herm Double -> L.GMatrix
+getMassesI xs
+  | n == 0 || m == 0 = error "getMassesI: Matrix empty."
+  | n /= m = error "getMassesI: Matrix not square."
+  | sign /= 1.0 = error "getMassesI: Determinant of matrix is negative."
+  | otherwise = toGMatrix $ cleanMatrix $ xsI
+  where
+    xs' = L.unSym $ xs
+    n = L.rows xs'
+    m = L.cols xs'
+    (xsI, (_, sign)) = L.invlndet xs'
+
+getMus :: Masses -> L.Vector Double
+getMus xs
+  | n == 0 || m == 0 = error "getMu: Matrix empty."
+  | n /= m = error "getMu: Matrix not square."
+  | otherwise = L.fromList $ replicate n 0.0
+  where
+    xs' = L.unSym xs
+    n = L.rows xs'
+    m = L.cols xs'
+
+-- Dimension of the proposal.
+type Dimension = Int
+
+massesToVector :: Masses -> VU.Vector Double
+massesToVector = VU.convert . L.flatten . L.unSym
+
+vectorToMasses :: Dimension -> VU.Vector Double -> Masses
+vectorToMasses d = L.trustSym . L.reshape d . VU.convert
+
+-- If changed, also change help text of 'HTuneMasses'.
+massMin :: Double
+massMin = precision
+
+-- If changed, also change help text of 'HTuneMasses'.
+massMax :: Double
+massMax = recip precision
+
+-- Minimal number of unique samples required for tuning the diagonal entries of
+-- the mass matrix.
+--
+-- NOTE: If changed, also change help text of 'HTuneMasses'.
+samplesDiagonalMin :: Int
+samplesDiagonalMin = 61
+
+-- Minimal number of samples required for tuning all entries of the mass matrix.
+--
+-- NOTE: If changed, also change help text of 'HTuneMasses'.
+samplesAllMinWith :: Dimension -> Int
+samplesAllMinWith d = samplesDiagonalMin + max samplesDiagonalMin d
+
+getSampleSize :: VS.Vector Double -> Int
+getSampleSize = VS.length . VS.uniq . S.gsort
+
+getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double
+getNewMassDiagonalWithRescue sampleSize massOld massEstimate
+  | sampleSize < samplesDiagonalMin = massOld
+  -- Be permissive with NaN and negative diagonal masses. Diagonal masses are
+  -- variances which are strictly positive.
+  | isNaN massEstimate = massOld
+  | massEstimate <= 0 = massOld
+  | massMin > massNew = massMin
+  | massNew > massMax = massMax
+  | otherwise = massNew
+  where
+    massNewSqrt = recip 3 * (sqrt massOld + 2 * sqrt massEstimate)
+    massNew = massNewSqrt ** 2
+
+-- The Cholesky decomposition, which is performed when sampling new momenta with
+-- 'generateMomenta', requires a positive definite covariance matrix. The
+-- Graphical Lasso algorithm finds positive definite covariance matrices, but
+-- sometimes positive definiteness is violated because of numerical errors.
+-- Further, when non-diagonal masses are already non-zero, the tuning of
+-- diagonal masses only may violate positive definiteness.
+--
+-- Find the closest positive definite matrix of a given matrix.
+--
+-- See https://gist.github.com/fasiha/fdb5cec2054e6f1c6ae35476045a0bbd.
+findClosestPositiveDefiniteMatrix :: L.Matrix Double -> L.Matrix Double
+findClosestPositiveDefiniteMatrix a
+  | n == 0 || m == 0 = error "findClosestPositiveDefiniteMatrix: Matrix empty."
+  | n /= m = error "findClosestPositiveDefiniteMatrix: Matrix not square."
+  | isPositiveDefinite a = a
+  | otherwise = go a3 1
+  where
+    n = L.rows a
+    m = L.cols a
+    b = L.unSym $ L.sym a
+    (_, s, v) = L.svd b
+    h = (L.tr v) L.<> (L.diag s L.<> v)
+    a2 = L.scale 0.5 (b + h)
+    a3 = L.unSym $ L.sym a2
+    isPositiveDefinite = isJust . L.mbChol . L.trustSym
+    --
+    i = L.ident n
+    -- See https://hackage.haskell.org/package/ieee754-0.8.0/docs/src/Numeric-IEEE.html#line-177.
+    eps = 2.2204460492503131e-16
+    go x k
+      | isPositiveDefinite x = x
+      | otherwise =
+          let minEig = L.minElement $ L.cmap L.realPart $ L.eigenvalues x
+              nu = negate minEig * (k ** 2) + eps
+              x' = x + L.scale nu i
+           in go x' (k + 1)
+
+tuneDiagonalMassesOnly ::
+  -- Conversion from value to vector.
+  (a -> Positions) ->
+  -- Value vector.
+  VB.Vector a ->
+  -- Old mass matrix, and inverted mass matrix.
+  (Masses, MassesI) ->
+  -- new mass matrix, and inverted mass matrix.
+  (Masses, MassesI)
+-- NOTE: Here, we lose time because we convert the states to vectors again,
+-- something that has already been done. But then, auto tuning is not a runtime
+-- determining factor.
+tuneDiagonalMassesOnly toVec xs (ms, msI)
+  -- If not enough data is available, do not tune.
+  | VB.length xs < samplesDiagonalMin = (ms, msI)
+  | dimState /= dimMs = error "tuneDiagonalMassesOnly: Dimension mismatch."
+  -- Replace the diagonal.
+  | otherwise =
+      let msDirty = msOld - L.diag msDiagonalOld + L.diag msDiagonalNew
+          -- Positive definite matrices are symmetric.
+          ms' = L.trustSym $ findClosestPositiveDefiniteMatrix $ cleanMatrix msDirty
+          msI' = getMassesI ms'
+       in (ms', msI')
+  where
+    -- xs: Each element contains all parameters of one iteration.
+    -- xs': Each element is a vector containing all parameters changed by the
+    -- proposal of one iteration.
+    xs' = VB.map toVec xs
+    -- xs'': Matrix with each row containing all parameter values changed by the
+    -- proposal of one iteration.
+    xs'' = L.fromRows $ VB.toList xs'
+    -- We can safely use 'VB.head' here since the length of 'xs' must be larger
+    -- than 'samplesDiagonalMin'.
+    dimState = VS.length $ VB.head xs'
+    sampleSizes = VS.fromList $ map getSampleSize $ L.toColumns xs''
+    msOld = L.unSym ms
+    dimMs = L.rows msOld
+    msDiagonalOld = L.takeDiag msOld
+    msDiagonalEstimate = VS.fromList $ map (recip . S.variance) $ L.toColumns xs''
+    msDiagonalNew =
+      VS.zipWith3
+        getNewMassDiagonalWithRescue
+        sampleSizes
+        msDiagonalOld
+        msDiagonalEstimate
+
+-- This value was carefully tuned using the example "hamiltonian".
+defaultGraphicalLassoPenalty :: Double
+defaultGraphicalLassoPenalty = 0.4
+
+tuneAllMasses ::
+  -- Conversion from value to vector.
+  (a -> Positions) ->
+  -- Value vector.
+  VB.Vector a ->
+  -- Old mass matrix, and inverted mass matrix.
+  (Masses, MassesI) ->
+  -- New mass matrix, and inverted mass matrix.
+  (Masses, MassesI)
+-- NOTE: Here, we lose time because we convert the states to vectors again,
+-- something that has already been done. But then, auto tuning is not a runtime
+-- determining factor.
+tuneAllMasses toVec xs (ms, msI)
+  -- If not enough data is available, do not tune.
+  | VB.length xs < samplesDiagonalMin = (ms, msI)
+  -- If not enough data is available, only the diagonal masses are tuned.
+  | VB.length xs < samplesAllMinWith dimMs = fallbackDiagonal
+  | L.rank xs'' /= dimState = fallbackDiagonal
+  | dimState /= dimMs = error "tuneAllMasses: Dimension mismatch."
+  | otherwise = (msPD', msPDI')
+  where
+    fallbackDiagonal = tuneDiagonalMassesOnly toVec xs (ms, msI)
+    -- xs: Each element contains all parameters of one iteration.
+    -- xs': Each element is a vector containing all parameters changed by the
+    -- proposal of one iteration.
+    xs' = VB.map toVec xs
+    -- xs'': Matrix with each row containing all parameter values changed by the
+    -- proposal of one iteration.
+    xs'' = L.fromRows $ VB.toList xs'
+    -- We can safely use 'VB.head' here since the length of 'xs' must be larger
+    -- than 'samplesDiagonalMin'.
+    dimState = VS.length $ VB.head xs'
+    dimMs = L.rows $ L.unSym ms
+    (_, ss, xsNormalized) = S.scale xs''
+    sigmaNormalized =
+      L.unSym $
+        either error fst $
+          S.graphicalLasso defaultGraphicalLassoPenalty xsNormalized
+    -- Sigma is the inverted mass matrix.
+    sigma = S.rescaleSWith ss sigmaNormalized
+    msI' = toGMatrix sigma
+    -- The masses should be positive definite, but sometimes they happen to be
+    -- not because of numerical errors.
+    ms' = L.sym $ cleanMatrix $ L.inv sigma
+    -- Positive definite matrices are symmetric.
+    msPD' = L.trustSym $ findClosestPositiveDefiniteMatrix $ L.unSym ms'
+    msPDI' = if L.unSym ms' == L.unSym msPD' then msI' else getMassesI msPD'
diff --git a/src/Mcmc/Proposal/Hamiltonian/Nuts.hs b/src/Mcmc/Proposal/Hamiltonian/Nuts.hs
new file mode 100644
--- /dev/null
+++ b/src/Mcmc/Proposal/Hamiltonian/Nuts.hs
@@ -0,0 +1,305 @@
+-- |
+-- Module      :  Mcmc.Proposal.Hamiltonian.Nuts
+-- Description :  No-U-Turn sampler (NUTS)
+-- Copyright   :  2022 Dominik Schrempf
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  experimental
+-- Portability :  portable
+--
+-- Creation date: Fri May 27 09:58:23 2022.
+--
+-- For a general introduction to Hamiltonian proposals, see
+-- "Mcmc.Proposal.Hamiltonian.Hamiltonian".
+--
+-- This module implements the No-U-Turn Sampler (NUTS), as described in [4].
+--
+-- Work in progress.
+--
+-- References:
+--
+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using
+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),
+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.
+--
+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data
+--   analysis (2014), CRC Press.
+--
+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual
+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).
+--
+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:
+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of
+--   Machine Learning Research.
+module Mcmc.Proposal.Hamiltonian.Nuts
+  ( NParams (..),
+    defaultNParams,
+    nuts,
+  )
+where
+
+import Data.Bifunctor
+import Mcmc.Acceptance
+import Mcmc.Proposal
+import Mcmc.Proposal.Hamiltonian.Common
+import Mcmc.Proposal.Hamiltonian.Internal
+import Mcmc.Proposal.Hamiltonian.Masses
+import Numeric.AD.Double
+import qualified Numeric.LinearAlgebra as L
+import Numeric.Log
+import System.Random.Stateful
+
+-- Internal; Slice variable 'u'.
+type SliceVariable = Log Double
+
+-- Internal; Forward is True.
+type Direction = Bool
+
+-- Internal; Doubling step number 'j'.
+type DoublingStep = Int
+
+-- Internal; Number of leapfrog steps within the slice 'n'.
+type NStepsOk = Int
+
+-- Internal; Estimated acceptance rate \(\alpha\)'.
+type Alpha = Log Double
+
+-- Internal; Number of accepted steps.
+type NAlpha = Int
+
+-- Internal; Well, that's fun, isn't it? Have a look at Algorithm 3 in [4].
+type BuildTreeReturnType = (Positions, Momenta, Positions, Momenta, Positions, NStepsOk, Alpha, NAlpha)
+
+-- Constant determining largest allowed leapfrog integration error. See
+-- discussion around Equation (3) in [4].
+deltaMax :: Log Double
+deltaMax = Exp 1000
+
+-- Second function in Algorithm 3 and Algorithm 6, respectively in [4].
+buildTreeWith ::
+  -- The exponent of the total energy of the starting state is used to
+  -- calcaulate the expected acceptance rate 'Alpha'.
+  Log Double ->
+  MassesI ->
+  Target ->
+  IOGenM StdGen ->
+  --
+  Positions ->
+  Momenta ->
+  SliceVariable ->
+  Direction ->
+  DoublingStep ->
+  LeapfrogScalingFactor ->
+  IO (Maybe BuildTreeReturnType)
+buildTreeWith expETot0 msI tfun g q p u v j e
+  | j <= 0 =
+      -- Move backwards or forwards?
+      let e' = if v then e else negate e
+       in case leapfrog tfun msI 1 e' q p of
+            Nothing -> pure Nothing
+            Just (q', p', _, expEPot') ->
+              if errorIsSmall
+                then pure $ Just (q', p', q', p', q', n', alpha, 1)
+                else pure Nothing
+              where
+                expEKin' = exponentialKineticEnergy msI p'
+                expETot' = expEPot' * expEKin'
+                n' = if u <= expEPot' * expEKin' then 1 else 0
+                errorIsSmall = u < deltaMax * expETot'
+                alpha' = expETot' / expETot0
+                alpha = min 1.0 alpha'
+
+  -- Recursive case. This is complicated because the algorithm is written for an
+  -- imperative language, and because we have two stacked monads.
+  | otherwise = do
+      mr <- buildTree q p u v (j - 1) e
+      case mr of
+        Nothing -> pure Nothing
+        -- Here, the suffixes 'm' and 'p' stand for minus and plus, respectively.
+        Just (qm, pm, qp, pp, q', n', a', na') -> do
+          mr' <-
+            if v
+              then -- Forwards.
+              do
+                mr'' <- buildTree qp pp u v (j - 1) e
+                case mr'' of
+                  Nothing -> pure Nothing
+                  Just (_, _, qp', pp', q'', n'', a'', na'') ->
+                    pure $ Just (qm, pm, qp', pp', q'', n'', a'', na'')
+              else -- Backwards.
+              do
+                mr'' <- buildTree qm pm u v (j - 1) e
+                case mr'' of
+                  Nothing -> pure Nothing
+                  Just (qm', pm', _, _, q'', n'', a'', na'') ->
+                    pure $ Just (qm', pm', qp, pp, q'', n'', a'', na'')
+          case mr' of
+            Nothing -> pure Nothing
+            Just (qm'', pm'', qp'', pp'', q''', n''', a''', na''') -> do
+              b <- uniformRM (0, 1) g :: IO Double
+              let q'''' = if b < fromIntegral n''' / (fromIntegral $ n' + n''') then q''' else q'
+                  a'''' = a' + a'''
+                  na'''' = na' + na'''
+                  n'''' = n' + n'''
+                  -- Important: Check for U-turn. This formula differs from the
+                  -- formula using indicator functions in Algorithm 3. However,
+                  -- check Equation (4).
+                  isUTurn = let dq = (qp'' - qm'') in (dq * pm'' < 0) || (dq * pp'' < 0)
+              if isUTurn
+                then pure Nothing
+                else pure $ Just (qm'', pm'', qp'', pp'', q'''', n'''', a'''', na'''')
+  where
+    buildTree = buildTreeWith expETot0 msI tfun g
+
+-- | Paramters of the NUTS proposal.
+--
+-- Includes tuning parameters and tuning configuration.
+data NParams = NParams
+  { nLeapfrogScalingFactor :: Maybe LeapfrogScalingFactor,
+    nMasses :: Maybe Masses
+  }
+  deriving (Show)
+
+-- | Default parameters.
+--
+-- - Estimate a reasonable leapfrog scaling factor using Algorithm 4 [4]. If all
+--   fails, use 0.1.
+--
+-- - The mass matrix is set to the identity matrix.
+defaultNParams :: NParams
+defaultNParams = NParams Nothing Nothing
+
+nutsPFunctionWithTuningParameters ::
+  Traversable s =>
+  Dimension ->
+  HStructure s ->
+  (s Double -> Target) ->
+  TuningParameter ->
+  AuxiliaryTuningParameters ->
+  Either String (PFunction (s Double))
+nutsPFunctionWithTuningParameters d hstruct targetWith _ ts = do
+  hParamsI <- fromAuxiliaryTuningParameters d ts
+  pure $ nutsPFunction hParamsI hstruct targetWith
+
+data IsNew
+  = Old
+  | OldWith {_acceptanceCountsOld :: AcceptanceCounts}
+  | NewWith {_acceptanceCountsNew :: AcceptanceCounts}
+
+-- First function in Algorithm 3.
+nutsPFunction ::
+  HParamsI ->
+  HStructure s ->
+  (s Double -> Target) ->
+  PFunction (s Double)
+nutsPFunction hparamsi hstruct targetWith x g = do
+  p <- generateMomenta mus ms g
+  uZeroOne <- uniformRM (0, 1) g :: IO Double
+  -- NOTE (runtime): Here we need the target function value from the previous
+  -- step. For now, I just recalculate the value, but this is, of course, slow!
+  -- However, if other proposals have changed the state inbetween, we do need to
+  -- recalculate this value.
+  let q = toVec x
+      expEPot = fst $ target q
+      expEKin = exponentialKineticEnergy msI p
+      expETot = expEPot * expEKin
+      uZeroOneL = Exp $ log uZeroOne
+      u = expETot * uZeroOneL
+  let -- Recursive case. This is complicated because the algorithm is written for an
+      -- imperative language, and because we have two stacked monads.
+      --
+      -- Here, the suffixes 'm' and 'p' stand for minus and plus, respectively.
+      go qm pm qp pp j y n isNew = do
+        v <- uniformM g :: IO Direction
+        mr' <-
+          if v
+            then -- Forwards.
+            do
+              mr <- buildTreeWith expETot msI target g qp pp u v j e
+              case mr of
+                Nothing -> pure Nothing
+                Just (_, _, qp', pp', y', n', a, na) -> pure $ Just (qm, pm, qp', pp', y', n', a, na)
+            else -- Backwards.
+            do
+              mr <- buildTreeWith expETot msI target g qm pm u v j e
+              case mr of
+                Nothing -> pure Nothing
+                Just (qm', pm', _, _, y', n', a, na) -> pure $ Just (qm', pm', qp, pp, y', n', a, na)
+        case mr' of
+          Nothing -> pure (y, isNew)
+          Just (qm'', pm'', qp'', pp'', y'', n'', a, na) -> do
+            let r = fromIntegral n'' / fromIntegral n :: Double
+                ar = (exp $ ln a) / fromIntegral na :: Double
+                getCounts s = max 0 $ min 100 $ round $ s * 100
+                ac =
+                  if ar >= 0
+                    then let cs = getCounts ar in AcceptanceCounts cs (100 - cs)
+                    else error $ "nutsPFunction: Acceptance rate negative."
+            isAccept <-
+              if r > 1.0
+                then pure True
+                else do
+                  b <- uniformRM (0, 1) g
+                  pure $ b < r
+            let (y''', isNew') = if isAccept then (y'', NewWith ac) else (y, OldWith ac)
+                isUTurn = let dq = (qp'' - qm'') in (dq * pm'' < 0) || (dq * pp'' < 0)
+            if isUTurn
+              then pure (y''', isNew')
+              else go qm'' pm'' qp'' pp'' (j + 1) y''' (n + n'') isNew'
+  (x', isNew) <- go q p q p 0 q 1 Old
+  pure $ case isNew of
+    Old -> (ForceReject, Just $ AcceptanceCounts 0 100)
+    OldWith ac -> (ForceReject, Just $ ac)
+    NewWith ac -> (ForceAccept $ fromVec x', Just $ ac)
+  where
+    (HParamsI e _ ms _ _ msI mus) = hparamsi
+    (HStructure _ toVec fromVecWith) = hstruct
+    fromVec = fromVecWith x
+    target = targetWith x
+
+-- | No U-turn Hamiltonian Monte Carlo sampler (NUTS).
+--
+-- The structure of the state is denoted as @s@.
+--
+-- May call 'error' during initialization.
+nuts ::
+  Traversable s =>
+  NParams ->
+  HTuningConf ->
+  HStructure s ->
+  HTarget s ->
+  PName ->
+  PWeight ->
+  Proposal (s Double)
+nuts nparams htconf hstruct htarget n w =
+  let -- Misc.
+      desc = PDescription "No U-turn sampler (NUTS)"
+      (HStructure sample toVec fromVec) = hstruct
+      dim = L.size $ toVec sample
+      -- See bottom of page 1616 in [4].
+      pDim = PSpecial dim 0.6
+      -- Vectorize and derive the target function.
+      (HTarget mPrF lhF mJcF) = htarget
+      tF y = case (mPrF, mJcF) of
+        (Nothing, Nothing) -> lhF y
+        (Just prF, Nothing) -> prF y * lhF y
+        (Nothing, Just jcF) -> lhF y * jcF y
+        (Just prF, Just jcF) -> prF y * lhF y * jcF y
+      tFnG = grad' (ln . tF)
+      targetWith x = bimap Exp toVec . tFnG . fromVec x
+      (NParams mEps mMs) = nparams
+      hParamsI =
+        either error id $
+          hParamsIWith (targetWith sample) (toVec sample) mEps Nothing mMs
+      ps = nutsPFunction hParamsI hstruct targetWith
+      nutsWith = Proposal n desc PSlow pDim w ps
+      -- Tuning.
+      ts = toAuxiliaryTuningParameters hParamsI
+      tuner = do
+        tfun <- hTuningFunctionWith dim toVec htconf
+        let pfun = nutsPFunctionWithTuningParameters dim hstruct targetWith
+        pure $ Tuner 1.0 ts True tfun pfun
+   in case checkHStructureWith (hpsMasses hParamsI) hstruct of
+        Just err -> error err
+        Nothing -> nutsWith tuner
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
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Proposal.Scale
 -- Description :  Multiplicative proposals
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -26,8 +26,8 @@
 
 -- The actual proposal with tuning parameter. The tuning parameter does not
 -- change the mean.
-scaleSimple :: Shape Double -> Scale Double -> TuningParameter -> ProposalSimple Double
-scaleSimple k th t =
+scalePFunction :: Shape Double -> Scale Double -> TuningParameter -> PFunction Double
+scalePFunction k th t =
   genericContinuous
     (gammaDistr (k / t) (th * t))
     (*)
@@ -52,7 +52,7 @@
   PWeight ->
   Tune ->
   Proposal Double
-scale k th = createProposal description (scaleSimple k th) PFast (PDimension 1)
+scale k th = createProposal description (scalePFunction k th) PFast (PDimension 1)
   where
     description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th
 
@@ -66,16 +66,16 @@
   PWeight ->
   Tune ->
   Proposal Double
-scaleUnbiased k = createProposal description (scaleSimple k (1 / k)) PFast (PDimension 1)
+scaleUnbiased k = createProposal description (scalePFunction k (1 / k)) PFast (PDimension 1)
   where
     description = PDescription $ "Scale unbiased; shape: " ++ show k
 
-scaleContrarilySimple ::
+scaleContrarilyPFunction ::
   Shape Double ->
   Scale Double ->
   TuningParameter ->
-  ProposalSimple (Double, Double)
-scaleContrarilySimple k th t =
+  PFunction (Double, Double)
+scaleContrarilyPFunction k th t =
   genericContinuous
     (gammaDistr (k / t) (th * t))
     contra
@@ -96,6 +96,6 @@
   PWeight ->
   Tune ->
   Proposal (Double, Double)
-scaleContrarily k th = createProposal description (scaleContrarilySimple k th) PFast (PDimension 2)
+scaleContrarily k th = createProposal description (scaleContrarilyPFunction k th) PFast (PDimension 2)
   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
--- a/src/Mcmc/Proposal/Simplex.hs
+++ b/src/Mcmc/Proposal/Simplex.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Proposal.Simplex
 -- Description :  Proposals on simplices
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -33,8 +33,6 @@
 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.
@@ -100,16 +98,14 @@
 -- 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 :: TuningParameter -> ProposalSimple Simplex
-dirichletSimple t (SimplexUnsafe xs) g = do
+dirichletPFunction :: TuningParameter -> PFunction Simplex
+dirichletPFunction 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
@@ -118,9 +114,7 @@
         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)
+  pure (Propose (SimplexUnsafe ys) r 1.0, Nothing)
   where
     tf = getTuningFunction t
 
@@ -140,7 +134,7 @@
 -- For high dimensional simplices, this proposal may have low acceptance rates.
 -- In this case, please see the coordinate wise 'beta' proposal.
 dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex
-dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple PFast
+dirichlet = createProposal (PDescription "Dirichlet") dirichletPFunction PFast
 
 -- The tuning parameter is the inverted mean of the shape values.
 --
@@ -149,8 +143,8 @@
 -- analyses.
 --
 -- See also the 'dirichlet' proposal.
-betaSimple :: Dimension -> TuningParameter -> ProposalSimple Simplex
-betaSimple i t (SimplexUnsafe xs) g = do
+betaPFunction :: Dimension -> TuningParameter -> PFunction Simplex
+betaPFunction 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.
@@ -179,7 +173,8 @@
       -- 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)
+      y = either error id $ simplexFromVector ys
+  pure (Propose y r jac, Nothing)
   where
     xI = xs V.! i
     xsSum = V.sum xs
@@ -206,6 +201,6 @@
 -- This proposal has been assigned a dimension of 2. See the discussion at
 -- 'PDimension'.
 beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex
-beta i = createProposal description (betaSimple i) PFast (PDimension 2)
+beta i = createProposal description (betaPFunction i) PFast (PDimension 2)
   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
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  Mcmc.Proposal.Slide
 -- Description :  Additive proposals
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -26,8 +26,8 @@
 import Statistics.Distribution.Uniform
 
 -- The actual proposal with tuning parameter.
-slideSimple :: Mean Double -> StandardDeviation Double -> TuningParameter -> ProposalSimple Double
-slideSimple m s t =
+slidePFunction :: Mean Double -> StandardDeviation Double -> TuningParameter -> PFunction Double
+slidePFunction m s t =
   genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing
 
 -- | Additive proposal.
@@ -40,13 +40,13 @@
   PWeight ->
   Tune ->
   Proposal Double
-slide m s = createProposal description (slideSimple m s) PFast (PDimension 1)
+slide m s = createProposal description (slidePFunction m s) PFast (PDimension 1)
   where
     description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s
 
 -- The actual proposal with tuning parameter.
-slideSymmetricSimple :: StandardDeviation Double -> TuningParameter -> ProposalSimple Double
-slideSymmetricSimple s t =
+slideSymmetricPFunction :: StandardDeviation Double -> TuningParameter -> PFunction Double
+slideSymmetricPFunction s t =
   genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing
 
 -- | See 'slide'.
@@ -60,13 +60,13 @@
   PWeight ->
   Tune ->
   Proposal Double
-slideSymmetric s = createProposal description (slideSymmetricSimple s) PFast (PDimension 1)
+slideSymmetric s = createProposal description (slideSymmetricPFunction s) PFast (PDimension 1)
   where
     description = PDescription $ "Slide symmetric; sd: " ++ show s
 
 -- The actual proposal with tuning parameter.
-slideUniformSimple :: Size -> TuningParameter -> ProposalSimple Double
-slideUniformSimple d t =
+slideUniformPFunction :: Size -> TuningParameter -> PFunction Double
+slideUniformPFunction d t =
   genericContinuous (uniformDistr (-t * d) (t * d)) (+) Nothing Nothing
 
 -- | See 'slide'.
@@ -80,19 +80,19 @@
   PWeight ->
   Tune ->
   Proposal Double
-slideUniformSymmetric d = createProposal description (slideUniformSimple d) PFast (PDimension 1)
+slideUniformSymmetric d = createProposal description (slideUniformPFunction d) PFast (PDimension 1)
   where
     description = PDescription $ "Slide uniform symmetric; delta: " ++ show d
 
 contra :: (Double, Double) -> Double -> (Double, Double)
 contra (x, y) u = (x + u, y - u)
 
-slideContrarilySimple ::
+slideContrarilyPFunction ::
   Mean Double ->
   StandardDeviation Double ->
   TuningParameter ->
-  ProposalSimple (Double, Double)
-slideContrarilySimple m s t =
+  PFunction (Double, Double)
+slideContrarilyPFunction m s t =
   genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing
 
 -- | See 'slide'.
@@ -108,6 +108,6 @@
   PWeight ->
   Tune ->
   Proposal (Double, Double)
-slideContrarily m s = createProposal description (slideContrarilySimple m s) PFast (PDimension 2)
+slideContrarily m s = createProposal description (slideContrarilyPFunction m s) PFast (PDimension 2)
   where
     description = PDescription $ "Slide contrarily; mean: " ++ show m ++ ", sd: " ++ show s
diff --git a/src/Mcmc/Settings.hs b/src/Mcmc/Settings.hs
--- a/src/Mcmc/Settings.hs
+++ b/src/Mcmc/Settings.hs
@@ -5,7 +5,7 @@
 -- |
 -- Module      :  Mcmc.Settings
 -- Description :  Settings of Markov chain Monte Carlo samplers
--- Copyright   :  (c) Dominik Schrempf, 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/Mcmc/Statistics/Types.hs b/src/Mcmc/Statistics/Types.hs
--- a/src/Mcmc/Statistics/Types.hs
+++ b/src/Mcmc/Statistics/Types.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Mcmc.Statistics.Types
 -- Description :  Types indicating properties of distributions
--- Copyright   :  (c) 2021 Dominik Schrempf
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/test/Mcmc/ProposalSpec.hs b/test/Mcmc/ProposalSpec.hs
--- a/test/Mcmc/ProposalSpec.hs
+++ b/test/Mcmc/ProposalSpec.hs
@@ -1,7 +1,7 @@
 -- |
 --   Module      :  Mcmc.ProposalSpec
 --   Description :  Unit tests for Mcmc.Proposal
---   Copyright   :  (c) Dominik Schrempf, 2021
+--   Copyright   :  2021 Dominik Schrempf
 --   License     :  GPL-3.0-or-later
 --
 --   Maintainer  :  dominik.schrempf@gmail.com
@@ -17,7 +17,7 @@
 import Mcmc.Cycle
 import Mcmc.Proposal
 import Mcmc.Proposal.Slide
-import System.Random.MWC
+import System.Random.Stateful
 import Test.Hspec
 
 p1 :: Proposal Double
@@ -34,7 +34,7 @@
   describe "prepareProposals" $
     it "returns the correct number of proposals in a cycle" $
       do
-        g <- create
+        g <- newIOGenM $ mkStdGen 0
         l1 <- length <$> prepareProposals AllProposals c g
         l1 `shouldBe` 4
         l2 <- length <$> prepareProposals AllProposals (setOrder RandomReversibleO c) g
diff --git a/test/Mcmc/SaveSpec.hs b/test/Mcmc/SaveSpec.hs
--- a/test/Mcmc/SaveSpec.hs
+++ b/test/Mcmc/SaveSpec.hs
@@ -1,7 +1,7 @@
 -- |
 --   Module      :  Mcmc.SaveSpec
 --   Description :  Unit tests for Mcmc.Save
---   Copyright   :  (c) Dominik Schrempf, 2021
+--   Copyright   :  2021 Dominik Schrempf
 --   License     :  GPL-3.0-or-later
 --
 --   Maintainer  :  dominik.schrempf@gmail.com
@@ -14,13 +14,14 @@
   )
 where
 
+import Data.IORef (readIORef)
 import Mcmc
 import Mcmc.Chain.Chain
 import Mcmc.Chain.Save
 import Mcmc.Chain.Trace
 import Statistics.Distribution
 import Statistics.Distribution.Normal
-import qualified System.Random.MWC as R
+import System.Random.Stateful
 import Test.Hspec
 
 trueMean :: Double
@@ -65,7 +66,7 @@
   describe "save and load" $
     it "doesn't change the MCMC chain" $
       do
-        gen <- R.create
+        gen <- newIOGenM $ mkStdGen 0
         a <- mhg settings noPrior lh proposals mon 0 gen
         c <- fromMHG <$> mcmc settings a
         savedChain <- toSavedChain c
@@ -87,17 +88,17 @@
         frozenT2 <- freezeT (trace c)
         frozenT2' <- freezeT (trace c')
         frozenT2 `shouldBe` frozenT2'
-        g2 <- R.save $ generator r
-        g2' <- R.save $ generator r'
+        g2 <- readIORef $ unIOGenM $ generator r
+        g2' <- readIORef $ unIOGenM $ generator r'
         g2 `shouldBe` g2'
 
   describe "mhContinue" $
     it "mcmc 50 + mcmcContinue 50 == mcmc 100" $
       do
-        gen1 <- R.create
+        gen1 <- newIOGenM $ mkStdGen 0
         a1 <- mhg settings noPrior lh proposals mon 0 gen1
         r1 <- fromMHG <$> mcmc settings a1
-        gen2 <- R.create
+        gen2 <- newIOGenM $ mkStdGen 0
         let settings' = settings {sIterations = Iterations 50}
         a2 <- mhg settings' noPrior lh proposals mon 0 gen2
         r2' <- mcmc settings' a2
@@ -107,6 +108,6 @@
         frozenT1 <- freezeT (trace r1)
         frozenT2 <- freezeT (trace r2)
         frozenT1 `shouldBe` frozenT2
-        g <- R.save $ generator r1
-        g' <- R.save $ generator r2
+        g <- readIORef $ unIOGenM $ generator r1
+        g' <- readIORef $ unIOGenM $ generator r2
         g `shouldBe` g'