mcmc 0.6.1.0 → 0.6.2.0

raw patch · 17 files changed

+192/−186 lines, 17 files

Files

• ChangeLog.md +8/−1
• mcmc.cabal +2/−3
• src/Mcmc/Algorithm/MC3.hs +10/−10
• src/Mcmc/Algorithm/MHG.hs +2/−2
• src/Mcmc/Chain/Save.hs +11/−11
• src/Mcmc/Cycle.hs +9/−0
• src/Mcmc/Internal/Gamma.hs +2/−2
• src/Mcmc/Internal/Random.hs +2/−2
• src/Mcmc/Mcmc.hs +6/−6
• src/Mcmc/Monitor.hs +34/−34
• src/Mcmc/Monitor/Parameter.hs +3/−0
• src/Mcmc/Prior.hs +1/−1
• src/Mcmc/Proposal.hs +79/−86
• src/Mcmc/Proposal/Bactrian.hs +2/−2
• src/Mcmc/Proposal/Hamiltonian.hs +17/−22
• src/Mcmc/Proposal/Slide.hs +1/−1
• src/Mcmc/Settings.hs +3/−3

@@ -5,6 +5,13 @@ ## Unreleased changes  +## 0.6.2.0++-   Improve leapfrog integrator.+-   Update tooling.+-   Cleanup proposals.++ ## 0.6.1.0  -   Revamp Hamiltonian proposal (storable vectors).@@ -70,7 +77,7 @@ ## 0.2.4  -   **Change order of arguments for proposals**.--   'slideStem' was renamed to 'slideBranch'.+-   &rsquo;slideStem&rsquo; was renamed to &rsquo;slideBranch&rsquo;. -   Change ProposalSimple from newtype to type. -   Contravariant instances of parameter and batch monitors. Use `(>$<)` instead     of `(@.)` and `(@#)`.

@@ -1,6 +1,6 @@-cabal-version:  2.2+cabal-version:  3.0 name:           mcmc-version:        0.6.1.0+version:        0.6.2.0 synopsis:       Sample from a posterior using Markov chain Monte Carlo description:    Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme> category:       Math, Statistics@@ -116,7 +116,6 @@   build-depends:       base >=4.7 && <5     , hspec-    , log-domain     , mcmc     , mwc-random     , statistics

@@ -255,11 +255,11 @@   -- Only set the id for the cold chain.   | i == 0 = return $ MHG $ c {chainId = Just 0}   | 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-    -- pointing to the link of the cold chain, but this has no effect.-    t' <- pushT l t-    return $ MHG $ c {chainId = Just i, trace = t'}+      -- 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+      -- pointing to the link of the cold chain, but this has no effect.+      t' <- pushT l t+      return $ MHG $ c {chainId = Just i, trace = t'}   where     a' = setReciprocalTemperature prf lhf beta a     c = fromMHG a'@@ -291,11 +291,11 @@   | sp < 1 = error "mc3: The swap period must be strictly positive."   | 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-    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+      -- Split random number generators.+      gs <- V.fromList <$> splitGen n g+      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   where     n = fromNChains $ mc3NChains sMc3     sp = fromSwapPeriod $ mc3SwapPeriod sMc3

@@ -170,8 +170,8 @@ mhgAccept r g   | ln r >= 0.0 = return True   | otherwise = do-    b <- uniform g-    return $ b < exp (ln r)+      b <- uniform g+      return $ b < exp (ln r)  mhgPropose :: MHG a -> Proposal a -> IO (MHG a) mhgPropose (MHG c) p = do

@@ -91,19 +91,19 @@   IO (Chain a) fromSavedChain pr lh cc mn (SavedChain ci it i tr ac' g' ts)   | pr (state it) /= prior it =-    let msg =-          unlines-            [ "fromSave: Provided prior function does not match the saved prior.",-              "fromSave: Current prior:" <> show (prior it) <> ".",-              "fromSave: Given prior:" <> show (pr $ state it) <> "."-            ]-     in error msg+      let msg =+            unlines+              [ "fromSave: Provided prior function does not match the saved prior.",+                "fromSave: Current prior:" <> show (prior it) <> ".",+                "fromSave: Given prior:" <> show (pr $ state it) <> "."+              ]+       in error msg   | lh (state it) /= likelihood it =-    error "fromSave: Provided likelihood function does not match the saved likelihood."+      error "fromSave: Provided likelihood function does not match the saved likelihood."   | otherwise = do-    g <- loadGen g'-    tr' <- thawT tr-    return $ Chain ci it i tr' ac g i pr lh cc' mn+      g <- loadGen g'+      tr' <- thawT tr+      return $ Chain ci it i tr' ac g i pr lh cc' mn   where     ac = transformKeysA [0 ..] (ccProposals cc) ac'     tunePs mt p = case mt of

@@ -138,6 +138,15 @@   where     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+ -- | Calculate acceptance rates and auto tunes the 'Proposal's in the 'Cycle'. -- -- Do not change 'Proposal's that are not tuneable.

@@ -77,8 +77,8 @@ lgamma15_2 :: RealFloat a => a -> a -> a lgamma15_2 zm1 zm2 =   r * y + r-    * ( evaluatePolynomial (- zm2) tableLogGamma_15_2P-          / evaluatePolynomial (- zm2) tableLogGamma_15_2Q+    * ( evaluatePolynomial (-zm2) tableLogGamma_15_2P+          / evaluatePolynomial (-zm2) tableLogGamma_15_2Q       )   where     r = zm1 * zm2

@@ -33,8 +33,8 @@ splitGen n gen   | n <= 0 = return []   | otherwise = do-    seeds :: [V.Vector Word32] <- replicateM n $ uniformVector gen 256-    mapM initialize seeds+      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.

@@ -93,12 +93,12 @@   | n < 0 = error "mcmcIterate: Number of iterations is negative."   | n == 0 = return a   | otherwise = do-    e <- ask-    p <- sParallelizationMode . settings <$> ask-    -- NOTE: User interrupt is handled during iterations.-    a' <- liftIO $ catch (aIterate p a) (mcmcExceptionHandler e a)-    mcmcExecuteMonitors a'-    mcmcIterate (n -1) a'+      e <- ask+      p <- sParallelizationMode . settings <$> ask+      -- NOTE: User interrupt is handled during iterations.+      a' <- liftIO $ catch (aIterate p a) (mcmcExceptionHandler e a)+      mcmcExecuteMonitors a'+      mcmcIterate (n - 1) a'  mcmcNewRun :: Algorithm a => a -> MCMC a mcmcNewRun a = do

@@ -77,7 +77,7 @@ simpleMonitor p   | p < 1 = error "simpleMonitor: Monitor period must be 1 or larger."   | otherwise =-    Monitor (MonitorStdOut [] p) [] []+      Monitor (MonitorStdOut [] p) [] []  -- | Monitor to standard output; constructed with 'monitorStdOut'. data MonitorStdOut a = MonitorStdOut@@ -208,19 +208,19 @@ mfExec i (Link x p l) m   | i `mod` mfPeriod m /= 0 = return ()   | otherwise = case mfHandle m of-    Nothing ->-      error $-        "mfExec: No handle available for monitor with name "-          <> mfName m-          <> "."-    Just h ->-      BL.hPutStrLn h $-        mfRenderRow $-          BL.pack (show i) :-          renderLog p :-          renderLog l :-          renderLog (p * l) :-            [BB.toLazyByteString $ mpFunc mp x | mp <- mfParams m]+      Nothing ->+        error $+          "mfExec: No handle available for monitor with name "+            <> mfName m+            <> "."+      Just h ->+        BL.hPutStrLn h $+          mfRenderRow $+            BL.pack (show i) :+            renderLog p :+            renderLog l :+            renderLog (p * l) :+              [BB.toLazyByteString $ mpFunc mp x | mp <- mfParams m]  mfClose :: MonitorFile a -> IO () mfClose m = case mfHandle m of@@ -288,26 +288,26 @@ mbExec i t m   | (i `mod` mbSize m /= 0) || (i == 0) = return ()   | otherwise = case mbHandle m of-    Nothing ->-      error $-        "mbExec: No handle available for batch monitor with name "-          <> mbName m-          <> "."-    Just h -> do-      xs <- takeT (mbSize m) t-      let lps = VB.map prior xs-          lls = VB.map likelihood xs-          los = VB.zipWith (*) lps lls-          mlps = mean lps-          mlls = mean lls-          mlos = mean los-      BL.hPutStrLn h $-        mfRenderRow $-          BL.pack (show i) :-          renderLog mlps :-          renderLog mlls :-          renderLog mlos :-            [BB.toLazyByteString $ mbpFunc mbp (VB.map state xs) | mbp <- mbParams m]+      Nothing ->+        error $+          "mbExec: No handle available for batch monitor with name "+            <> mbName m+            <> "."+      Just h -> do+        xs <- takeT (mbSize m) t+        let lps = VB.map prior xs+            lls = VB.map likelihood xs+            los = VB.zipWith (*) lps lls+            mlps = mean lps+            mlls = mean lls+            mlos = mean los+        BL.hPutStrLn h $+          mfRenderRow $+            BL.pack (show i) :+            renderLog mlps :+            renderLog mlls :+            renderLog mlos :+              [BB.toLazyByteString $ mbpFunc mbp (VB.map state xs) | mbp <- mbParams m]  mbClose :: MonitorBatch a -> IO () mbClose m = case mbHandle m of

@@ -26,6 +26,9 @@ import qualified Data.Double.Conversion.ByteString as BC import Data.Functor.Contravariant +-- XXX: 'MonitorParameter' has a drawback. Extracting and monitoring multiple+-- parameters in one go is impossible.+ -- | Instruction about a parameter to monitor. -- -- Convert a parameter monitor from one data type to another with '(>$<)'.

@@ -161,7 +161,7 @@ normal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a normal m s x   | s <= 0 = error "normal: Standard deviation is zero or negative."-  | otherwise = Exp $ (- xm * xm / (2 * s * s)) - denom+  | otherwise = Exp $ (-xm * xm / (2 * s * s)) - denom   where     xm = x - m     denom = mLnSqrt2Pi + log s

@@ -14,7 +14,7 @@ -- -- Creation date: Wed May 20 13:42:53 2020. module Mcmc.Proposal-  ( -- * Proposals and types+  ( -- * Proposals     PName (..),     PDescription (..),     PWeight (fromPWeight),@@ -28,18 +28,21 @@     liftProposal,     liftProposalWith,     ProposalSimple,+    createProposal,++    -- * Tuners     Tuner (..),     Tune (..),     TuningParameter,+    TuningFunction,     AuxiliaryTuningParameters,+    AuxiliaryTuningFunction,     defaultTuningFunctionWith,     noTuningFunction,     noAuxiliaryTuningFunction,-    createProposal,     tuningParameterMin,     tuningParameterMax,     tuneWithTuningParameters,-    tuneWithChainParameters,     getOptimalRate,      -- * Output@@ -55,7 +58,6 @@ import qualified Data.Vector as VB import qualified Data.Vector.Unboxed as VU import Lens.Micro-import Lens.Micro.Extras import Mcmc.Acceptance import Mcmc.Internal.ByteString import Numeric.Log hiding (sum)@@ -71,27 +73,30 @@   deriving (Show, Eq, Ord)  -- | The positive weight determines how often a 'Proposal' is executed per--- iteration of the Markov chain.+-- iteration of the Markov chain. Abstract data type; for construction, see+-- 'pWeight'. newtype PWeight = PWeight {fromPWeight :: Int}   deriving (Show, Eq, Ord)  -- | Check if the weight is positive.+--+-- Call 'error' if weight is zero or negative. pWeight :: Int -> PWeight pWeight n-  | n <= 0 = error "pWeight: Proposal weight is zero or negative."+  | n <= 0 = error $ "pWeight: Proposal weight is zero or negative: " <> show n <> "."   | otherwise = PWeight n  -- | Proposal dimension. -- -- The number of affected, independent parameters. ----- The optimal acceptance rate of low dimensional proposals is higher than for--- high dimensional ones.+-- Usually, the optimal acceptance rate of low dimensional proposals is higher+-- than for high dimensional ones. However, this is not always true (see below). ----- Optimal acceptance rates are still subject to controversies. As far as I--- know, research has focused on random walk proposals with multivariate normal--- distributions of dimension @d@. In this case, the following acceptance rates--- are desired:+-- Further, optimal acceptance rates are still subject to controversies. To my+-- knowledge, research has focused on random walk proposals with multivariate+-- normal distributions of dimension @d@. In this case, the following acceptance+-- rates are desired: -- -- - one dimension: 0.44 (numerical results); --@@ -130,8 +135,9 @@ -- state space @a@. Essentially, it is a probability mass or probability density -- conditioned on the current state (i.e., a Markov kernel). ----- A 'Proposal' may be tuneable in that it contains information about how to enlarge--- or shrink the step size to tune the acceptance rate.+-- A 'Proposal' may be tuneable in that it contains information about how to+-- enlarge or shrink the proposal size to decrease or increase the acceptance+-- rate. -- -- Predefined proposals are provided. To create custom proposals, one may use -- the convenience function 'createProposal'.@@ -165,8 +171,8 @@ -- See also 'Jacobian'. -- -- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,--- it is more declarative to have them separate. It is a constant reminder: Is--- the Jacobian modifier different from 1.0?+-- it is more declarative to have them separate. Like so, we are constantly+-- reminded: Is the Jacobian modifier different from 1.0? type KernelRatio = Log Double  -- | Absolute value of the determinant of the Jacobian matrix.@@ -189,7 +195,7 @@ -- scaleFirstEntryOfTuple = _1 @~ scale -- @ ----- See also 'liftProposal' and 'liftProposalWith'.+-- See also 'liftProposalWith'. infixl 7 @~  (@~) :: Lens' b a -> Proposal a -> Proposal b@@ -201,14 +207,13 @@  -- | Lift a proposal from one data type to another. ----- A function to calculate the Jacobian has to be provided (see also--- 'liftProposal').+-- A function to calculate the Jacobian has to be provided (but see '(@~)'). -- -- 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 w s t) =-  Proposal n r d w (liftProposalSimple jf l s) (liftTuner jf l <$> t)+  Proposal n r d w (liftProposalSimpleWith jf l s) (liftTunerWith jf l <$> t)  -- | Simple proposal without tuning information. --@@ -231,8 +236,8 @@ type ProposalSimple a = a -> GenIO -> IO (a, KernelRatio, Jacobian)  -- Lift a simple proposal from one data type to another.-liftProposalSimple :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b-liftProposalSimple jf l s = s'+liftProposalSimpleWith :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b+liftProposalSimpleWith jf l s = s'   where     s' y g = do       (x', r, j) <- s (y ^. l) g@@ -242,19 +247,40 @@           j' = j * jyx / jxy       return (y', r, j') +-- | Create a proposal with a single tuning parameter.+--+-- Proposals with arbitrary 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) ->+  -- | Dimension.+  PDimension ->+  -- | Name.+  PName ->+  -- | Weight.+  PWeight ->+  -- | Activate tuning?+  Tune ->+  Proposal a+createProposal r f d n w Tune =+  Proposal n r d w (f 1.0) (Just tuner)+  where+    fT = defaultTuningFunctionWith d+    fTs = noAuxiliaryTuningFunction+    g t _ = Right $ f t+    tuner = Tuner 1.0 fT VU.empty fTs g+createProposal r f d n w NoTune =+  Proposal n r d w (f 1.0) Nothing+ -- | Required information to tune 'Proposal's. data Tuner a = Tuner   { tTuningParameter :: TuningParameter,-    -- | Instruction about how to compute new tuning parameter from a given-    -- acceptance rate and the old tuning parameter.-    tComputeTuningParameter :: AcceptanceRate -> TuningParameter -> TuningParameter,+    tTuningFunction :: TuningFunction,     tAuxiliaryTuningParameters :: AuxiliaryTuningParameters,-    -- | Instruction about how to compute new auxiliary tuning parameters from a-    -- given trace and the old auxiliary tuning parameters.-    tComputeAuxiliaryTuningParameters ::-      VB.Vector a ->-      AuxiliaryTuningParameters ->-      AuxiliaryTuningParameters,+    tAuxiliaryTuningFunction :: AuxiliaryTuningFunction a,     -- | Given the tuning parameter, and the auxiliary tuning parameters, get     -- the tuned simple proposal.     --@@ -267,11 +293,11 @@   }  -- Lift tuner from one data type to another.-liftTuner :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b-liftTuner jf l (Tuner p fP ps fPs g) = Tuner p fP ps fPs' g'+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'   where-    fPs' = fPs . VB.map (view l)-    g' x xs = liftProposalSimple jf l <$> g x xs+    fPs' = fPs . VB.map (^. l)+    g' x xs = liftProposalSimpleWith jf l <$> g x xs  -- | Tune proposal? data Tune = Tune | NoTune@@ -283,75 +309,49 @@ -- 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+ -- | Auxiliary tuning parameters; vector may be empty.+--+-- Auxiliary tuning parameters are not shown in proposal summaries. 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+ -- | Default tuning function. -- -- Subject to change. defaultTuningFunctionWith ::   -- Optimal acceptance rate.   PDimension ->-  AcceptanceRate ->-  TuningParameter ->-  TuningParameter+  TuningFunction defaultTuningFunctionWith d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t --- | Do no tune.+-- | Do not tune. -- -- Useful if auxiliary tuning parameters are tuned, but not the main tuning -- parameter.-noTuningFunction :: AcceptanceRate -> TuningParameter -> TuningParameter+noTuningFunction :: TuningFunction noTuningFunction _ = id  -- | Do not tune auxiliary parameters.-noAuxiliaryTuningFunction :: VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters+noAuxiliaryTuningFunction :: AuxiliaryTuningFunction a noAuxiliaryTuningFunction _ = id --- | Create a proposal with a single tuning parameter.------ Proposals with arbitrary 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) ->-  -- | Dimension.-  PDimension ->-  -- | Name.-  PName ->-  -- | Weight.-  PWeight ->-  -- | Activate tuning?-  Tune ->-  Proposal a-createProposal r f d n w Tune =-  Proposal n r d w (f 1.0) (Just tuner)-  where-    fT = defaultTuningFunctionWith d-    fTs = noAuxiliaryTuningFunction-    g t _ = Right $ f t-    tuner = Tuner 1.0 fT VU.empty fTs g-createProposal r f d n w NoTune =-  Proposal n r d w (f 1.0) Nothing- -- 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 -- tuning parameters are very different from one, a different base proposal -- should be chosen. --- | Minimal tuning parameter; @1e-5@, subject to change.------ >>> tuningParameterMin--- 1e-5+-- | Minimal tuning parameter; subject to change. tuningParameterMin :: TuningParameter tuningParameterMin = 1e-5 --- | Maximal tuning parameter; @1e3@, subject to change.------ >>> tuningParameterMax--- 1e3+-- | Maximal tuning parameter; subject to change. tuningParameterMax :: TuningParameter tuningParameterMax = 1e3 @@ -368,6 +368,8 @@ -- - the 'Proposal' is not tuneable; -- -- - the auxiliary tuning parameters are invalid.+--+-- Used by 'Mcmc.Chain.Save.fromSavedChain'. tuneWithTuningParameters ::   TuningParameter ->   AuxiliaryTuningParameters ->@@ -383,15 +385,6 @@      in case psE of           Left err -> Left $ "tune: " <> err           Right ps -> Right $ p {prSimple = ps, prTuner = Just $ Tuner t'' fT ts fTs g}---- | 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  -- | See 'PDimension'. getOptimalRate :: PDimension -> Double

@@ -40,14 +40,14 @@       d = normalDistr mn sd   x <- genContVar d g   b <- bernoulli 0.5 g-  return $ if b then x else - x+  return $ if b then x else -x  logDensityBactrian :: SpikeParameter -> StandardDeviation Double -> Double -> Log Double logDensityBactrian m s x = Exp $ log $ kernel1 + kernel2   where     mn = m * s     sd = sqrt (1 - m * m) * s-    dist1 = normalDistr (- mn) sd+    dist1 = normalDistr (-mn) sd     dist2 = normalDistr mn sd     kernel1 = density dist1 x     kernel2 = density dist2 x

@@ -333,14 +333,14 @@   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+      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+      phiL = leapfrogStepMomenta (0.5 * eps) grad thetaL phiLM1Half   return (thetaL, phiL)   where     valF x = case mVal of@@ -350,14 +350,13 @@     go 0 (Just (t, p)) = Just (t, p)     go n (Just (t, p)) =       let t' = leapfrogStepPositions hMassesInvEps t p-          p' = leapfrogStepMomenta 1.0 eps grad t' p+          p' = leapfrogStepMomenta eps grad t' p           r = (,p') <$> valF t'        in go (n - 1) r  leapfrogStepMomenta ::-  -- Size of step (half or full step).+  -- Step size (usually a multiple of the leapfrog scaling factor).   Double ->-  LeapfrogScalingFactor ->   Gradient Positions ->   -- Current positions.   Positions ->@@ -365,7 +364,7 @@   Momenta ->   -- New momenta.   Momenta-leapfrogStepMomenta xi eps grad theta phi = phi + L.scale (xi * eps) (grad theta)+leapfrogStepMomenta eps grad theta phi = phi + L.scale eps (grad theta)  leapfrogStepPositions ::   HMassesInvEps ->@@ -394,14 +393,14 @@   Either String (HSettings a) hTuningParametersToSettings s t ts   | nTsNotOK =-    Left "hTuningParametersToSettings: Auxiliary variables do not have correct dimension."+      Left "hTuningParametersToSettings: Auxiliary variables do not have correct dimension."   | otherwise =-    Right $-      s-        { hMasses = msTuned,-          hLeapfrogTrajectoryLength = lTuned,-          hLeapfrogScalingFactor = eTuned-        }+      Right $+        s+          { hMasses = msTuned,+            hLeapfrogTrajectoryLength = lTuned,+            hLeapfrogScalingFactor = eTuned+          }   where     ms = hMasses s     d = L.rows $ L.unSym ms@@ -519,16 +518,14 @@ tuneDiagonalMassesOnly ::   Int ->   (a -> Positions) ->-  VB.Vector a ->-  AuxiliaryTuningParameters ->-  AuxiliaryTuningParameters+  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+      -- 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.@@ -550,9 +547,7 @@ tuneAllMasses ::   Int ->   (a -> Positions) ->-  VB.Vector a ->-  AuxiliaryTuningParameters ->-  AuxiliaryTuningParameters+  AuxiliaryTuningFunction a tuneAllMasses dim toVec xs ts   -- If not enough data is available, do not tune.   | VB.length xs < samplesMinDiagonal = ts

@@ -67,7 +67,7 @@ -- The actual proposal with tuning parameter. slideUniformSimple :: Size -> TuningParameter -> ProposalSimple Double slideUniformSimple d t =-  genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing+  genericContinuous (uniformDistr (-t * d) (t * d)) (+) Nothing Nothing  -- | See 'slide'. --

@@ -316,12 +316,12 @@   | not $ burnInValid bi = serr $ "Burn in setting invalid: " <> show bi <> "."   | fromIterations i < 0 = serr "Number of iterations is negative."   | burnInIterations bi + fromIterations i - iCurrent < 0 =-    serr "Current iteration is larger than the total number of iterations."+      serr "Current iteration is larger than the total number of iterations."   | not $ validTraceLength tl = serr $ "Trace length invalid: " <> show tl <> "."   | iCurrent /= 0 && em /= Continue =-    serr "Current iteration is non-zero but execution mode is not 'Continue'."+      serr "Current iteration is non-zero but execution mode is not 'Continue'."   | iCurrent == 0 && em == Continue =-    serr "Current iteration is zero but execution mode is 'Continue'."+      serr "Current iteration is zero but execution mode is 'Continue'."   | otherwise = return ()   where     serr = settingsError s iCurrent