diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -5,6 +5,14 @@
 ## Unreleased changes
 
 
+## 0.8.1.0
+
+-   Automatic intermediate tuning for HMC and NUTS.
+-   Fix documentation for generic proposals.
+-   Tooling updates.
+-   Slight runtime improvements (strictness annotations).
+
+
 ## 0.8.0.1
 
 -   Improve exception handling (also during execution of monitors; also improve
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -22,7 +22,7 @@
 import System.Random.Stateful
 
 gammaBenchG :: RealFloat a => (a -> a) -> [a] -> a
-gammaBenchG f = foldl' (\acc x -> acc + (f x)) 0
+gammaBenchG f = foldl' (\acc x -> acc + f x) 0
 {-# SPECIALIZE gammaBenchG :: (Double -> Double) -> [Double] -> Double #-}
 
 gammaVals :: [Double]
diff --git a/bench/Poisson.hs b/bench/Poisson.hs
--- a/bench/Poisson.hs
+++ b/bench/Poisson.hs
@@ -1,7 +1,6 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE TypeSynonymInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# OPTIONS_GHC -Wno-orphans #-}
 
@@ -60,10 +59,10 @@
 lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]
 
 proposalAlpha :: Proposal I
-proposalAlpha = (VB.element 0) @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune
+proposalAlpha = VB.element 0 @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune
 
 proposalBeta :: Proposal I
-proposalBeta = (VB.element 1) @~ 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]
@@ -103,7 +102,7 @@
   void $ mcmc s a
 
 toVec :: I -> VS.Vector Double
-toVec xs = VS.generate 2 (\i -> xs VB.! i)
+toVec xs = VS.generate 2 (xs VB.!)
 
 fromVec :: I -> VS.Vector Double -> I
 fromVec _ xs = VB.mk2 (xs VS.! 0) (xs VS.! 1)
diff --git a/mcmc.cabal b/mcmc.cabal
--- a/mcmc.cabal
+++ b/mcmc.cabal
@@ -1,6 +1,6 @@
 cabal-version:      3.0
 name:               mcmc
-version:            0.8.0.1
+version:            0.8.1.0
 synopsis:           Sample from a posterior using Markov chain Monte Carlo
 description:
   Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme>
diff --git a/src/Mcmc/Acceptance.hs b/src/Mcmc/Acceptance.hs
--- a/src/Mcmc/Acceptance.hs
+++ b/src/Mcmc/Acceptance.hs
@@ -15,11 +15,13 @@
   ( -- * Acceptance rates
     AcceptanceRate,
     AcceptanceCounts (..),
-    Acceptance (fromAcceptance),
+    AcceptanceRates (..),
+    Acceptance,
+    Acceptances (fromAcceptances),
     emptyA,
     pushAccept,
     pushReject,
-    pushAcceptanceCounts,
+    ResetAcceptance (..),
     resetA,
     transformKeysA,
     acceptanceRate,
@@ -44,48 +46,74 @@
 
 $(deriveJSON defaultOptions ''AcceptanceCounts)
 
-addAccept :: AcceptanceCounts -> AcceptanceCounts
-addAccept (AcceptanceCounts a r) = AcceptanceCounts (a + 1) r
+-- | Proposals based on Hamiltonian dynamics use expected acceptance rates, not counts.
+data AcceptanceRates = AcceptanceRates
+  { totalAcceptanceRate :: !Double,
+    nAcceptanceRates :: !Int
+  }
+  deriving (Show, Eq)
 
-addReject :: AcceptanceCounts -> AcceptanceCounts
-addReject (AcceptanceCounts a r) = AcceptanceCounts a (r + 1)
+$(deriveJSON defaultOptions ''AcceptanceRates)
 
-addAcceptanceCounts :: AcceptanceCounts -> AcceptanceCounts -> AcceptanceCounts
-addAcceptanceCounts (AcceptanceCounts al rl) (AcceptanceCounts ar rr) =
-  AcceptanceCounts (al + ar) (rl + rr)
+-- | Stored actual acceptance counts and maybe expected acceptance rates.
+data Acceptance = A AcceptanceCounts (Maybe AcceptanceRates)
+  deriving (Show, Eq)
 
+$(deriveJSON defaultOptions ''Acceptance)
+
+addAccept :: Maybe AcceptanceRates -> Acceptance -> Acceptance
+addAccept mr' (A (AcceptanceCounts a r) mr) = A (AcceptanceCounts (a + 1) r) (addAcceptanceRates mr' mr)
+
+addReject :: Maybe AcceptanceRates -> Acceptance -> Acceptance
+addReject mr' (A (AcceptanceCounts a r) mr) = A (AcceptanceCounts a (r + 1)) (addAcceptanceRates mr' mr)
+
+addAcceptanceRates :: Maybe AcceptanceRates -> Maybe AcceptanceRates -> Maybe AcceptanceRates
+addAcceptanceRates Nothing Nothing = Nothing
+addAcceptanceRates (Just r) Nothing = Just r
+addAcceptanceRates Nothing (Just r) = Just r
+addAcceptanceRates (Just (AcceptanceRates al rl)) (Just (AcceptanceRates ar rr)) =
+  Just $ AcceptanceRates (al + ar) (rl + rr)
+
 -- | For each key @k@, store the number of accepted and rejected proposals.
-newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k AcceptanceCounts}
+newtype Acceptances k = Acceptances {fromAcceptances :: M.Map k Acceptance}
   deriving (Eq, Show)
 
-instance ToJSONKey k => ToJSON (Acceptance k) where
-  toJSON (Acceptance m) = toJSON m
-  toEncoding (Acceptance m) = toEncoding m
+instance ToJSONKey k => ToJSON (Acceptances k) where
+  toJSON (Acceptances m) = toJSON m
+  toEncoding (Acceptances m) = toEncoding m
 
-instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where
-  parseJSON v = Acceptance <$> parseJSON v
+instance (Ord k, FromJSONKey k) => FromJSON (Acceptances k) where
+  parseJSON v = Acceptances <$> parseJSON v
 
 -- | In the beginning there was the Word.
 --
 -- Initialize an empty storage of accepted/rejected values.
-emptyA :: Ord k => [k] -> Acceptance k
-emptyA ks = Acceptance $ M.fromList [(k, AcceptanceCounts 0 0) | k <- ks]
+emptyA :: Ord k => [k] -> Acceptances k
+emptyA ks = Acceptances $ M.fromList [(k, A noCounts Nothing) | k <- ks]
+  where
+    noCounts = AcceptanceCounts 0 0
 
 -- | For key @k@, add an accept.
-pushAccept :: Ord k => k -> Acceptance k -> Acceptance k
-pushAccept k = Acceptance . M.adjust addAccept k . fromAcceptance
+pushAccept :: Ord k => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k
+pushAccept mr k = Acceptances . M.adjust (addAccept mr) k . fromAcceptances
 
 -- | For key @k@, add a reject.
-pushReject :: Ord k => k -> Acceptance k -> Acceptance k
-pushReject k = Acceptance . M.adjust addReject k . fromAcceptance
+pushReject :: Ord k => Maybe AcceptanceRates -> k -> Acceptances k -> Acceptances k
+pushReject mr k = Acceptances . M.adjust (addReject mr) k . fromAcceptances
 
--- | 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 specification.
+data ResetAcceptance
+  = -- | Reset actual acceptance counts and expected acceptance rates.
+    ResetEverything
+  | -- | Only reset expected acceptance rates.
+    ResetExpectedRatesOnly
 
 -- | Reset acceptance counts.
-resetA :: Ord k => Acceptance k -> Acceptance k
-resetA = emptyA . M.keys . fromAcceptance
+resetA :: Ord k => ResetAcceptance -> Acceptances k -> Acceptances k
+resetA ResetEverything = emptyA . M.keys . fromAcceptances
+resetA ResetExpectedRatesOnly = Acceptances . M.map f . fromAcceptances
+  where
+    f (A cs _) = A cs Nothing
 
 transformKeys :: (Ord k1, Ord k2) => [(k1, k2)] -> M.Map k1 v -> M.Map k2 v
 transformKeys ks m = foldl' insrt M.empty ks
@@ -94,10 +122,11 @@
 
 -- | Transform keys using the given lists. Keys not provided will not be present
 -- in the new 'Acceptance' variable.
-transformKeysA :: (Ord k1, Ord k2) => [(k1, k2)] -> Acceptance k1 -> Acceptance k2
-transformKeysA ks = Acceptance . transformKeys ks . fromAcceptance
+transformKeysA :: (Ord k1, Ord k2) => [(k1, k2)] -> Acceptances k1 -> Acceptances k2
+transformKeysA ks = Acceptances . transformKeys ks . fromAcceptances
 
--- | Acceptance counts and rate for a specific proposal.
+-- | Compute acceptance counts, and actual and expected acceptances rates for a
+-- specific proposal.
 --
 -- Return @Just (accepts, rejects, acceptance rate)@.
 --
@@ -106,23 +135,28 @@
 acceptanceRate ::
   Ord k =>
   k ->
-  Acceptance k ->
-  Maybe (Int, Int, AcceptanceRate)
-acceptanceRate k a = case fromAcceptance a M.!? k of
-  Just (AcceptanceCounts 0 0) -> Nothing
-  Just (AcceptanceCounts as rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))
+  Acceptances k ->
+  -- | (nAccepts, nRejects, actualRate, expectedRate)
+  (Int, Int, Maybe AcceptanceRate, Maybe AcceptanceRate)
+acceptanceRate k a = case fromAcceptances a M.!? k of
+  Just (A (AcceptanceCounts as rs) mrs) -> (as, rs, mar, mtr)
+    where
+      s = as + rs
+      mar = if s <= 0 then Nothing else Just $ fromIntegral as / fromIntegral s
+      mtr = case mrs of
+        Nothing -> Nothing
+        Just (AcceptanceRates xs n) -> Just $ xs / fromIntegral n
   Nothing -> error "acceptanceRate: Key not found in map."
 
--- | Acceptance rates for all proposals.
+-- | Compute actual acceptance rates for all proposals.
 --
 -- Set rate to 'Nothing' if no proposals have been accepted or rejected
 -- (division by zero).
-acceptanceRates :: Acceptance k -> M.Map k (Maybe AcceptanceRate)
-acceptanceRates =
-  M.map
-    ( \(AcceptanceCounts as rs) ->
-        if as + rs == 0
-          then Nothing
-          else Just $ fromIntegral as / fromIntegral (as + rs)
-    )
-    . fromAcceptance
+acceptanceRates :: Acceptances k -> M.Map k (Maybe AcceptanceRate)
+acceptanceRates = M.map getRate . fromAcceptances
+  where
+    getRate (A (AcceptanceCounts as rs) _) =
+      let s = as + rs
+       in if s <= 0
+            then Nothing
+            else Just $ fromIntegral as / fromIntegral s
diff --git a/src/Mcmc/Algorithm.hs b/src/Mcmc/Algorithm.hs
--- a/src/Mcmc/Algorithm.hs
+++ b/src/Mcmc/Algorithm.hs
@@ -16,6 +16,7 @@
 
 import qualified Data.ByteString.Lazy.Char8 as BL
 import Data.Time
+import Mcmc.Acceptance
 import Mcmc.Cycle
 import Mcmc.Proposal
 import Mcmc.Settings
@@ -42,7 +43,7 @@
   aAutoTune :: TuningType -> Int -> a -> IO a
 
   -- | Reset acceptance counts.
-  aResetAcceptance :: a -> a
+  aResetAcceptance :: ResetAcceptance -> a -> a
 
   -- | Clean after burn in. In particular, this is used to reduce the length of
   -- the trace, if required.
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
@@ -127,7 +127,7 @@
     savedMC3Chains :: V.Vector (SavedChain a),
     savedMC3ReciprocalTemperatures :: ReciprocalTemperatures,
     savedMC3Iteration :: Int,
-    savedMC3SwapAcceptance :: Acceptance Int,
+    savedMC3SwapAcceptance :: Acceptances Int,
     savedMC3Generator :: (Word64, Word64)
   }
   deriving (Eq, Show)
@@ -172,7 +172,7 @@
     -- | Current iteration.
     mc3Iteration :: Int,
     -- | Number of accepted and rejected swaps.
-    mc3SwapAcceptance :: Acceptance Int,
+    mc3SwapAcceptances :: Acceptances Int,
     mc3Generator :: IOGenM StdGen
   }
 
@@ -252,7 +252,7 @@
 initMHG prf lhf i beta a
   | i < 0 = error "initMHG: Chain index negative."
   -- Do nothing for the cold chain.
-  | i == 0 = return $ MHG $ c
+  | i == 0 = return $ MHG c
   | otherwise = do
       -- We have to push the current link in the trace, since it is not set by
       -- 'setReciprocalTemperature'. The other links in the trace are still
@@ -411,11 +411,11 @@
       -- 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' = pushAccept i (mc3SwapAcceptance a)
-      return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'}
+      let !ac' = pushAccept Nothing i (mc3SwapAcceptances a)
+      return $ a {mc3MHGChains = y, mc3SwapAcceptances = ac'}
     else do
-      let !ac' = pushReject i (mc3SwapAcceptance a)
-      return $ a {mc3SwapAcceptance = ac'}
+      let !ac' = pushReject Nothing i (mc3SwapAcceptances a)
+      return $ a {mc3SwapAcceptances = ac'}
   where
     g = mc3Generator a
 
@@ -483,7 +483,7 @@
   mhgs' <- V.mapM (aAutoTune b l) $ mc3MHGChains a
   -- 2. Auto tune temperatures.
   let optimalRate = getOptimalRate PDimensionUnknown
-      mCurrentRates = acceptanceRates $ mc3SwapAcceptance a
+      mCurrentRates = acceptanceRates $ mc3SwapAcceptances a
       -- We assume that the acceptance rate of state swaps between two chains is
       -- roughly proportional to the ratio of the temperatures of the chains.
       -- Hence, we focus on temperature ratios, actually reciprocal temperature
@@ -512,15 +512,15 @@
             (V.tail mhgs')
   return $ a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'}
 
-mc3ResetAcceptance :: ToJSON a => MC3 a -> MC3 a
-mc3ResetAcceptance a = a'
+mc3ResetAcceptance :: ToJSON a => ResetAcceptance -> MC3 a -> MC3 a
+mc3ResetAcceptance x a = a'
   where
     -- 1. Reset acceptance of all chains.
-    mhgs' = V.map aResetAcceptance (mc3MHGChains a)
+    mhgs' = V.map (aResetAcceptance x) (mc3MHGChains a)
     -- 2. Reset acceptance of swaps.
-    ac' = resetA $ mc3SwapAcceptance a
+    ac' = resetA x $ mc3SwapAcceptances a
     --
-    a' = a {mc3MHGChains = mhgs', mc3SwapAcceptance = ac'}
+    a' = a {mc3MHGChains = mhgs', mc3SwapAcceptances = ac'}
 
 mc3CleanAfterBurnIn :: ToJSON a => TraceLength -> MC3 a -> IO (MC3 a)
 mc3CleanAfterBurnIn tl a = do
@@ -573,14 +573,14 @@
     -- Acceptance rates may be 'Nothing' when no proposals have been undertaken.
     -- The 'sequence' operations pull the 'Nothing's out of the inner
     -- structures.
-    as = sequence $ V.map (sequence . acceptanceRates . acceptance) cs
+    as = sequence $ V.map (sequence . acceptanceRates . acceptances) cs
     mVecAr = V.map (\mp -> sum mp / fromIntegral (length mp)) <$> as
     mAr = (\vec -> V.sum vec / fromIntegral (V.length vec)) <$> mVecAr
     bs = mc3ReciprocalTemperatures a
     bsB = map (BB.toLazyByteString . BB.formatDouble (BB.standard 2)) $ U.toList bs
     swapPeriod = fromSwapPeriod $ mc3SwapPeriod $ mc3Settings a
     swapPeriodB = BB.toLazyByteString $ BB.intDec swapPeriod
-    swapAcceptance = mc3SwapAcceptance a
+    swapAcceptance = mc3SwapAcceptances a
     n = fromNChains $ mc3NChains $ mc3Settings a
     proposalHLine = BL.replicate (BL.length proposalHeader) '-'
 
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
@@ -233,7 +233,7 @@
 mhgPropose :: MHG a -> Proposal a -> IO (MHG a)
 mhgPropose (MHG c) p = do
   -- 1. Sample new state.
-  !(pres, mcs) <- liftIO $ s x g
+  (!pres, !mcs) <- liftIO $ s x g
   -- 2. Define new prior and likelihood calculation functions. Avoid actual
   -- calculation of the values.
   --
@@ -242,15 +242,11 @@
   -- https://stackoverflow.com/a/46603680/3536806.
   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'}
+        let !ac' = pushAccept mcs p ac
+         in pure $ MHG $ c {link = Link y pr lh, acceptances = ac'}
       reject =
-        let !ac' = case mcs of
-              Nothing -> pushReject p ac
-              Just cs -> pushAcceptanceCounts p cs ac
-         in pure $ MHG $ c {acceptance = ac'}
+        let !ac' = pushReject mcs p ac
+         in pure $ MHG $ c {acceptances = ac'}
   -- 3. Accept or reject.
   --
   -- 3a. When rejection is inevitable, avoid calculation of the prior, the
@@ -274,7 +270,7 @@
     (Link x pX lX) = link c
     pF = priorFunction c
     lF = likelihoodFunction c
-    ac = acceptance c
+    ac = acceptances c
     g = generator c
 
 mhgPush :: MHG a -> IO (MHG a)
@@ -314,21 +310,31 @@
     g = generator c
 
 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}
+mhgAutoTune tt n (MHG c)
+  | isIntermediate =
+      pure . MHG $
+        if ccHasIntermediateTuners cc
+          then -- Do not provide trace when tuning intermediately.
+            c {cycle = autoTuneCycle tt ac Nothing cc}
+          else -- Skip intermediate tuning completely when unnecessary.
+            c
+  | otherwise = do
+      mxs <-
+        -- Provide the trace if required.
+        if ccRequireTrace cc
+          then Just . VB.map state <$> takeT n tr
+          else pure Nothing
+      pure $ MHG c {cycle = autoTuneCycle tt ac mxs cc}
   where
-    ac = acceptance c
+    isIntermediate = tt == IntermediateTuningFastProposalsOnly || tt == IntermediateTuningAllProposals
+    ac = acceptances c
     cc = cycle c
     tr = trace c
 
-mhgResetAcceptance :: MHG a -> MHG a
-mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}
+mhgResetAcceptance :: ResetAcceptance -> MHG a -> MHG a
+mhgResetAcceptance a (MHG c) = MHG $ c {acceptances = resetA a ac}
   where
-    ac = acceptance c
+    ac = acceptances c
 
 mhgCleanAfterBurnIn :: TraceLength -> MHG a -> IO (MHG a)
 mhgCleanAfterBurnIn tl (MHG c) = do
@@ -346,7 +352,7 @@
 mhgSummarizeCycle m (MHG c) = summarizeCycle m ac cc
   where
     cc = cycle c
-    ac = acceptance c
+    ac = acceptances c
 
 mhgOpenMonitors ::
   AnalysisName ->
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
@@ -92,10 +92,9 @@
     -- | The 'Trace' of the Markov chain. In contrast to the link, the trace is
     -- updated only after all proposals in the cycle have been executed.
     trace :: Trace a,
-    -- | For each 'Proposal', store the list of accepted (True) and rejected (False)
-    -- proposals; for reasons of efficiency, the list is also stored in reverse
-    -- order.
-    acceptance :: Acceptance (Proposal a),
+    -- | For each 'Proposal', store actual acceptance counts and for some
+    -- proposals also expected acceptance rates.
+    acceptances :: Acceptances (Proposal a),
     -- | The random number generator.
     generator :: IOGenM StdGen,
     --
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
@@ -51,7 +51,7 @@
   { savedLink :: Link a,
     savedIteration :: Int,
     savedTrace :: C.Stack VB.Vector (Link a),
-    savedAcceptance :: Acceptance Int,
+    savedAcceptances :: Acceptances Int,
     savedSeed :: (Word64, Word64),
     savedTuningParameters :: [Maybe (TuningParameter, AuxiliaryTuningParameters)]
   }
@@ -109,18 +109,18 @@
                 "fromSave: Given likelihood:" <> show (lh $ state it) <> "."
               ]
        in error msg
-  | length (fromAcceptance ac) /= length (ccProposals cc) =
+  | length (fromAcceptances ac) /= length (ccProposals cc) =
       let msg =
             unlines
               [ "fromSave: The number of proposals does not match.",
-                "fromSave: Number of saved proposals:" <> show (length $ fromAcceptance ac) <> ".",
+                "fromSave: Number of saved proposals:" <> show (length $ fromAcceptances ac) <> ".",
                 "fromSave: Number of given proposals:" <> show (length $ ccProposals cc) <> "."
               ]
        in error msg
   | otherwise = fromSavedChainUnsafe pr lh cc mn sv
   where
     it = savedLink sv
-    ac = savedAcceptance sv
+    ac = savedAcceptances sv
 
 -- | See 'fromSavedChain' but do not perform sanity checks. Useful when
 -- restarting a run with changed prior function, likelihood function or
diff --git a/src/Mcmc/Cycle.hs b/src/Mcmc/Cycle.hs
--- a/src/Mcmc/Cycle.hs
+++ b/src/Mcmc/Cycle.hs
@@ -15,7 +15,7 @@
 module Mcmc.Cycle
   ( -- * Cycles
     Order (..),
-    Cycle (ccProposals, ccRequireTrace),
+    Cycle (ccProposals, ccRequireTrace, ccHasIntermediateTuners),
     cycleFromList,
     setOrder,
     IterationMode (..),
@@ -27,10 +27,12 @@
   )
 where
 
+import Control.Applicative
 import qualified Data.ByteString.Builder as BB
 import qualified Data.ByteString.Lazy.Char8 as BL
 import Data.List
 import qualified Data.Map.Strict as M
+import Data.Maybe
 import qualified Data.Vector as VB
 import Mcmc.Acceptance
 import Mcmc.Internal.Shuffle
@@ -91,7 +93,10 @@
   { ccProposals :: [Proposal a],
     ccOrder :: Order,
     -- | Does the cycle require the trace when auto tuning? See 'tRequireTrace'.
-    ccRequireTrace :: Bool
+    ccRequireTrace :: Bool,
+    -- | Does the cycle include proposals that can be tuned every iterations?
+    -- See 'tSuitableForIntermediateTuning'.
+    ccHasIntermediateTuners :: Bool
   }
 
 -- | Create a 'Cycle' from a list of 'Proposal's; use 'RandomO', but see 'setOrder'.
@@ -100,7 +105,7 @@
   error "cycleFromList: Received an empty list but cannot create an empty Cycle."
 cycleFromList xs =
   if length uniqueXs == length xs
-    then Cycle xs RandomO (any needsTrace xs)
+    then Cycle xs RandomO (any needsTrace xs) (any isIntermediate xs)
     else error $ "\n" ++ msg ++ "cycleFromList: Proposals are not unique."
   where
     uniqueXs = nub xs
@@ -109,9 +114,8 @@
     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
+    needsTrace p = maybe False tRequireTrace (prTuner p)
+    isIntermediate p = maybe False tSuitableForIntermediateTuning (prTuner p)
 
 -- | Set the order of 'Proposal's in a 'Cycle'.
 setOrder :: Order -> Cycle a -> Cycle a
@@ -123,9 +127,13 @@
 
 -- | Replicate 'Proposal's according to their weights and possibly shuffle them.
 prepareProposals :: StatefulGen g m => IterationMode -> Cycle a -> g -> m [Proposal a]
-prepareProposals m (Cycle xs o _) g =
+prepareProposals m (Cycle xs o _ _) g =
   if null ps
-    then error "prepareProposals: No proposals found."
+    then
+      let msg = case m of
+            FastProposals -> "no fast proposals found"
+            AllProposals -> "no proposals found"
+       in error $ "prepareProposals: " <> msg
     else case o of
       RandomO -> shuffle ps g
       SequentialO -> return ps
@@ -146,7 +154,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
@@ -160,40 +168,61 @@
 -- See 'tuneWithTuningParameters' and 'Tuner'.
 tuneWithChainParameters ::
   TuningType ->
-  AcceptanceRate ->
+  Maybe AcceptanceRate ->
   Maybe (VB.Vector a) ->
   Proposal a ->
   Either String (Proposal a)
-tuneWithChainParameters b ar mxs p = case prTuner p of
+tuneWithChainParameters tt mar 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
+  Just (Tuner t ts rt it fT _) -> case (tt, it, prSpeed p) of
+    (IntermediateTuningFastProposalsOnly, True, PFast) -> tuneIntermediate
+    (IntermediateTuningAllProposals, True, _) -> tuneIntermediate
+    (NormalTuningFastProposalsOnly, _, PFast) -> tuneNormally
+    (NormalTuningAllProposals, _, _) -> tuneNormally
+    (LastTuningFastProposalsOnly, _, _) -> tuneNormally
+    (LastTuningAllProposals, _, _) -> tuneNormally
+    _ -> Right p
+    where
+      hasTrace = isJust mxs
+      err m = Left $ "tuneWithChainParameters: " <> m
+      tuneIntermediate =
+        if hasTrace
+          then err "intermediate tuning but trace provided"
+          else tune
+      tuneNormally =
+        if rt && not hasTrace
+          then err "trace required"
+          else tune
+      tune =
+        let (t', ts') = fT tt (prDimension p) mar mxs (t, ts)
+         in tuneWithTuningParameters t' ts' p
 
+-- (_, False, Just _) ->
+-- (True, _, Just _) ->
+-- (False, True, Nothing) ->
+-- _ ->
+
 -- | Calculate acceptance rates and auto tunes the 'Proposal's in the 'Cycle'.
 --
 -- Do not change 'Proposal's that are not tuneable.
-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
+autoTuneCycle :: TuningType -> Acceptances (Proposal a) -> Maybe (VB.Vector a) -> Cycle a -> Cycle a
+autoTuneCycle tt a mxs c
+  | isJust mxs && not (ccRequireTrace c) = err "trace provided but not required"
+  | otherwise =
+      if sort (M.keys $ fromAcceptances a) == sort ps
+        then c {ccProposals = map tuneF ps}
+        else err "proposals in map and cycle do not match"
+  where
+    err msg = error $ "autoTuneCycle: " <> msg
+    ps = ccProposals c
+    tuneF p =
+      let (_, _, mar, mtr) = acceptanceRate p a
+          -- Favor the expected rate, if available.
+          mr = mtr <|> mar
+       in either error id $ tuneWithChainParameters tt mr mxs p
 
 -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates.
-summarizeCycle :: IterationMode -> Acceptance (Proposal a) -> Cycle a -> BL.ByteString
+summarizeCycle :: IterationMode -> Acceptances (Proposal a) -> Cycle a -> BL.ByteString
 summarizeCycle m a c =
   BL.intercalate "\n" $
     [ "Summary of proposal(s) in cycle.",
diff --git a/src/Mcmc/Internal/SpecFunctions.hs b/src/Mcmc/Internal/SpecFunctions.hs
--- a/src/Mcmc/Internal/SpecFunctions.hs
+++ b/src/Mcmc/Internal/SpecFunctions.hs
@@ -35,7 +35,7 @@
 -- 'Numeric.SpecFunctions.logGamma'.
 logGammaG :: (Typeable a, RealFloat a) => a -> a
 logGammaG z
-  | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z
+  | typeOf z == typeRep (Proxy :: Proxy Double) = unsafeCoerce logGamma z
   | otherwise = logGammaNonDouble z
 {-# SPECIALIZE logGammaG :: Double -> Double #-}
 
@@ -214,7 +214,7 @@
 -- '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
+  | typeRep (Proxy :: Proxy b) == typeRep (Proxy :: Proxy Double) = unsafeCoerce $ logFactorial n
   | otherwise = logFactorialNonDouble n
 {-# SPECIALIZE logFactorialG :: Int -> Double #-}
 
diff --git a/src/Mcmc/MarginalLikelihood.hs b/src/Mcmc/MarginalLikelihood.hs
--- a/src/Mcmc/MarginalLikelihood.hs
+++ b/src/Mcmc/MarginalLikelihood.hs
@@ -147,12 +147,12 @@
 sampleAtPoint x ss lhf a = do
   a'' <- liftIO $ mcmc ss' a'
   let ch'' = fromMHG a''
-      ac = acceptance ch''
+      ac = acceptances ch''
       mAr = sequence $ acceptanceRates ac
   logDebugB "sampleAtPoint: Summarize cycle."
   logDebugB $ summarizeCycle AllProposals ac $ cycle ch''
   case mAr of
-    Nothing -> logWarnB "Some acceptance rates are unavailable. The tuning period may be too small."
+    Nothing -> logWarnB "Some acceptance rates are unavailable."
     Just ar -> do
       unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1."
       unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9."
diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs
--- a/src/Mcmc/Mcmc.hs
+++ b/src/Mcmc/Mcmc.hs
@@ -25,11 +25,13 @@
 import Control.Monad
 import Control.Monad.IO.Class
 import Control.Monad.Trans.Reader
+import Data.Functor
+import Mcmc.Acceptance (ResetAcceptance (ResetEverything, ResetExpectedRatesOnly))
 import Mcmc.Algorithm
 import Mcmc.Cycle
 import Mcmc.Environment
 import Mcmc.Logger
-import Mcmc.Proposal (TuningType (LastTuningStep, NormalTuningStep))
+import Mcmc.Proposal
 import Mcmc.Settings
 import System.IO
 import Prelude hiding (cycle)
@@ -52,7 +54,7 @@
 mcmcResetAcceptance :: Algorithm a => a -> MCMC a
 mcmcResetAcceptance a = do
   logDebugB "Reset acceptance rates."
-  return $ aResetAcceptance a
+  return $ aResetAcceptance ResetEverything a
 
 mcmcExceptionHandler :: Algorithm a => Environment Settings -> a -> AsyncException -> IO b
 mcmcExceptionHandler e a err = do
@@ -77,8 +79,16 @@
   mStdLog <- liftIO $ aExecuteMonitors vb t0 iTotal a
   forM_ mStdLog (logOutB "   ")
 
-mcmcIterate :: Algorithm a => IterationMode -> Int -> a -> MCMC a
-mcmcIterate m n a
+-- When intermediate tuning is activated, specific proposals get tuned every
+-- iterations.
+data IntermediateTuningSpec
+  = IntermediateTuningFastProposalsOnlyOn
+  | IntermediateTuningAllProposalsOn
+  | IntermediateTuningOff
+  deriving (Eq)
+
+mcmcIterate :: Algorithm a => IntermediateTuningSpec -> IterationMode -> Int -> a -> MCMC a
+mcmcIterate t m n a
   | n < 0 = error "mcmcIterate: Number of iterations is negative."
   | n == 0 = return a
   | otherwise = do
@@ -87,7 +97,20 @@
       -- NOTE: Handle interrupts during iterations, before writing monitors,
       -- using the old algorithm state @a@.
       let handlerOld = mcmcExceptionHandler e a
-          actionIterate = aIterate m p a
+          maybeIntermediateAutoTune x =
+            -- Do not perform intermediate tuning at the last step, because a
+            -- normal tuning will be performed.
+            case t of
+              IntermediateTuningFastProposalsOnlyOn
+                | n > 1 ->
+                    aAutoTune IntermediateTuningFastProposalsOnly 1 x
+                      <&> aResetAcceptance ResetExpectedRatesOnly
+              IntermediateTuningAllProposalsOn
+                | n > 1 ->
+                    aAutoTune IntermediateTuningAllProposals 1 x
+                      <&> aResetAcceptance ResetExpectedRatesOnly
+              _ -> pure x
+          actionIterate = aIterate m p a >>= maybeIntermediateAutoTune
       a' <- liftIO $ actionIterate `catch` handlerOld
       -- NOTE: Mask asynchronous exceptions while writing monitor files. Handle
       -- interrupts after writing monitors; use the new state @a'@.
@@ -100,8 +123,8 @@
       -- recover from partly written monitor files.
       let handlerNew = mcmcExceptionHandler e a'
           actionWrite = runReaderT (mcmcExecuteMonitors a') e
-      liftIO $ (uninterruptibleMask_ actionWrite) `catch` handlerNew
-      mcmcIterate m (n - 1) a'
+      liftIO $ uninterruptibleMask_ actionWrite `catch` handlerNew
+      mcmcIterate t m (n - 1) a'
 
 mcmcNewRun :: Algorithm a => a -> MCMC a
 mcmcNewRun a = do
@@ -112,15 +135,15 @@
   mcmcExecuteMonitors a
   when (aIsInvalidState a) (logWarnB "The initial state is invalid!")
   a' <- mcmcBurnIn a
-  logInfoS $ "Cleaning chain after burn in."
+  logInfoS "Cleaning chain after burn in."
   let tl = sTraceLength s
   a'' <- liftIO $ aCleanAfterBurnIn tl a'
-  logInfoS $ "Saving chain after burn in."
+  logInfoS "Saving chain after burn in."
   mcmcSave a''
   let i = fromIterations $ sIterations s
   logInfoS $ "Running chain for " ++ show i ++ " iterations."
   logInfoB $ aStdMonitorHeader a''
-  mcmcIterate AllProposals i a''
+  mcmcIterate IntermediateTuningOff AllProposals i a''
 
 mcmcContinueRun :: Algorithm a => a -> MCMC a
 mcmcContinueRun a = do
@@ -139,7 +162,7 @@
   logInfoB $ aSummarizeCycle AllProposals a
   logInfoS $ "Running chain for " ++ show di ++ " iterations."
   logInfoB $ aStdMonitorHeader a
-  mcmcIterate AllProposals di a
+  mcmcIterate IntermediateTuningOff AllProposals di a
 
 mcmcBurnIn :: Algorithm a => a -> MCMC a
 mcmcBurnIn a = do
@@ -154,7 +177,7 @@
       logInfoS $ "Burning in for " <> show n <> " iterations."
       logInfoS "Auto tuning is disabled."
       logInfoB $ aStdMonitorHeader a
-      a' <- mcmcIterate AllProposals n a
+      a' <- mcmcIterate IntermediateTuningOff AllProposals n a
       logInfoB $ aSummarizeCycle AllProposals a'
       a'' <- mcmcResetAcceptance a'
       logInfoB "Burn in finished."
@@ -191,25 +214,34 @@
 
 -- Auto tune the proposals.
 mcmcAutotune :: Algorithm a => TuningType -> Int -> a -> MCMC a
-mcmcAutotune NormalTuningStep n a = do
-  logDebugB "Intermediate auto tune."
-  liftIO $ aAutoTune NormalTuningStep n a
-mcmcAutotune LastTuningStep n a = do
-  logDebugB "Last auto tune."
-  liftIO $ aAutoTune LastTuningStep n a
+mcmcAutotune t n a = do
+  case t of
+    NormalTuningFastProposalsOnly -> logDebugB "Normal auto tune; fast proposals only."
+    IntermediateTuningFastProposalsOnly -> pure ()
+    LastTuningFastProposalsOnly -> logDebugB "Last auto tune; fast proposals only."
+    NormalTuningAllProposals -> logDebugB "Normal auto tune; all proposals."
+    IntermediateTuningAllProposals -> pure ()
+    LastTuningAllProposals -> logDebugB "Last auto tune; all proposals."
+  liftIO $ aAutoTune t 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 LastTuningStep x a'
+  let (tti, ttl) = case m of
+        FastProposals -> (IntermediateTuningFastProposalsOnlyOn, LastTuningFastProposalsOnly)
+        AllProposals -> (IntermediateTuningAllProposalsOn, LastTuningAllProposals)
+  a' <- mcmcIterate tti m x a
+  a'' <- mcmcAutotune ttl 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 NormalTuningStep x a'
+  let (tti, ttn) = case m of
+        FastProposals -> (IntermediateTuningFastProposalsOnlyOn, NormalTuningFastProposalsOnly)
+        AllProposals -> (IntermediateTuningAllProposalsOn, NormalTuningAllProposals)
+  a' <- mcmcIterate tti m x a
+  a'' <- mcmcAutotune ttn x a'
   logDebugB $ aSummarizeCycle m a''
   logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations."
   logDebugB $ aStdMonitorHeader a''
diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs
--- a/src/Mcmc/Proposal.hs
+++ b/src/Mcmc/Proposal.hs
@@ -38,8 +38,6 @@
     TuningFunction,
     AuxiliaryTuningParameters,
     tuningFunction,
-    tuningFunctionWithAux,
-    tuningFunctionOnlyAux,
     tuningParameterMin,
     tuningParameterMax,
     tuneWithTuningParameters,
@@ -207,16 +205,14 @@
     Propose !a !KernelRatio !Jacobian
   deriving (Show, Eq)
 
--- 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.
 --
--- Maybe report acceptance counts internal to the proposal (e.g., used by
+-- Maybe report acceptance rates internal to the proposal (e.g., used by
 -- proposals based on Hamiltonian dynamics).
-type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceCounts)
+type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceRates)
 
 -- | Create a proposal with a single tuning parameter.
 --
@@ -243,7 +239,7 @@
   where
     fT = tuningFunction
     g t _ = Right $ f t
-    tuner = Tuner 1.0 VU.empty False fT g
+    tuner = Tuner 1.0 VU.empty False False fT g
 createProposal r f s d n w NoTune =
   Proposal n r s d w (f 1.0) Nothing
 
@@ -253,6 +249,8 @@
     tAuxiliaryTuningParameters :: AuxiliaryTuningParameters,
     -- | Does the tuner require the trace over the last tuning period?
     tRequireTrace :: Bool,
+    -- | Can the tuner be used for intermediate tuning (see 'TuningType')?
+    tSuitableForIntermediateTuning :: Bool,
     tTuningFunction :: TuningFunction a,
     -- | Given the tuning parameter, and the auxiliary tuning parameters, get
     -- the tuned propose function.
@@ -275,16 +273,35 @@
 -- expected acceptance rate; and vice versa.
 type TuningParameter = Double
 
--- | The last tuning step may be special.
-data TuningType = NormalTuningStep | LastTuningStep
+-- | Tuning type. To distinguish between fast and slow proposals, see
+-- 'Mcmc.Cycle.IterationMode'.
+data TuningType
+  = -- | Normal tuning step with fast proposals only.
+    NormalTuningFastProposalsOnly
+  | -- | Intermediate tuning step executed after each iteration with fast
+    -- proposals only. Only suitable for proposals which can calculate expected
+    -- acceptance rates.
+    IntermediateTuningFastProposalsOnly
+  | -- | The last tuning step with fast proposals only may be special.
+    LastTuningFastProposalsOnly
+  | -- | Normal tuning step of all proposals.
+    NormalTuningAllProposals
+  | -- | Intermediate tuning step of all proposals.
+    IntermediateTuningAllProposals
+  | -- | The last tuning step with all proposals.
+    LastTuningAllProposals
+  deriving (Eq)
 
 -- | 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.
+  -- | Acceptance rate of last tuning period. May not always be available
+  -- because proposals may be skipped.
+  Maybe AcceptanceRate ->
+  -- | Trace of last tuning period. Not available for intermediate tuning' steps
+  -- (see 'TuningType'), and only available for other tuning types when
+  -- requested by proposal.
   Maybe (VB.Vector a) ->
   (TuningParameter, AuxiliaryTuningParameters) ->
   (TuningParameter, AuxiliaryTuningParameters)
@@ -301,27 +318,15 @@
 
 -- | Default tuning function.
 --
--- 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.
+-- The default tuning function only uses the actual acceptance rate. In
+-- particular, it does not handle auxiliary tuning parameters, ignores
+-- intermediate tuning steps, and ignores the actual samples attained during the
+-- last tuning period.
 tuningFunction :: TuningFunction a
-tuningFunction _ d r _ (!t, !ts) = bimap (tuningFunctionSimple d r) id (t, ts)
-
--- | 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) (!t, !ts) = bimap (tuningFunctionSimple d r) (f b xs) (t, ts)
-
--- | 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) (!t, !ts) = bimap id (f b xs) (t, ts)
+tuningFunction IntermediateTuningFastProposalsOnly _ _ _ t = t
+tuningFunction IntermediateTuningAllProposals _ _ _ t = t
+tuningFunction _ _ Nothing _ t = t
+tuningFunction _ d (Just r) _ (!t, !ts) = first (tuningFunctionSimple d r) (t, ts)
 
 -- 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
@@ -330,7 +335,7 @@
 
 -- | Minimal tuning parameter; subject to change.
 tuningParameterMin :: TuningParameter
-tuningParameterMin = 1e-5
+tuningParameterMin = 1e-6
 
 -- | Maximal tuning parameter; subject to change.
 tuningParameterMax :: TuningParameter
@@ -358,14 +363,14 @@
   Either String (Proposal a)
 tuneWithTuningParameters t ts p = case prTuner p of
   Nothing -> Left "tuneWithTuningParameters: Proposal is not tunable."
-  Just (Tuner _ _ nt fT g) ->
+  Just (Tuner _ _ reqTr inTn 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 {prFunction = ps, prTuner = Just $ Tuner t'' ts nt fT g}
+          Right ps -> Right $ p {prFunction = ps, prTuner = Just $ Tuner t'' ts reqTr inTn fT g}
 
 -- | See 'PDimension'.
 getOptimalRate :: PDimension -> Double
@@ -431,7 +436,7 @@
 
 -- Lift tuner from one data type to another.
 liftTunerWith :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b
-liftTunerWith jf l (Tuner p ps nt fP g) = Tuner p ps nt fP' g'
+liftTunerWith jf l (Tuner p ps reqTr inTn fP g) = Tuner p ps reqTr inTn fP' g'
   where
     fP' b d r = fP b d r . fmap (VB.map (^. l))
     g' x xs = liftPFunctionWith jf l <$> g x xs
@@ -455,14 +460,15 @@
   BL.ByteString ->
   BL.ByteString ->
   BL.ByteString
-renderRow name ptype weight nAccept nReject acceptRate optimalRate tuneParam manualAdjustment = nm <> pt <> wt <> na <> nr <> ra <> ro <> tp <> mt
+renderRow name ptype weight nAccept nReject acceptRateActual optimalRate tuneParam manualAdjustment =
+  nm <> pt <> wt <> na <> nr <> ra <> ro <> tp <> mt
   where
     nm = alignLeft 30 name
     pt = alignLeft 50 ptype
     wt = alignRight 8 weight
     na = alignRight 14 nAccept
     nr = alignRight 14 nReject
-    ra = alignRight 14 acceptRate
+    ra = alignRight 14 acceptRateActual
     ro = alignRight 14 optimalRate
     tp = alignRight 20 tuneParam
     mt = alignRight 30 manualAdjustment
@@ -476,7 +482,7 @@
     "Weight"
     "Accepted"
     "Rejected"
-    "Rate"
+    "Actual rate"
     "Optimal rate"
     "Tuning parameter"
     "Consider manual adjustment"
@@ -488,7 +494,7 @@
   PWeight ->
   Maybe TuningParameter ->
   PDimension ->
-  Maybe (Int, Int, Double) ->
+  (Int, Int, Maybe Double, Maybe Double) ->
   BL.ByteString
 summarizeProposal name description weight tuningParameter dimension ar =
   renderRow
@@ -497,18 +503,19 @@
     weightStr
     nAccept
     nReject
-    acceptRate
+    acceptRateActual
+    -- acceptRateExpected
     optimalRate
     tuneParamStr
     manualAdjustmentStr
   where
     fN n = BB.formatDouble (BB.standard n)
     weightStr = BB.toLazyByteString $ BB.intDec $ fromPWeight weight
-    nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) ar
-    nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) ar
-    acceptRate = BB.toLazyByteString $ maybe "" (fN 2 . (^. _3)) ar
+    nAccept = BB.toLazyByteString $ BB.intDec $ ar ^. _1
+    nReject = BB.toLazyByteString $ BB.intDec $ ar ^. _2
+    acceptRateActual = BB.toLazyByteString $ maybe "" (fN 2) (ar ^. _3)
     optimalRate = BB.toLazyByteString $ fN 2 $ getOptimalRate dimension
-    tuneParamStr = BB.toLazyByteString $ maybe "" (fN 4) tuningParameter
+    tuneParamStr = BB.toLazyByteString $ maybe "" (fN 6) tuningParameter
     checkRate rate
       | rate < rateMin = Just "rate too low"
       | rate > rateMax = Just "rate too high"
@@ -518,7 +525,8 @@
       | tp >= (0.9 * tuningParameterMax) = Just "tuning parameter too high"
       | otherwise = Nothing
     tps = checkTuningParam =<< tuningParameter
-    ars = (checkRate . (^. _3)) =<< ar
+    -- Use actual acceptance rate.
+    ars = checkRate =<< (ar ^. _3)
     manualAdjustmentStr =
       let
        in case (ars, tps) of
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
@@ -19,26 +19,25 @@
 import Numeric.Log
 import Statistics.Distribution
 
--- | Generic function to create proposals for continuous parameters (e.g.,
--- 'Double').
+-- | Generic function to create proposals using a continuous auxiliary variable
+-- of type 'Double'.
 --
 -- The procedure is as follows: Let \(\mathbb{X}\) be the state space and \(x\)
 -- be the current state.
 --
 -- 1. Let \(D\) be a continuous probability distribution on \(\mathbb{D}\);
---    sample an auxiliary variable \(epsilon \sim D\).
+--    sample an auxiliary variable \(u \sim D\).
 --
--- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). PFunction a
---    new state \(x' = x \odot \epsilon\).
+-- 2. Let \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). Propose a
+--    new state \(x' = x \odot u\).
 --
--- 3. If the proposal is unbiased, the Metropolis-Hastings-Green ratio can
---    directly be calculated using the posterior function.
+-- If the proposal is unbiased, the Metropolis-Hastings-Green ratio can directly
+-- be calculated using the posterior function.
 --
--- 4. However, if the proposal is biased: Suppose \(g : \mathbb{D} \to
---    \mathbb{D}\) inverses the auxiliary variable \(\epsilon\) such that \(x =
---    x' \odot g(\epsilon)\). Calculate the Metropolis-Hastings-Green ratio
---    using the posterior function, \(g\), \(D\), \(\epsilon\), and possibly a
---    Jacobian function.
+-- However, if the proposal is biased: Let \(g : \mathbb{D} \to \mathbb{D}\);
+-- \(g\) inverses the auxiliary variable \(u\) such that \(x = x' \odot g(u)\).
+-- Calculate the Metropolis-Hastings-Green ratio using the posterior function,
+-- \(g\), \(D\), \(u\), and possibly a Jacobian function.
 genericContinuous ::
   (ContDistr d, ContGen d) =>
   -- | Probability distribution
@@ -46,12 +45,12 @@
   -- | Forward operator \(\odot\).
   --
   -- For example, for a multiplicative proposal on one variable the forward
-  -- operator is @(*)@, so that @x * u = y@.
+  -- operator is @(*)@, so that \(x' = x * u\).
   (a -> Double -> a) ->
   -- | Inverse operator \(g\) of the auxiliary variable.
   --
   -- For example, 'recip' for a multiplicative proposal on one variable, since
-  -- @y * (recip u) = x * u * (recip u) = x@.
+  -- \(x' * u^{-1} = x * u * u^{-1} = x\).
   --
   -- Required for biased proposals.
   Maybe (Double -> Double) ->
@@ -86,7 +85,8 @@
   pure (Propose (x `f` u) r j, Nothing)
 {-# INLINEABLE genericContinuous #-}
 
--- | Generic function to create proposals for discrete parameters (e.g., 'Int').
+-- | Generic function to create proposals using a discrete auxiliary variable of
+-- type 'Int'.
 --
 -- See 'genericContinuous'.
 genericDiscrete ::
@@ -95,11 +95,11 @@
   d ->
   -- | Forward operator.
   --
-  -- For example, (+), so that x + dx = x'.
+  -- For example, (+), so that \(x + dx = x'\).
   (a -> Int -> a) ->
   -- | Inverse operator \(g\) of the auxiliary variable.
   --
-  -- For example, 'negate', so that x' + (negate dx) = x.
+  -- For example, 'negate', so that \(x' - dx = x + dx - dx = x\).
   --
   -- Only required for biased proposals.
   Maybe (Int -> Int) ->
diff --git a/src/Mcmc/Proposal/Hamiltonian/Common.hs b/src/Mcmc/Proposal/Hamiltonian/Common.hs
--- a/src/Mcmc/Proposal/Hamiltonian/Common.hs
+++ b/src/Mcmc/Proposal/Hamiltonian/Common.hs
@@ -59,7 +59,8 @@
 import qualified Numeric.LinearAlgebra as L
 
 -- NOTE: Implementing the Riemannian adaptation (state-dependent mass matrix).
--- seems a little bit of an overkill.
+-- seems a little bit of an overkill. See also
+-- https://discourse.mc-stan.org/t/riemann-manifold-hmc-in-stan/19466/5.
 
 -- | The Hamiltonian proposal acts on a vector of floating point values referred
 -- to as positions.
diff --git a/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs b/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs
--- a/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs
+++ b/src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs
@@ -108,14 +108,6 @@
   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 ::
@@ -130,11 +122,13 @@
   -- 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]
+      lL = max (1 :: Int) (floor $ 0.9 * lM)
+      lR = max 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)
+    -- NOTE: I am not sure if it is correct to set the expected acceptance rate
+    -- to 0 when the leapfrog integrator fails.
+    Nothing -> pure (ForceReject, Just $ AcceptanceRates 0 1)
     -- 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
@@ -149,13 +143,9 @@
       -- 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)
+          -- Limit expected acceptance rate between 0 and 1.
+          ar = max 0 $ min 1 (exp $ ln r)
+      pure (pr, Just $ AcceptanceRates ar 1)
   where
     (HParamsI e la ms _ _ msI mus) = hparamsi
     (HStructure _ toVec fromVec) = hstruct
@@ -204,7 +194,7 @@
       tuner = do
         tfun <- hTuningFunctionWith dim toVec htconf
         let pfun = hamiltonianPFunctionWithTuningParameters dim hstruct targetWith
-        pure $ Tuner 1.0 ts True tfun pfun
+        pure $ Tuner 1.0 ts True 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
--- a/src/Mcmc/Proposal/Hamiltonian/Internal.hs
+++ b/src/Mcmc/Proposal/Hamiltonian/Internal.hs
@@ -55,6 +55,7 @@
 import Control.Monad
 import Control.Monad.ST
 import Data.Foldable
+import Data.Maybe
 import qualified Data.Vector.Storable as VS
 import qualified Data.Vector.Unboxed as VU
 import Mcmc.Proposal
@@ -92,31 +93,48 @@
   }
   deriving (Show)
 
--- The default tuning parameters in [4] are:
+-- The default tuning parameters in [4] which have been tweaked for tuning the
+-- proposals after every iteration 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.
+-- For reference, I used the following default parameters with longer auto
+-- tuning intervals.
 --
--- 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.
+--   mu = log $ 10 * eps
+--   ga = 0.1
+--   t0 = 3
+--   ka = 0.5
 --
--- NOTE: In theory, these we could expose these internal tuning parameters to
--- the user.
+-- Another good resource:
+-- https://mc-stan.org/docs/2_29/reference-manual/hmc-algorithm-parameters.html.
+--
+-- NOTE: In theory, we could expose these internal tuning parameters to the
+-- user.
 tParamsFixedWith :: LeapfrogScalingFactor -> TParamsFixed
 tParamsFixedWith eps = TParamsFixed eps mu ga t0 ka
   where
+    -- "Mu is a freely chosen point that the iterators are shrunk towards". I am
+    -- not exactly sure what this means. The parameter does not seem to have
+    -- much of an effect.
     mu = log $ 10 * eps
-    ga = 0.1
-    t0 = 3
-    ka = 0.5
+    -- Gamma "controls the amount of shrinkage towards mu". The larger gamma is,
+    -- the less variant epsilon is.
+    --
+    -- I changed this parameter from 0.05 to get better results in test runs.
+    ga = 0.15
+    -- "Free parameter that stabilizes the initial iterations". The larger t0
+    -- is, the stabler epsilon is in the first iterations.
+    t0 = 10
+    -- "Setting the parameter ka < 1 allows us to give higher weight to more
+    -- recent iterates and to more quickly forget the iterates produced during
+    -- the early warmup stages."
+    --
+    -- I changed this parameter from 0.75 to get better results in test runs.
+    ka = 0.75
 
 -- All internal parameters.
 data HParamsI = HParamsI
@@ -235,7 +253,7 @@
           a :: Double
           a = if rI > 0.5 then 1 else (-1)
           go e r =
-            if r ** a > 2 ** (negate a)
+            if r ** a > 2 ** negate a
               then case leapfrog t msI 1 e q p of
                 Nothing -> e
                 Just (_, p'', _, prQ'') ->
@@ -257,38 +275,55 @@
   (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)
+hTuningFunctionWith _ _ (HTuningConf HNoTuneLeapfrog HNoTuneMasses) = Nothing
+hTuningFunctionWith n toVec (HTuningConf lc mc) = Just $ \tt pdim mar mxs (_, !ts) ->
+  case tt of
+    IntermediateTuningFastProposalsOnly -> err "fast intermediate tuning step but slow proposal"
+    NormalTuningFastProposalsOnly -> err "fast normal tuning step but slow proposal"
+    _ ->
+      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
+          m' = SmoothingParameter $ round m
+          (ms', msI') = case tt of
+            IntermediateTuningAllProposals -> (ms, msI)
+            _ ->
+              let xs = fromMaybe (err "empty trace") mxs
+               in case mc of
+                    HNoTuneMasses -> (ms, msI)
+                    HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly m' toVec xs (ms, msI)
+                    HTuneAllMasses -> tuneAllMasses m' toVec xs (ms, msI)
+          (eps'', epsMean'', h'') = case tt of
+            LastTuningFastProposalsOnly -> (eps, epsMean, h)
+            _ -> case lc of
+              HNoTuneLeapfrog -> (eps, epsMean, h)
+              HTuneLeapfrog ->
+                let ar = fromMaybe (err "no acceptance rate") mar
+                    delta = getOptimalRate pdim
+                    -- Algorithm 6; explained in Section 3.2.
+                    --
+                    -- Another good resource is the Tensorflow API
+                    -- documentation:
+                    -- https://www.tensorflow.org/probability/api_docs/python/tfp/mcmc/DualAveragingStepSizeAdaptation.
+                    --
+                    -- See also Nesterov (2007) Primal-dual subgradient methods
+                    -- for convex problems, Mathematical Programming.
+                    c = recip $ m + t0
+                    h' = (1.0 - c) * h + c * (delta - ar)
+                    eps' = exp $ mu - (sqrt m / ga) * h'
+                    mMKa = m ** negate ka
+                    -- Original formula is:
+                    -- epsMean' = exp $ mMKa * logEps' + (1 - mMKa) * log epsMean
+                    -- Which is the same as:
+                    epsMean' = (eps' ** mMKa) * (epsMean ** (1 - mMKa))
+                    epsF = if tt == LastTuningAllProposals then epsMean' else eps'
+                 in (epsF, epsMean', h')
+          tpv' = TParamsVar epsMean'' h'' (m + 1.0)
+       in (eps'' / eps0, toAuxiliaryTuningParameters $ HParamsI eps'' la ms' tpv' tpf msI' mus)
+  where
+    err msg = error $ "hTuningFunctionWith: " <> msg
 
 checkHStructureWith :: Foldable s => Masses -> HStructure s -> Maybe String
 checkHStructureWith ms (HStructure x toVec fromVec)
@@ -361,7 +396,7 @@
           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)
+  (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.
@@ -379,7 +414,7 @@
           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')
+                then go (n - 1) $ Just (qs', ps')
                 else Nothing
 
 leapfrogStepMomenta ::
diff --git a/src/Mcmc/Proposal/Hamiltonian/Masses.hs b/src/Mcmc/Proposal/Hamiltonian/Masses.hs
--- a/src/Mcmc/Proposal/Hamiltonian/Masses.hs
+++ b/src/Mcmc/Proposal/Hamiltonian/Masses.hs
@@ -21,6 +21,7 @@
     massesToVector,
 
     -- * Tuning
+    SmoothingParameter (..),
     tuneDiagonalMassesOnly,
     tuneAllMasses,
   )
@@ -32,6 +33,7 @@
 import qualified Data.Vector.Unboxed as VU
 import Mcmc.Proposal.Hamiltonian.Common
 import qualified Numeric.LinearAlgebra as L
+import Numeric.Natural
 import qualified Statistics.Covariance as S
 import qualified Statistics.Function as S
 import qualified Statistics.Sample as S
@@ -108,7 +110,7 @@
 -- 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
+  L.diag (L.cmap cleanDiag xsDiag) + L.cmap cleanOffDiag xsOffDiag
   where
     xsDiag = L.takeDiag xs
     cleanDiag x
@@ -129,9 +131,9 @@
   | 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
+  | otherwise = toGMatrix $ cleanMatrix xsI
   where
-    xs' = L.unSym $ xs
+    xs' = L.unSym xs
     n = L.rows xs'
     m = L.cols xs'
     (xsI, (_, sign)) = L.invlndet xs'
@@ -179,20 +181,68 @@
 getSampleSize :: VS.Vector Double -> Int
 getSampleSize = VS.length . VS.uniq . S.gsort
 
-getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double
-getNewMassDiagonalWithRescue sampleSize massOld massEstimate
+rescueAfter :: Double -> Double
+rescueAfter y
+  | massMin > abs y = signum y * massMin
+  | abs y > massMax = signum y * massMax
+  | otherwise = y
+
+rescueBefore :: (Double -> Double -> Double) -> Double -> Double -> Double
+rescueBefore f old new
+  -- Be permissive with NaN and infinite values.
+  | isNaN new = old
+  | isInfinite new = old
+  | otherwise = f old new
+
+rescue :: (Double -> Double -> Double) -> Double -> Double -> Double
+rescue f old new = rescueAfter $ rescueBefore f old new
+
+-- -- I do not know where I got this function from, but it works pretty well!
+-- interpolate :: Double -> Double -> Double
+-- interpolate old new = interSqrt ** 2
+--   where
+--     interSqrt = recip 3 * (sqrt old + 2 * sqrt new)
+
+-- | This parameter plays the same role as @m@ in the dual averaging algorithm.
+-- In the beginning, when @m@ is zero, the newly estimated masses will take full
+-- precedence over the current masses. After some auto tuning steps, when @m@ is
+-- larger, the newly estimated masses influence the current masses only
+-- slightly.
+newtype SmoothingParameter = SmoothingParameter Natural
+
+-- Another interpolation function I came up with. It is pretty cool, because
+-- (similar to above) it interpolates the square roots, which is what we want.
+interpolate' :: SmoothingParameter -> Double -> Double -> Double
+interpolate' (SmoothingParameter m) oldSquared newSquared = finSign * (fin ** 2)
+  where
+    oldSign = signum oldSquared
+    newSign = signum newSquared
+    sqrt' = sqrt . abs
+    old = oldSign * sqrt' oldSquared
+    new = newSign * sqrt' newSquared
+    -- The new mass will be the second last boundary. That is, the larger the
+    -- number of bins is, the more informative the new mass is compared to the
+    -- old mass.
+    nbins = max (100 - fromIntegral m) 2
+    fin = recip nbins * (old + (nbins - 1) * new)
+    finSign = signum fin
+
+getNewMassDiagonalWithSampleSize :: SmoothingParameter -> Int -> Double -> Double -> Double
+getNewMassDiagonalWithSampleSize m 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
+  -- Diagonal masses are variances which are strictly positive.
   | massEstimate <= 0 = massOld
-  | massMin > massNew = massMin
-  | massNew > massMax = massMax
-  | otherwise = massNew
-  where
-    massNewSqrt = recip 3 * (sqrt massOld + 2 * sqrt massEstimate)
-    massNew = massNewSqrt ** 2
+  | otherwise = rescue (interpolate' m) massOld massEstimate
 
+getNewMassDiagonal :: SmoothingParameter -> Double -> Double -> Double
+getNewMassDiagonal m massOld massEstimate
+  -- Diagonal masses are variances which are strictly positive.
+  | massEstimate <= 0 = massOld
+  | otherwise = rescue (interpolate' m) massOld massEstimate
+
+getNewMassOffDiagonal :: SmoothingParameter -> Double -> Double -> Double
+getNewMassOffDiagonal m = rescue (interpolate' m)
+
 -- 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
@@ -214,7 +264,7 @@
     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)
+    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
@@ -231,6 +281,7 @@
            in go x' (k + 1)
 
 tuneDiagonalMassesOnly ::
+  SmoothingParameter ->
   -- Conversion from value to vector.
   (a -> Positions) ->
   -- Value vector.
@@ -242,7 +293,7 @@
 -- 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)
+tuneDiagonalMassesOnly m 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."
@@ -271,16 +322,40 @@
     msDiagonalEstimate = VS.fromList $ map (recip . S.variance) $ L.toColumns xs''
     msDiagonalNew =
       VS.zipWith3
-        getNewMassDiagonalWithRescue
+        (getNewMassDiagonalWithSampleSize m)
         sampleSizes
         msDiagonalOld
         msDiagonalEstimate
 
 -- This value was carefully tuned using the example "hamiltonian".
 defaultGraphicalLassoPenalty :: Double
-defaultGraphicalLassoPenalty = 0.4
+defaultGraphicalLassoPenalty = 0.3
 
+interpolateElementWise :: SmoothingParameter -> L.Matrix Double -> L.Matrix Double -> L.Matrix Double
+interpolateElementWise m old new
+  | mO /= mN = err "different number of rows"
+  | nO /= nN = err "different number of columns"
+  | mO /= nO = err "not square"
+  | mO < 1 = err "empty matrix"
+  | otherwise = L.build (mO, nO) f
+  where
+    mO = L.rows old
+    nO = L.cols old
+    mN = L.rows new
+    nN = L.cols new
+    err msg = error $ "interpolateElementWise: " <> msg
+    f iD jD =
+      let -- This sucks a bit, because we need a function (e -> e -> e), and
+          -- since the return type is Double, the indices are also Doubble.
+          i = round iD
+          j = round jD
+          g = if i == j then getNewMassDiagonal m else getNewMassOffDiagonal m
+          eO = old `L.atIndex` (i, j)
+          eN = new `L.atIndex` (i, j)
+       in g eO eN
+
 tuneAllMasses ::
+  SmoothingParameter ->
   -- Conversion from value to vector.
   (a -> Positions) ->
   -- Value vector.
@@ -292,16 +367,16 @@
 -- 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)
+tuneAllMasses m 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')
+  | otherwise = (msNew, msINew)
   where
-    fallbackDiagonal = tuneDiagonalMassesOnly toVec xs (ms, msI)
+    fallbackDiagonal = tuneDiagonalMassesOnly m 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.
@@ -314,16 +389,18 @@
     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
+    -- The first value is the covariance matrix sigma, which the inverted mass
+    -- matrix (precision matrix). However, we interpolate the new mass matrix
+    -- using the old one and the new estimate, so we have to recalculate the
+    -- covariance matrix anyways.
+    (_, precNormalized) =
+      either error id $
+        S.graphicalLasso defaultGraphicalLassoPenalty xsNormalized
+    ms' = S.rescalePWith ss (L.unSym precNormalized)
+    -- Clean NaNs, infinities; ensure positive definiteness. The masses should
+    -- be positive definite, but sometimes they happen to be not because of
+    -- numerical errors.
+    msNewDirty = interpolateElementWise m (L.unSym ms) ms'
     -- 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'
+    msNew = L.trustSym $ findClosestPositiveDefiniteMatrix $ cleanMatrix msNewDirty
+    msINew = getMassesI msNew
diff --git a/src/Mcmc/Proposal/Hamiltonian/Nuts.hs b/src/Mcmc/Proposal/Hamiltonian/Nuts.hs
--- a/src/Mcmc/Proposal/Hamiltonian/Nuts.hs
+++ b/src/Mcmc/Proposal/Hamiltonian/Nuts.hs
@@ -79,7 +79,7 @@
 -- 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'.
+  -- calculate the expected acceptance rate 'Alpha'.
   Log Double ->
   MassesI ->
   Target ->
@@ -108,8 +108,8 @@
                 n' = if u <= expEPot' * expEKin' then 1 else 0
                 errorIsSmall = u < deltaMax * expETot'
                 alpha' = expETot' / expETot0
-                alpha = min 1.0 alpha'
-
+                -- Limit expected acceptance rate between 0 and 1.
+                alpha = max 0 $ min 1 alpha'
   -- Recursive case. This is complicated because the algorithm is written for an
   -- imperative language, and because we have two stacked monads.
   | otherwise = do
@@ -138,10 +138,10 @@
             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'
+              let n'''' = n' + n'''
+                  q'''' = if b < fromIntegral n''' / fromIntegral 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).
@@ -152,7 +152,7 @@
   where
     buildTree = buildTreeWith expETot0 msI tfun g
 
--- | Paramters of the NUTS proposal.
+-- | Parameters of the NUTS proposal.
 --
 -- Includes tuning parameters and tuning configuration.
 data NParams = NParams
@@ -184,8 +184,8 @@
 
 data IsNew
   = Old
-  | OldWith {_acceptanceCountsOld :: AcceptanceCounts}
-  | NewWith {_acceptanceCountsNew :: AcceptanceCounts}
+  | OldWith {_acceptanceRatesOld :: AcceptanceRates}
+  | NewWith {_acceptanceRatesNew :: AcceptanceRates}
 
 -- First function in Algorithm 3.
 nutsPFunction ::
@@ -230,28 +230,28 @@
           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."
+                -- Individual expected acceptance rates are limited in
+                -- 'buildTreeWith'.
+                ar = max 0 $ exp $ ln a
+                ars = AcceptanceRates ar na
             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)
+            let (y''', isNew') = if isAccept then (y'', NewWith ars) else (y, OldWith ars)
                 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)
+    -- NOTE: I am not sure if it is correct to set the expected acceptance rate
+    -- to 0 when the no u-turn leapfrog integrator fails.
+    Old -> (ForceReject, Just $ AcceptanceRates 0 1)
+    OldWith ac -> (ForceReject, Just ac)
+    NewWith ac -> (ForceAccept $ fromVec x', Just ac)
   where
     (HParamsI e _ ms _ _ msI mus) = hparamsi
     (HStructure _ toVec fromVecWith) = hstruct
@@ -299,7 +299,7 @@
       tuner = do
         tfun <- hTuningFunctionWith dim toVec htconf
         let pfun = nutsPFunctionWithTuningParameters dim hstruct targetWith
-        pure $ Tuner 1.0 ts True tfun pfun
+        pure $ Tuner 1.0 ts True True tfun pfun
    in case checkHStructureWith (hpsMasses hParamsI) hstruct of
         Just err -> error err
         Nothing -> nutsWith tuner