mcmc 0.4.0.0 → 0.5.0.0

raw patch · 39 files changed

+1808/−849 lines, 39 filesdep ~aesondep ~bytestringdep ~circularPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: aeson, bytestring, circular, containers, criterion, data-default, deepseq, directory, dirichlet, double-conversion, hspec, log-domain, microlens, monad-parallel, mwc-random, pretty-show, primitive, statistics, time, transformers, vector, zlib

API changes (from Hackage documentation)

Files

• ChangeLog.md +10/−0
• bench/Bench.hs +1/−1
• bench/Normal.hs +12/−8
• bench/Poisson.hs +5/−5
• mcmc.cabal +36/−31
• src/Mcmc.hs +28/−15
• src/Mcmc/Algorithm.hs +5/−1
• src/Mcmc/Algorithm/MC3.hs +86/−30
• src/Mcmc/Algorithm/MHG.hs +283/−0
• src/Mcmc/Algorithm/Metropolis.hs +0/−254
• src/Mcmc/Chain/Chain.hs +12/−17
• src/Mcmc/Chain/Link.hs +7/−5
• src/Mcmc/Chain/Save.hs +20/−16
• src/Mcmc/Chain/Trace.hs +1/−1
• src/Mcmc/Environment.hs +53/−16
• src/Mcmc/Internal/ByteString.hs +1/−1
• src/Mcmc/Internal/Random.hs +1/−1
• src/Mcmc/Internal/Shuffle.hs +1/−1
• src/Mcmc/Likelihood.hs +29/−0
• src/Mcmc/Logger.hs +170/−0
• src/Mcmc/MarginalLikelihood.hs +404/−0
• src/Mcmc/Mcmc.hs +87/−133
• src/Mcmc/Monitor.hs +31/−22
• src/Mcmc/Monitor/Log.hs +1/−1
• src/Mcmc/Monitor/Parameter.hs +1/−1
• src/Mcmc/Monitor/ParameterBatch.hs +1/−1
• src/Mcmc/Monitor/Time.hs +1/−1
• src/Mcmc/Posterior.hs +20/−0
• src/Mcmc/Prior.hs +69/−42
• src/Mcmc/Proposal.hs +229/−124
• src/Mcmc/Proposal/Bactrian.hs +24/−33
• src/Mcmc/Proposal/Generic.hs +2/−2
• src/Mcmc/Proposal/Scale.hs +9/−13
• src/Mcmc/Proposal/Simplex.hs +9/−7
• src/Mcmc/Proposal/Slide.hs +16/−29
• src/Mcmc/Settings.hs +72/−22
• src/Mcmc/Statistics/Types.hs +56/−0
• test/Mcmc/ProposalSpec.hs +7/−7
• test/Mcmc/SaveSpec.hs +8/−8

@@ -5,6 +5,16 @@ ## Unreleased changes  +## 0.5.0.0++-   Marginal likelihood estimation using thermodynamic integration or stepping+    stone sampling.+-   Various changes of function names (e.g., metropologis -> mhg).+-   Updated examples.+-   Proper but minimal logging framework.+-   Various other changes.++ ## 0.4.0.0  -   Greatly improve documentation.

@@ -1,6 +1,6 @@ -- | -- Module      :  Main--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Normal -- Description :  Benchmark Metropolis-Hastings-Green algorithm--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -31,7 +31,7 @@ lh = normal trueMean stdDev  cc :: Cycle Double-cc = cycleFromList [slideSymmetric 1.0 (PName "Medium") (PWeight 1) Tune]+cc = cycleFromList [slideSymmetric 1.0 (PName "Medium") (pWeight 1) Tune]  mons :: [MonitorParameter Double] mons = [monitorDouble "mu"]@@ -52,14 +52,15 @@           Overwrite           Sequential           NoSave+          LogStdOutOnly           Quiet-  a <- mhg noPrior lh cc mon 0 g+  a <- mhg noPrior lh cc mon TraceAuto 0 g   void $ mcmc s a  ccLarge :: Cycle Double ccLarge =   cycleFromList-    [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (PWeight 1) Tune | i <- [0 .. 100 :: Int]]+    [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (pWeight 1) Tune | i <- [0 .. 100 :: Int]]  -- Should have the same run time as 'normalSlide'. normalLargeCycleBench :: GenIO -> IO ()@@ -72,12 +73,13 @@           Overwrite           Sequential           NoSave+          LogStdOutOnly           Quiet-  a <- mhg noPrior lh ccLarge mon 0 g+  a <- mhg noPrior lh ccLarge mon TraceAuto 0 g   void $ mcmc s a  ccBactrian :: Cycle Double-ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (PWeight 1) Tune]+ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (pWeight 1) Tune]  normalBactrianBench :: GenIO -> IO () normalBactrianBench g = do@@ -89,8 +91,9 @@           Overwrite           Sequential           NoSave+          LogStdOutOnly           Quiet-  a <- mhg noPrior lh ccBactrian mon 0 g+  a <- mhg noPrior lh ccBactrian mon TraceAuto 0 g   void $ mcmc s a  normalMC3 :: GenIO -> Int -> IO ()@@ -103,7 +106,8 @@           Overwrite           Sequential           NoSave+          LogStdOutOnly           Quiet       mc3S = MC3Settings (NChains n) (SwapPeriod 2) (NSwaps 1)-  a <- mc3 mc3S noPrior lh cc mon 0 g+  a <- mc3 mc3S noPrior lh cc mon TraceAuto 0 g   void $ mcmc mcmcS a

@@ -1,7 +1,7 @@ -- | -- Module      :  Poisson -- Description :  Poisson regression model for airline fatalities--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -19,7 +19,6 @@ import Control.Monad import Lens.Micro import Mcmc-import Numeric.Log hiding (sum) import System.Random.MWC  type I = (Double, Double)@@ -42,10 +41,10 @@ lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]  proposalAlpha :: Proposal I-proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (PWeight 1) NoTune+proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune  proposalBeta :: Proposal I-proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (PWeight 1) NoTune+proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune  proposals :: Cycle I proposals = cycleFromList [proposalAlpha, proposalBeta]@@ -75,6 +74,7 @@           Overwrite           Sequential           NoSave+          LogStdOutOnly           Quiet-  a <- mhg noPrior lh proposals mon initial g+  a <- mhg noPrior lh proposals mon TraceAuto initial g   void $ mcmc s a

@@ -1,8 +1,8 @@ cabal-version:      2.2 name:               mcmc-version:            0.4.0.0+version:            0.5.0.0 license:            GPL-3.0-or-later-copyright:          Dominik Schrempf (2020)+copyright:          Dominik Schrempf (2021) maintainer:         dominik.schrempf@gmail.com author:             Dominik Schrempf homepage:           https://github.com/dschrempf/mcmc#readme@@ -25,19 +25,23 @@     exposed-modules:         Mcmc         Mcmc.Algorithm-        Mcmc.Algorithm.Metropolis         Mcmc.Algorithm.MC3+        Mcmc.Algorithm.MHG         Mcmc.Chain.Chain         Mcmc.Chain.Link         Mcmc.Chain.Save         Mcmc.Chain.Trace         Mcmc.Environment+        Mcmc.Likelihood+        Mcmc.Logger+        Mcmc.MarginalLikelihood         Mcmc.Mcmc         Mcmc.Monitor         Mcmc.Monitor.Log         Mcmc.Monitor.Parameter         Mcmc.Monitor.ParameterBatch         Mcmc.Monitor.Time+        Mcmc.Posterior         Mcmc.Prior         Mcmc.Proposal         Mcmc.Proposal.Bactrian@@ -46,6 +50,7 @@         Mcmc.Proposal.Slide         Mcmc.Proposal.Simplex         Mcmc.Settings+        Mcmc.Statistics.Types      hs-source-dirs:   src     other-modules:@@ -58,27 +63,27 @@     default-language: Haskell2010     ghc-options:      -Wall -Wunused-packages     build-depends:-        aeson >=1.5.4.1 && <1.6,+        aeson >=1.5.6.0,         base >=4.7 && <5,-        bytestring >=0.10.10.0 && <0.11,-        circular >=0.3.1.1 && <0.4,-        containers >=0.6.2.1 && <0.7,-        data-default >=0.7.1.1 && <0.8,-        deepseq >=1.4.4.0 && <1.5,-        directory >=1.3.6.0 && <1.4,-        dirichlet >=0.1.0.0 && <0.2,-        double-conversion >=2.0.2.0 && <2.1,-        log-domain ==0.13.*,-        microlens >=0.4.11.2 && <0.5,-        mwc-random >=0.14.0.0 && <0.15,-        monad-parallel >=0.7.2.3 && <0.8,-        pretty-show ==1.10.*,-        primitive >=0.7.1.0 && <0.8,-        statistics >=0.15.2.0 && <0.16,-        time >=1.9.3 && <1.10,-        transformers >=0.5.6.2 && <0.6,-        vector >=0.12.1.2 && <0.13,-        zlib >=0.6.2.2 && <0.7+        bytestring >=0.10.12.0,+        circular >=0.4.0.0,+        containers >=0.6.2.1,+        data-default >=0.7.1.1,+        deepseq >=1.4.4.0,+        directory >=1.3.6.0,+        dirichlet >=0.1.0.4,+        double-conversion >=2.0.2.0,+        log-domain >=0.13.1,+        microlens >=0.4.12.0,+        mwc-random >=0.15.0.1,+        monad-parallel >=0.7.2.4,+        pretty-show >=1.10,+        primitive >=0.7.1.0,+        statistics >=0.15.2.0,+        time >=1.9.3,+        transformers >=0.5.6.2,+        vector >=0.12.3.0,+        zlib >=0.6.2.3  test-suite mcmc-test     type:             exitcode-stdio-1.0@@ -93,11 +98,11 @@     ghc-options:      -Wall -Wunused-packages     build-depends:         base >=4.7 && <5,-        hspec >=2.7.4 && <2.8,-        log-domain ==0.13.*,+        hspec >=2.7.10,+        log-domain >=0.13.1,         mcmc -any,-        mwc-random >=0.14.0.0 && <0.15,-        statistics >=0.15.2.0 && <0.16+        mwc-random >=0.15.0.1,+        statistics >=0.15.2.0  benchmark mcmc-bench     type:             exitcode-stdio-1.0@@ -112,8 +117,8 @@     ghc-options:      -Wall -Wunused-packages     build-depends:         base >=4.7 && <5,-        criterion >=1.5.7.0 && <1.6,-        log-domain ==0.13.*,+        criterion >=1.5.9.0,+        log-domain >=0.13.1,         mcmc -any,-        microlens >=0.4.11.2 && <0.5,-        mwc-random >=0.14.0.0 && <0.15+        microlens >=0.4.12.0,+        mwc-random >=0.15.0.1

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc -- Description :  Markov chain Monte Carlo samplers, batteries included--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -76,9 +76,13 @@     -- Gibbs sampler, or a random sequence scan Gibbs sampler, respectively. See     -- 'Order'.     ---    -- Note that it is of utter importance that the given 'Cycle' enables-    -- traversal of the complete state space. Otherwise, the Markov chain will-    -- not converge to the correct stationary posterior distribution.+    -- Notes:+    -- - It is important that the given 'Cycle' enables traversal of the+    --   complete state space. Otherwise, the Markov chain will not converge to+    --   the correct stationary posterior distribution.+    -- - Be careful when assigning proposals because acceptance ratios may have+    --   to be amended using Jacobians. Please see an [example involving a pair+    --   of numbers](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs).     --     -- Proposals are named according to what they do, i.e., how they change the     -- state of a Markov chain, and not according to the intrinsically used@@ -114,8 +118,12 @@     -- http://dx.doi.org/10.1093/sysbio/syw021     PName (..),     PWeight (..),+    pWeight,     Proposal,+    JacobianFunction,     (@~),+    liftProposal,+    liftProposalWith,     Tune (..),     scale,     scaleUnbiased,@@ -157,36 +165,39 @@     monitorFile,     MonitorBatch,     monitorBatch,+    simpleMonitor,     module Mcmc.Monitor.Parameter,     module Mcmc.Monitor.ParameterBatch, -    -- * Prior distributions--    -- | Convenience functions for computing priors.+    -- * Prior, likelihood, and posterior values and functions     module Mcmc.Prior,+    module Mcmc.Likelihood,      -- * MCMC samplers     mcmc,     mcmcContinue,-     -- | See also 'settingsLoad', 'mhgLoad', and 'mc3Load'.      -- * Algorithms-    module Mcmc.Algorithm.Metropolis,+    module Mcmc.Algorithm.MHG,     module Mcmc.Algorithm.MC3, +    -- * Marginal likelihood calculation+    module Mcmc.MarginalLikelihood, -    -- * Useful type synonyms-    PriorFunction,-    noPrior,-    LikelihoodFunction,-    noLikelihood,+    -- * Useful types+    module Mcmc.Statistics.Types,++    -- * Useful re-exports+    module Numeric.Log,   ) where  import Mcmc.Algorithm.MC3-import Mcmc.Algorithm.Metropolis+import Mcmc.Algorithm.MHG import Mcmc.Chain.Chain+import Mcmc.Likelihood+import Mcmc.MarginalLikelihood import Mcmc.Mcmc import Mcmc.Monitor import Mcmc.Monitor.Parameter@@ -198,3 +209,5 @@ import Mcmc.Proposal.Simplex import Mcmc.Proposal.Slide import Mcmc.Settings+import Mcmc.Statistics.Types+import Numeric.Log hiding (sum)

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Algorithm -- Description :  Algortihms for Markov chain Monte Carlo samplers--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -25,6 +25,10 @@    -- | Current iteration.   aIteration :: a -> Int++  -- | Check if the current state is invalid. At the moment this should just+  -- check if the posterior probability is zero or NaN.+  aIsInValidState :: a -> Bool    -- | Sample the next state.   aIterate :: ParallelizationMode -> a -> IO a

@@ -5,7 +5,7 @@ -- | -- Module      :  Mcmc.Algorithm.MC3 -- Description :  Metropolis-coupled Markov chain Monte Carlo algorithm--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -19,8 +19,8 @@ -- Also known as parallel tempering. -- -- Like any other parallel MCMC algorithm, the 'MC3' algorithm is essentially an--- 'Mcmc.Algorithm.Metropolis.MHG' algorithm on the product space of all--- parallel chains.+-- 'Mcmc.Algorithm.MHG.MHG' algorithm on the product space of all parallel+-- chains. -- -- For example, see --@@ -64,7 +64,7 @@ import Data.Word -- import Debug.Trace hiding (trace) import Mcmc.Algorithm-import Mcmc.Algorithm.Metropolis+import Mcmc.Algorithm.MHG import Mcmc.Chain.Chain import Mcmc.Chain.Link import Mcmc.Chain.Save@@ -72,11 +72,14 @@ import Mcmc.Internal.Random import Mcmc.Internal.Shuffle import Mcmc.Monitor+import Mcmc.Posterior import Mcmc.Proposal import Mcmc.Settings import Numeric.Log hiding (sum) import System.Random.MWC +-- import Debug.Trace hiding (trace)+ -- | Total number of parallel chains. -- -- Must be two or larger.@@ -146,9 +149,17 @@   SavedMC3 a ->   IO (MC3 a) fromSavedMC3 pr lh cc mn (SavedMC3 s scs bs i ac g') = do-  mhgs <- V.mapM (fmap MHG . fromSavedChain pr lh cc mn) scs+  mhgs <-+    V.fromList+      <$> sequence+        [ MHG <$> fromSavedChain pf lf cc mn sc+          | (sc, pf, lf) <- zip3 (V.toList scs) prs lhs+        ]   g <- loadGen g'   return $ MC3 s mhgs bs i ac g+  where+    prs = map (heatFunction pr) $ U.toList bs+    lhs = map (heatFunction lh) $ U.toList bs  -- | The MC3 algorithm. data MC3 a = MC3@@ -167,6 +178,7 @@ instance ToJSON a => Algorithm (MC3 a) where   aName = const "Metropolis-coupled Markov chain Monte Carlo (MC3)"   aIteration = mc3Iteration+  aIsInValidState = mc3IsInValidState   aIterate = mc3Iterate   aAutoTune = mc3AutoTune   aResetAcceptance = mc3ResetAcceptance@@ -177,6 +189,20 @@   aCloseMonitors = mc3CloseMonitors   aSave = mc3Save +heatFunction ::+  -- Cold Function.+  (a -> Log Double) ->+  -- Reciprocal temperature.+  Double ->+  -- The heated prior or likelihood function+  (a -> Log Double)+heatFunction f b+  | b <= 0 = error "heatFunction: Reciprocal temperature is zero or negative."+  | b == 1.0 = f+  | otherwise = (** b') . f+  where+    b' = Exp $ log b+ --  The prior and likelihood values of the current link are updated. -- -- NOTE: The trace is not changed! In particular, the prior and likelihood@@ -191,7 +217,7 @@   Double ->   MHG a ->   MHG a-setReciprocalTemperature prf lhf beta a =+setReciprocalTemperature coldPrf coldLhf b a =   MHG $     c       { priorFunction = prf',@@ -200,15 +226,14 @@       }   where     c = fromMHG a-    b' = Exp $ log beta     -- We need twice the amount of computations compared to taking the power     -- after calculating the posterior (pr x * lh x) ** b'. However, I don't     -- think this is a serious problem.     --     -- To minimize computations, it is key to avoid modification of the     -- reciprocal temperature for the cold chain.-    prf' = (** b') . prf-    lhf' = (** b') . lhf+    prf' = heatFunction coldPrf b+    lhf' = heatFunction coldLhf b     x = state $ link c  initMHG ::@@ -224,14 +249,14 @@   IO (MHG a) initMHG prf lhf i beta a   | i < 0 = error "initMHG: Chain index negative."-  -- Do not temper with the cold chain.-  | i == 0 = return a+  -- 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 = i, trace = t'}+    return $ MHG $ c {chainId = Just i, trace = t'}   where     a' = setReciprocalTemperature prf lhf beta a     c = fromMHG a'@@ -254,17 +279,18 @@   LikelihoodFunction a ->   Cycle a ->   Monitor a ->-  a ->+  TraceLength ->+  InitialState a ->   GenIO ->   IO (MC3 a)-mc3 s pr lh cc mn i0 g+mc3 s pr lh cc mn tr i0 g   | n < 2 = error "mc3: The number of chains must be two or larger."   | 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 pr lh cc mn i0) gs+    cs <- V.mapM (mhg pr lh cc mn tr i0) gs     hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs     return $ MC3 s hcs bs 0 (emptyA [0 .. n - 2]) g   where@@ -276,10 +302,10 @@     --     -- NOTE: Have to 'take n' elements, because vectors are not as lazy as     -- lists.-    bs = U.fromList $ take n $ iterate (* 0.92) 1.0+    bs = U.fromList $ take n $ iterate (* 0.97) 1.0  mc3Fn :: AnalysisName -> FilePath-mc3Fn (AnalysisName nm) = nm ++ ".mc3"+mc3Fn (AnalysisName nm) = nm ++ ".mcmc.mc3"  -- | Save an MC3 algorithm. mc3Save ::@@ -315,7 +341,7 @@   -- Index j>=0, j/=i of right chain.   Int ->   MHGChains a ->-  (MHGChains a, Log Double)+  (MHGChains a, Posterior) swapWith i j xs   | i < 0 = error "swapWith: Left index is negative."   | j < 0 = error "swapWith: Right index is negative."@@ -356,12 +382,28 @@   Int ->   IO (MC3 a) mc3ProposeSwap a i = do+  let cs = mc3MHGChains a+  -- -- Debug.+  -- prL = prior $ link $ fromMHG $ cs V.! i+  -- prR = prior $ link $ fromMHG $ cs V.! (i+1)+  -- lhL = likelihood $ link $ fromMHG $ cs V.! i+  -- lhR = likelihood $ link $ fromMHG $ cs V.! (i+1)   -- 1. Sample new state and get the Metropolis ratio.-  let (!y, !r) = swapWith i (i + 1) $ mc3MHGChains a+  let (!y, !r) = swapWith i (i + 1) cs   -- 2. Accept or reject.   accept <- mhgAccept r g   if accept     then do+      -- -- Debug.+      -- traceIO $ "Swap accepted: " <> show i <> " <-> " <> show (i+1)+      -- let prL' = prior $ link $ fromMHG $ y V.! i+      --     prR' = prior $ link $ fromMHG $ y V.! (i+1)+      --     lhL' = likelihood $ link $ fromMHG $ y V.! i+      --     lhR' = likelihood $ link $ fromMHG $ y V.! (i+1)+      -- traceIO $ "Log priors (left, right, before swap): " <> show (ln prL) <> " " <> show (ln prR)+      -- traceIO $ "Log priors (left, right, after swap): " <> show (ln prL') <> " " <> show (ln prR')+      -- traceIO $ "Log likelihoods (left, right, before swap): " <> show (ln lhL) <> " " <> show (ln lhR)+      -- traceIO $ "Log likelihood (left, right, after swap): " <> show (ln lhL') <> " " <> show (ln lhR')       let !ac' = pushA i True (mc3SwapAcceptance a)       return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'}     else do@@ -370,9 +412,15 @@   where     g = mc3Generator a +mc3IsInValidState :: ToJSON a => MC3 a -> Bool+mc3IsInValidState a = V.any aIsInValidState mhgs+  where+    mhgs = mc3MHGChains a+ -- TODO: Splimix. 'mc3Iterate' is actually not parallel, but concurrent because -- of the IO constraint. Use pure parallel code when we have a pure generator.--- However, we have honor the mutable traces.+--+-- However, we to have honor the mutable traces. mc3Iterate ::   ToJSON a =>   ParallelizationMode ->@@ -399,7 +447,7 @@     Parallel ->       -- See 'Control.Monad.Parallel' of package 'monad-parallel'. Go via a       -- list, and use 'forkIO'.-      V.fromList <$> P.mapM (aIterate pm) (V.toList (mc3MHGChains a'))+      V.fromList <$> P.mapM (aIterate pm) (V.toList $ mc3MHGChains a')   let i = mc3Iteration a'   return $ a' {mc3MHGChains = mhgs, mc3Iteration = succ i} @@ -432,7 +480,7 @@     mhgs' = V.map aAutoTune mhgs     -- 2. Auto tune temperatures.     optimalRate = getOptimalRate PDimensionUnknown-    currentRates = acceptanceRates $ mc3SwapAcceptance a+    mCurrentRates = acceptanceRates $ mc3SwapAcceptance 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@@ -442,7 +490,9 @@     --     -- The factor (1/2) was determined by a few tests and is otherwise     -- absolutely arbitrary.-    xi i = exp $ (/ 2) $ (currentRates M.! i) - optimalRate+    xi i = case mCurrentRates M.! i of+      Nothing -> 1.0+      Just currentRate -> exp $ (/ 2) $ currentRate - optimalRate     bs = mc3ReciprocalTemperatures a     n = fromNChains $ mc3NChains $ mc3Settings a     -- Do not change the temperature, and the prior and likelihood functions of@@ -481,9 +531,12 @@     [ "MC3: Cycle of cold chain.",       coldMHGCycleSummary     ]-      ++ [ "MC3: Average acceptance rate across all chains: " <> BL.fromStrict (BC.toFixed 2 ar)-           | not $ isNaN ar-         ]+      ++ case mAr of+        Nothing -> []+        Just ar ->+          [ "MC3: Average acceptance rate across all chains: "+              <> BL.fromStrict (BC.toFixed 2 ar)+          ]       ++ [ "MC3: Reciprocal temperatures of the chains: " <> BL.intercalate ", " bsB <> ".",            "MC3: Summary of state swaps.",            "MC3: The swap period is " <> swapPeriodB <> ".",@@ -494,7 +547,7 @@       ++ [ summarizeProposal              (PName $ show i ++ " <-> " ++ show (i + 1))              (PDescription "Swap states between chains")-             (PWeight 1)+             (pWeight 1)              (Just $ bs U.! (i + 1))              PDimensionUnknown              (acceptanceRate i swapAcceptance)@@ -505,9 +558,12 @@     mhgs = mc3MHGChains a     coldMHGCycleSummary = aSummarizeCycle $ V.head mhgs     cs = V.map fromMHG mhgs-    as = V.map (acceptanceRates . acceptance) cs-    vAr = V.map (\m -> sum m / fromIntegral (length m)) as-    ar = V.sum vAr / fromIntegral (V.length vAr)+    -- 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+    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 (BL.fromStrict . BC.toFixed 2) $ U.toList bs     swapPeriod = fromSwapPeriod $ mc3SwapPeriod $ mc3Settings a

@@ -0,0 +1,283 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  Mcmc.Algorithm.MHG+-- Description :  Metropolis-Hastings-Green algorithm+-- Copyright   :  (c) Dominik Schrempf 2021+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue May  5 20:11:30 2020.+--+-- The Metropolis-Hastings-Green ('MHG') algorithm.+--+-- For example, see Geyer, C. J., Introduction to Markov chain Monte Carlo, In+-- Handbook of Markov Chain Monte Carlo (pp. 45) (2011). CRC press.+module Mcmc.Algorithm.MHG+  ( MHG (..),+    mhg,+    mhgSave,+    mhgLoad,+    MHGRatio,+    mhgAccept,+  )+where++import Codec.Compression.GZip+import Control.Monad+import Control.Monad.IO.Class+import Data.Aeson+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.Time+import Mcmc.Algorithm+import Mcmc.Chain.Chain+import Mcmc.Chain.Link+import Mcmc.Chain.Save+import Mcmc.Chain.Trace+import Mcmc.Monitor+import Mcmc.Posterior+import Mcmc.Proposal+import Mcmc.Settings+import Numeric.Log+import System.Random.MWC+import Text.Printf+import Prelude hiding (cycle)++-- | The MHG algorithm.+newtype MHG a = MHG {fromMHG :: Chain a}++instance ToJSON a => Algorithm (MHG a) where+  aName = const "Metropolis-Hastings-Green (MHG)"+  aIteration = iteration . fromMHG+  aIsInValidState = mhgIsInValidState+  aIterate = mhgIterate+  aAutoTune = mhgAutoTune+  aResetAcceptance = mhgResetAcceptance+  aSummarizeCycle = mhgSummarizeCycle+  aOpenMonitors = mhgOpenMonitors+  aExecuteMonitors = mhgExecuteMonitors+  aStdMonitorHeader = mhgStdMonitorHeader+  aCloseMonitors = mhgCloseMonitors+  aSave = mhgSave++-- NOTE: IO is required because the trace is mutable.++-- | Initialize an MHG algorithm.+mhg ::+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  TraceLength ->+  InitialState a ->+  GenIO ->+  IO (MHG a)+mhg pr lh cc mn trLen i0 g = do+  -- The trace is a mutable vector and the mutable state needs to be handled by+  -- a monad.+  tr <- replicateT traceLength l0+  return $ MHG $ Chain Nothing l0 0 tr ac g 0 pr lh cc mn+  where+    l0 = Link i0 (pr i0) (lh i0)+    ac = emptyA $ ccProposals cc+    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn+    minimumTraceLength = case trLen of+      TraceAuto -> 1+      TraceMinimum n -> n+    traceLength = maximum $ minimumTraceLength : batchMonitorSizes++mhgFn :: AnalysisName -> FilePath+mhgFn (AnalysisName nm) = nm ++ ".mcmc.mhg"++-- | Save an MHG algorithm.+mhgSave ::+  ToJSON a =>+  AnalysisName ->+  MHG a ->+  IO ()+mhgSave nm (MHG c) = do+  savedChain <- toSavedChain c+  BL.writeFile (mhgFn nm) $ compress $ encode savedChain++-- | Load an MHG algorithm.+--+-- See 'Mcmc.Mcmc.mcmcContinue'.+mhgLoad ::+  FromJSON a =>+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  AnalysisName ->+  IO (MHG a)+mhgLoad pr lh cc mn nm = do+  savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm)+  chain <- either error (fromSavedChain pr lh cc mn) savedChain+  return $ MHG chain++-- | MHG ratios are stored in log domain.+type MHGRatio = Log Double++-- The MHG ratio. This implementation has the following properties:+--+-- - The ratio is 'Infinity' if fX is zero. In this case, the proposal is always+--   accepted.+--+-- - The ratio 'NaN' if (fY or q or j) and fX are zero. In this case, the+--   proposal is always rejected.+--+-- This means that a chain in a state with posterior probability zero (fX=0) can+-- only move if a state with non-zero posterior probability is proposed.+-- Otherwise it is stuck. Therefore, I print a warning when the posterior+-- probability is zero in the beginning of the MCMC run. This is probably not+-- the best behavior, but see below.+--+-- There is a discrepancy between authors saying that one should (a) always+-- accept the new state when the current posterior is zero (Chapter 4 of [1],+-- [2]), or (b) almost surely reject the proposal when either fY or q are zero+-- (Chapter 1 of [1]).+--+-- Since I trust the author of Chapter 1 (Charles Geyer) I choose to follow+-- option (b). However, Option (a) is more user-friendly.+--+-- [1] Handbook of markov chain monte carlo (2011), CRC press.+--+-- [2] Dellaportas, P., & Roberts, G. O., An introduction to mcmc, Lecture Notes+-- in Statistics, (), 1–41 (2003).+-- http://dx.doi.org/10.1007/978-0-387-21811-3_1.+mhgRatio :: Posterior -> Posterior -> KernelRatio -> Jacobian -> MHGRatio+-- q = qYX / qXY * jXY; see 'ProposalSimple'.+-- j = Jacobian.+mhgRatio fX fY q j = fY / fX * q * j+{-# INLINE mhgRatio #-}++-- | Accept or reject a proposal with given MHG ratio?+mhgAccept :: MHGRatio -> GenIO -> IO Bool+mhgAccept r g+  | ln r >= 0.0 = return True+  | otherwise = do+    b <- uniform g+    return $ b < exp (ln r)++mhgPropose :: MHG a -> Proposal a -> IO (MHG a)+mhgPropose (MHG c) p = do+  -- 1. Sample new state.+  (!y, !q, !j) <- liftIO $ s x g+  -- 2. Calculate Metropolis-Hastings-Green ratio.+  let !pY = pF y+      !lY = lF y+      !r = mhgRatio (pX * lX) (pY * lY) q j+  -- 3. Accept or reject.+  -- if ln r >= 0.0+  --   then do+  --     let !ac' = pushA p True ac+  --     return $ MHG $ c {link = Link y pY lY, acceptance = ac'}+  --   else do+  --     b <- uniform g+  --     if b < exp (ln r)+  --       then do+  --         let !ac' = pushA p True ac+  --         return $ MHG $ c {link = Link y pY lY, acceptance = ac'}+  --       else do+  --         let !ac' = pushA p False ac+  --         return $ MHG $ c {acceptance = pushA p False ac'}+  accept <- mhgAccept r g+  if accept+    then do+      let !ac' = pushA p True ac+      return $ MHG $ c {link = Link y pY lY, acceptance = ac'}+    else do+      let !ac' = pushA p False ac+      return $ MHG $ c {acceptance = pushA p False ac'}+  where+    s = prSimple p+    (Link x pX lX) = link c+    pF = priorFunction c+    lF = likelihoodFunction c+    ac = acceptance c+    g = generator c++mhgPush :: MHG a -> IO (MHG a)+mhgPush (MHG c) = do+  t' <- pushT i t+  return $ MHG c {trace = t', iteration = succ n}+  where+    i = link c+    t = trace c+    n = iteration c++mhgIsInValidState :: MHG a -> Bool+mhgIsInValidState a = ((p * l) == 0) || (p * l == (0 / 0))+  where+    x = link $ fromMHG a+    p = prior x+    l = likelihood x++-- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG+-- algorithm is just inherently sequential. Parallelization can be achieved by+-- having parallel prior and/or likelihood functions, or by using algorithms+-- running parallel chains such as 'MC3'.+mhgIterate :: ParallelizationMode -> MHG a -> IO (MHG a)+mhgIterate _ a = do+  ps <- prepareProposals cc g+  a' <- foldM mhgPropose a ps+  mhgPush a'+  where+    c = fromMHG a+    cc = cycle c+    g = generator c++mhgAutoTune :: MHG a -> MHG a+mhgAutoTune (MHG c) = MHG $ c {cycle = autoTuneCycle ac cc}+  where+    ac = acceptance c+    cc = cycle c++mhgResetAcceptance :: MHG a -> MHG a+mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}+  where+    ac = acceptance c++mhgSummarizeCycle :: MHG a -> BL.ByteString+mhgSummarizeCycle (MHG c) = summarizeCycle ac cc+  where+    cc = cycle c+    ac = acceptance c++mhgOpenMonitors :: AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a)+mhgOpenMonitors nm em (MHG c) = do+  m' <- mOpen pre suf em m+  return $ MHG c {monitor = m'}+  where+    m = monitor c+    pre = fromAnalysisName nm+    suf = maybe "" (printf "%02d") $ chainId c++mhgExecuteMonitors ::+  Verbosity ->+  -- Starting time.+  UTCTime ->+  -- Total number of iterations.+  Int ->+  MHG a ->+  IO (Maybe BL.ByteString)+mhgExecuteMonitors vb t0 iTotal (MHG c) = mExec vb i i0 t0 tr iTotal m+  where+    i = iteration c+    i0 = start c+    tr = trace c+    m = monitor c++mhgStdMonitorHeader :: MHG a -> BL.ByteString+mhgStdMonitorHeader (MHG c) = msHeader (mStdOut $ monitor c)++mhgCloseMonitors :: MHG a -> IO (MHG a)+mhgCloseMonitors (MHG c) = do+  m' <- mClose m+  return $ MHG $ c {monitor = m'}+  where+    m = monitor c

@@ -1,254 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE OverloadedStrings #-}---- |--- Module      :  Mcmc.Algorithm.Metropolis--- Description :  Metropolis-Hastings-Green algorithm--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue May  5 20:11:30 2020.------ The Metropolis-Hastings-Green ('MHG') algorithm.------ For example, see Geyer, C. J., Introduction to Markov chain Monte Carlo, In--- Handbook of Markov Chain Monte Carlo (pp. 45) (2011). CRC press.-module Mcmc.Algorithm.Metropolis-  ( MHG (..),-    mhg,-    mhgSave,-    mhgLoad,-    mhgAccept,-  )-where--import Codec.Compression.GZip-import Control.Monad-import Control.Monad.IO.Class-import Data.Aeson-import qualified Data.ByteString.Lazy.Char8 as BL-import Data.Time-import Mcmc.Algorithm-import Mcmc.Chain.Chain-import Mcmc.Chain.Link-import Mcmc.Chain.Save-import Mcmc.Chain.Trace-import Mcmc.Monitor-import Mcmc.Proposal-import Mcmc.Settings-import Numeric.Log-import System.Random.MWC-import Text.Printf-import Prelude hiding (cycle)---- | The MHG algorithm.-newtype MHG a = MHG {fromMHG :: Chain a}--instance ToJSON a => Algorithm (MHG a) where-  aName = const "Metropolis-Hastings-Green (MHG)"-  aIteration = iteration . fromMHG-  aIterate = mhgIterate-  aAutoTune = mhgAutoTune-  aResetAcceptance = mhgResetAcceptance-  aSummarizeCycle = mhgSummarizeCycle-  aOpenMonitors = mhgOpenMonitors-  aExecuteMonitors = mhgExecuteMonitors-  aStdMonitorHeader = mhgStdMonitorHeader-  aCloseMonitors = mhgCloseMonitors-  aSave = mhgSave---- NOTE: IO is required because the trace is mutable.---- | Initialize an MHG algorithm.-mhg ::-  PriorFunction a ->-  LikelihoodFunction a ->-  Cycle a ->-  Monitor a ->-  -- | The initial state in the state space @a@.-  a ->-  -- | A source of randomness. For reproducible runs, make sure to use-  -- generators with the same seed.-  GenIO ->-  IO (MHG a)-mhg pr lh cc mn i0 g = do-  -- The trace is a mutable vector and the mutable state needs to be handled by-  -- a monad.-  tr <- replicateT traceLength l0-  return $ MHG $ Chain 0 l0 0 tr ac g 0 pr lh cc mn-  where-    l0 = Link i0 (pr i0) (lh i0)-    ac = emptyA $ ccProposals cc-    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn-    traceLength = maximum $ 1 : batchMonitorSizes--mhgFn :: AnalysisName -> FilePath-mhgFn (AnalysisName nm) = nm ++ ".mhg"---- | Save an MHG algorithm.-mhgSave ::-  ToJSON a =>-  AnalysisName ->-  MHG a ->-  IO ()-mhgSave nm (MHG c) = do-  savedChain <- toSavedChain c-  BL.writeFile (mhgFn nm) $ compress $ encode savedChain---- | Load an MHG algorithm.------ See 'Mcmc.Mcmc.mcmcContinue'.-mhgLoad ::-  FromJSON a =>-  PriorFunction a ->-  LikelihoodFunction a ->-  Cycle a ->-  Monitor a ->-  AnalysisName ->-  IO (MHG a)-mhgLoad pr lh cc mn nm = do-  savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm)-  chain <- either error (fromSavedChain pr lh cc mn) savedChain-  return $ MHG chain---- The MHG ratio.------ 'Infinity' if fX is zero. In this case, the proposal is always accepted.------ 'NaN' if (fY or q) and fX are zero. In this case, the proposal is always--- rejected.---- There is a discrepancy between authors saying that one should (a) always--- accept the new state when the current posterior is zero (Chapter 4 of the--- Handbook of Markov Chain Monte Carlo), or (b) almost surely reject the--- proposal when either fY or q are zero (Chapter 1). Since I trust the author--- of Chapter 1 (Charles Geyer) I choose to follow option (b).-mhgRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double--- q = qYX / qXY * jXY; see 'ProposalSimple'.--- j = Jacobian.-mhgRatio fX fY q j = fY / fX * q * j-{-# INLINE mhgRatio #-}---- | Accept or reject a proposal with given MHG ratio?-mhgAccept :: Log Double -> GenIO -> IO Bool-mhgAccept r g-  | ln r >= 0.0 = return True-  | otherwise = do-    b <- uniform g-    return $ b < exp (ln r)--mhgPropose :: MHG a -> Proposal a -> IO (MHG a)-mhgPropose (MHG c) p = do-  -- 1. Sample new state.-  (!y, !q, !j) <- liftIO $ s x g-  -- 2. Calculate Metropolis-Hastings-Green ratio.-  let !pY = pF y-      !lY = lF y-      !r = mhgRatio (pX * lX) (pY * lY) q j-  -- 3. Accept or reject.-  -- if ln r >= 0.0-  --   then do-  --     let !ac' = pushA p True ac-  --     return $ MHG $ c {link = Link y pY lY, acceptance = ac'}-  --   else do-  --     b <- uniform g-  --     if b < exp (ln r)-  --       then do-  --         let !ac' = pushA p True ac-  --         return $ MHG $ c {link = Link y pY lY, acceptance = ac'}-  --       else do-  --         let !ac' = pushA p False ac-  --         return $ MHG $ c {acceptance = pushA p False ac'}-  accept <- mhgAccept r g-  if accept-    then do-      let !ac' = pushA p True ac-      return $ MHG $ c {link = Link y pY lY, acceptance = ac'}-    else do-      let !ac' = pushA p False ac-      return $ MHG $ c {acceptance = pushA p False ac'}-  where-    s = pSimple p-    (Link x pX lX) = link c-    pF = priorFunction c-    lF = likelihoodFunction c-    ac = acceptance c-    g = generator c--mhgPush :: MHG a -> IO (MHG a)-mhgPush (MHG c) = do-  t' <- pushT i t-  return $ MHG c {trace = t', iteration = succ n}-  where-    i = link c-    t = trace c-    n = iteration c---- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG--- algorithm is just inherently sequential. Parallelization can be achieved by--- having parallel prior and/or likelihood functions, or by using algorithms--- running parallel chains such as 'MC3'.-mhgIterate :: ParallelizationMode -> MHG a -> IO (MHG a)-mhgIterate _ a = do-  ps <- orderProposals cc g-  a' <- foldM mhgPropose a ps-  mhgPush a'-  where-    c = fromMHG a-    cc = cycle c-    g = generator c--mhgAutoTune :: MHG a -> MHG a-mhgAutoTune (MHG c) = MHG $ c {cycle = autoTuneCycle ac cc}-  where-    ac = acceptance c-    cc = cycle c--mhgResetAcceptance :: MHG a -> MHG a-mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}-  where-    ac = acceptance c--mhgSummarizeCycle :: MHG a -> BL.ByteString-mhgSummarizeCycle (MHG c) = summarizeCycle ac cc-  where-    cc = cycle c-    ac = acceptance c--mhgOpenMonitors :: AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a)-mhgOpenMonitors nm em (MHG c) = do-  m' <- mOpen pre suf em m-  return $ MHG c {monitor = m'}-  where-    m = monitor c-    pre = fromAnalysisName nm-    suf = printf "%02d" $ chainId c--mhgExecuteMonitors ::-  Verbosity ->-  -- Starting time.-  UTCTime ->-  -- Total number of iterations.-  Int ->-  MHG a ->-  IO (Maybe BL.ByteString)-mhgExecuteMonitors vb t0 iTotal (MHG c) = mExec vb i i0 t0 tr iTotal m-  where-    i = iteration c-    i0 = start c-    tr = trace c-    m = monitor c--mhgStdMonitorHeader :: MHG a -> BL.ByteString-mhgStdMonitorHeader (MHG c) = msHeader (mStdOut $ monitor c)--mhgCloseMonitors :: MHG a -> IO (MHG a)-mhgCloseMonitors (MHG c) = do-  m' <- mClose m-  return $ MHG $ c {monitor = m'}-  where-    m = monitor c

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Chain.Chain -- Description :  Simple representation of a Markov chain--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -14,35 +14,29 @@     noPrior,     LikelihoodFunction,     noLikelihood,+    InitialState,     Chain (..),   ) where --- Note: It is not necessary to add another type @b@ to store supplementary+-- NOTE: It is not necessary to add another type @b@ to store supplementary -- information about the chain. The information can just be stored in @a@ -- equally well. +-- TODO: First class parameters. Make a type class for parameter types (name,+-- lens, proposals, monitors).+ import Mcmc.Chain.Link import Mcmc.Chain.Trace+import Mcmc.Likelihood import Mcmc.Monitor+import Mcmc.Prior import Mcmc.Proposal-import Numeric.Log import System.Random.MWC hiding (save) import Prelude hiding (cycle) --- | Prior function.-type PriorFunction a = a -> Log Double---- | Flat prior function. Useful for testing and debugging.-noPrior :: PriorFunction a-noPrior = const 1.0---- | Likelihood function.-type LikelihoodFunction a = a -> Log Double---- | Flat likelihood function. Useful for testing and debugging.-noLikelihood :: LikelihoodFunction a-noLikelihood = const 1.0+-- | Type synonym to indicate the initial state.+type InitialState a = a  -- | The chain contains all information to run an MCMC sampler. --@@ -60,8 +54,9 @@ -- The 'Mcmc.Environment.Environment' of the chain is not stored externally. data Chain a = Chain   { -- Variables; saved.+     -- | Chain index; useful if more chains are run.-    chainId :: Int,+    chainId :: Maybe Int,     -- | The current 'Link' of the chain combines the current state and the     -- current likelihood. The link is updated after a proposal has been     -- executed.

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Chain.Link -- Description :  The state combined with auxiliary variables--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -18,6 +18,8 @@  import Data.Aeson import Data.Aeson.Types+import Mcmc.Likelihood+import Mcmc.Prior import Numeric.Log  -- | Link of a Markov chain. For reasons of computational efficiency, each state@@ -25,10 +27,10 @@ data Link a = Link   { -- | The current state in the state space @a@.     state :: a,-    -- | The current prior.-    prior :: Log Double,-    -- | The current likelihood.-    likelihood :: Log Double+    -- | The current prior value.+    prior :: Prior,+    -- | The current likelihood value.+    likelihood :: Likelihood   }   deriving (Eq, Ord, Show, Read) 

@@ -1,13 +1,6 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-} --- |--- Module      :  Mcmc.Chain.Save--- Description :  Save and load a Markov chain--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later---- -- Maintainer  :  dominik.schrempf@gmail.com -- Stability   :  unstable -- Portability :  portable@@ -17,6 +10,12 @@ -- Save and load chains. It is easy to save and restore the current state and -- likelihood (or the trace), but it is not feasible to store all the proposals -- and so on, so they have to be provided again when continuing a run.++-- |+-- Module      :  Mcmc.Chain.Save+-- Description :  Save and load a Markov chain+-- Copyright   :  (c) Dominik Schrempf, 2021+-- License     :  GPL-3.0-or-later module Mcmc.Chain.Save   ( SavedChain (..),     toSavedChain,@@ -30,6 +29,7 @@ import Data.List hiding (cycle) import qualified Data.Map as M import Data.Maybe+import qualified Data.Stack.Circular as C import qualified Data.Vector as VB import qualified Data.Vector.Unboxed as VU import Data.Word@@ -40,20 +40,18 @@ import Mcmc.Monitor import Mcmc.Proposal import Prelude hiding (cycle)-import qualified Data.Stack.Circular as C  -- | Storable values of a Markov chain. -- -- See 'toSavedChain'. data SavedChain a = SavedChain-  {-    savedId :: Int,+  { savedId :: Maybe Int,     savedLink :: Link a,     savedIteration :: Int,     savedTrace :: C.Stack VB.Vector (Link a),     savedAcceptance :: Acceptance Int,     savedSeed :: VU.Vector Word32,-    savedTuningParameters :: [Maybe Double]+    savedTuningParameters :: [Maybe TuningParameter]   }   deriving (Eq, Read, Show) @@ -70,7 +68,7 @@   where     ps = ccProposals cc     ac' = transformKeysA ps [0 ..] ac-    ts = [fmap tParam mt | mt <- map pTuner ps]+    ts = [fmap tParam mt | mt <- map prTuner ps]  -- | Load a saved chain. --@@ -87,13 +85,19 @@   IO (Chain a) fromSavedChain pr lh cc mn (SavedChain ci it i tr ac' g' ts)   | pr (state it) /= prior it =-    error "fromSave: Provided prior function does not match the saved prior."+    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."   | 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'     getTuningF mt = case mt of

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Chain.Trace -- Description :  History of a Markov chain--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Environment -- Description :  Runtime environment--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -12,38 +12,75 @@ module Mcmc.Environment   ( Environment (..),     initializeEnvironment,+    closeEnvironment,   ) where +import Control.Concurrent.MVar+import Control.Monad import Data.Time.Clock+import Mcmc.Logger import Mcmc.Settings import System.IO  -- | The environment of an MCMC run.-data Environment = Environment-  { settings :: Settings,-    -- | We have to use 'Maybe' here, because we do not want to open any log-    -- file when being 'Quiet'.-    logHandle :: Maybe Handle,+data Environment s = Environment+  { settings :: s,+    -- | List will be empty if using 'Quiet'. If 'LogStdOutAndFile' is used+    -- 'logHandles' contains two handles to the standard output and the log+    -- file.+    logHandles :: [Handle],+    -- | MVar blocking output.+    outLock :: MVar (),     -- | Used to calculate the ETA.     startingTime :: UTCTime   }-  deriving (Eq, Show)+  deriving (Eq) +instance HasExecutionMode s => HasExecutionMode (Environment s) where+  getExecutionMode = getExecutionMode . settings++instance HasLock (Environment s) where+  getLock = outLock++instance HasLogHandles (Environment s) where+  getLogHandles = logHandles++instance HasStartingTime (Environment s) where+  getStartingTime = startingTime++instance HasLogMode s => HasLogMode (Environment s) where+  getLogMode = getLogMode . settings++instance HasVerbosity s => HasVerbosity (Environment s) where+  getVerbosity = getVerbosity . settings+ -- | Initialize the environment. -- -- Open log file, get current time. initializeEnvironment ::-  Settings ->-  IO Environment+  (HasAnalysisName s, HasExecutionMode s, HasLogMode s, HasVerbosity s) =>+  s ->+  IO (Environment s) initializeEnvironment s = do   t <- getCurrentTime-  mh <- case sVerbosity s of-    Quiet -> return Nothing-    _ -> do+  mh <- case (getLogMode s, getVerbosity s) of+    (_, Quiet) -> return []+    (LogStdOutAndFile, _) -> do       h <- openWithExecutionMode em fn-      return $ Just h-  return $ Environment s mh t+      return [stdout, h]+    (LogFileOnly, _) -> do+      h <- openWithExecutionMode em fn+      return [h]+    (LogStdOutOnly, _) -> return [stdout]+  lock <- newMVar ()+  return $ Environment s mh lock t   where-    fn = fromAnalysisName (sAnalysisName s) ++ ".log"-    em = sExecutionMode s+    fn = fromAnalysisName (getAnalysisName s) ++ ".mcmc.log"+    em = getExecutionMode s++-- | Close file handles.+closeEnvironment :: Environment s -> IO ()+closeEnvironment e = forM_ hs hClose+  where+    hs = filter (/= stdout) $ logHandles e

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Internal.ByteString -- Description :  ByteString tools--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Internal.Random -- Description :  Tools for random calculations--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Internal.Shuffle -- Description :  Shuffle a list--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -0,0 +1,29 @@+-- |+-- Module      :  Mcmc.Likelihood+-- Description :  Types and convenience functions for computing likelihoods+-- Copyright   :  (c) 2021 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Wed Mar  3 11:39:04 2021.+module Mcmc.Likelihood+  ( Likelihood,+    LikelihoodFunction,+    noLikelihood,+  )+where++import Numeric.Log++-- | Likelihood values are stored in log domain.+type Likelihood = Log Double++-- | Likelihood function.+type LikelihoodFunction a = a -> Log Double++-- | Flat likelihood function. Useful for testing and debugging.+noLikelihood :: LikelihoodFunction a+noLikelihood = const 1.0

@@ -0,0 +1,170 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TemplateHaskell #-}++-- |+-- Module      :  Mcmc.Logger+-- Description :  Minimal monad logger+-- Copyright   :  (c) Dominik Schrempf 2021+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Jan 12 09:03:04 2021.+module Mcmc.Logger+  ( LogMode (..),+    Verbosity (..),+    HasLock (..),+    HasLogHandles (..),+    HasStartingTime (..),+    HasLogMode (..),+    HasVerbosity (..),+    Logger,+    logOutB,+    logDebugB,+    logDebugS,+    logWarnB,+    logWarnS,+    logInfoB,+    logInfoS,+    logInfoHeader,+    logInfoStartingTime,+    logInfoEndTime,+  )+where++import Control.Concurrent.MVar+import Control.Monad+import Control.Monad.IO.Class+import Control.Monad.Trans.Reader+import Data.Aeson+import Data.Aeson.TH+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.Time.Clock+import Mcmc.Monitor.Time+import System.IO+import Paths_mcmc (version)+import Data.Version (showVersion)++-- TODO: Combine LogMode and Verbosity to:+--+-- data Verbosity = Quiet | Warn LogMode | Info LogMode | Debug LogMode++-- | Define where the log output should be directed to.+data LogMode = LogStdOutAndFile | LogStdOutOnly | LogFileOnly+  deriving (Eq, Read, Show)++$(deriveJSON defaultOptions ''LogMode)++-- | Not much to say here.+data Verbosity = Quiet | Warn | Info | Debug+  deriving (Eq, Ord, Read, Show)++$(deriveJSON defaultOptions ''Verbosity)++-- | Types with an output lock for concurrent output.+class HasLock e where+  getLock :: e -> MVar ()++-- | Types with logging information.+class HasLogHandles e where+  getLogHandles :: e -> [Handle]++-- | Types with starting time.+class HasStartingTime s where+  getStartingTime :: s -> UTCTime++-- | Types with a log mode.+class HasLogMode s where+  getLogMode :: s -> LogMode++-- | Types with verbosity.+class HasVerbosity s where+  getVerbosity :: s -> Verbosity++-- | Reader transformer used for logging to a file and to standard output.+type Logger e a = ReaderT e IO a++msgPrepare :: BL.ByteString -> BL.ByteString -> BL.ByteString+msgPrepare pref msg = BL.intercalate "\n" $ map (BL.append pref) $ BL.lines msg++-- Make sure that concurrent output is not scrambled.+atomicAction :: HasLock e => IO () -> Logger e ()+atomicAction a = do+  l <- reader getLock+  liftIO $ withMVar l (const a)++-- | Write to standard output and maybe to log file.+logOutB ::+  (HasLogHandles e, HasLock e) =>+  -- | Prefix.+  BL.ByteString ->+  -- | Message.+  BL.ByteString ->+  Logger e ()+logOutB pref msg = do+  hs <- reader getLogHandles+  mapM_ (atomicAction . (`BL.hPutStrLn` msg')) hs+  where+    msg' = msgPrepare pref msg++-- Perform debug action.+logDebugA :: (HasLock e, HasLogHandles e, HasVerbosity e) => Logger e () -> Logger e ()+logDebugA a = reader getVerbosity >>= \v -> when (v >= Debug) a++-- | Log debug message.+logDebugB :: (HasLock e, HasLogHandles e, HasVerbosity e) => BL.ByteString -> Logger e ()+logDebugB = logDebugA . logOutB "D: "++-- | Log debug message.+logDebugS :: (HasLock e, HasLogHandles e, HasVerbosity e) => String -> Logger e ()+logDebugS = logDebugB . BL.pack++-- Perform warning action.+logWarnA :: (HasLogHandles e, HasVerbosity e) => Logger e () -> Logger e ()+logWarnA a = reader getVerbosity >>= \v -> when (v >= Warn) a++-- | Log warning message.+logWarnB :: (HasLock e, HasLogHandles e, HasVerbosity e) => BL.ByteString -> Logger e ()+logWarnB = logWarnA . logOutB "W: "++-- | Log warning message.+logWarnS :: (HasLock e, HasLogHandles e, HasVerbosity e) => String -> Logger e ()+logWarnS = logWarnB . BL.pack++-- Perform info action.+logInfoA :: (HasLogHandles e, HasVerbosity e) => Logger e () -> Logger e ()+logInfoA a = reader getVerbosity >>= \v -> when (v >= Info) a++-- | Log info message.+logInfoB :: (HasLock e, HasLogHandles e, HasVerbosity e) => BL.ByteString -> Logger e ()+logInfoB = logInfoA . logOutB "   "++-- | Log info message.+logInfoS :: (HasLock e, HasLogHandles e, HasVerbosity e) => String -> Logger e ()+logInfoS = logInfoB . BL.pack++-- | Log info header.+logInfoHeader :: (HasLock e, HasLogHandles e, HasVerbosity e) => Logger e ()+logInfoHeader = do+  logInfoS (replicate 70 '-')+  logInfoS ("MCMC sampler; version " ++ showVersion version <> ".")+  logInfoS "Developed by: Dominik Schrempf."+  logInfoS "License: GPL-3.0-or-later."+  logInfoS (replicate 70 '-')++-- | Log starting time.+logInfoStartingTime :: (HasLock e, HasLogHandles e, HasStartingTime e, HasVerbosity e) => Logger e ()+logInfoStartingTime = do+  ti <- reader getStartingTime+  logInfoS $ "Starting time: " <> renderTime ti++-- | Log end time.+logInfoEndTime :: (HasLock e, HasLogHandles e, HasStartingTime e, HasVerbosity e) => Logger e ()+logInfoEndTime = do+  ti <- reader getStartingTime+  te <- liftIO getCurrentTime+  let dt = te `diffUTCTime` ti+  logInfoB $ "Wall clock run time: " <> renderDuration dt <> "."+  logInfoS $ "End time: " <> renderTime te

@@ -0,0 +1,404 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  Mcmc.MarginalLikelihood+-- Description :  Calculate the marginal likelihood+-- Copyright   :  (c) Dominik Schrempf 2021+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon Jan 11 16:34:18 2021.+module Mcmc.MarginalLikelihood+  ( MarginalLikelihood,+    NPoints (..),+    MLAlgorithm (..),+    MLSettings (..),+    marginalLikelihood,+  )+where++import Control.Monad+import Control.Monad.IO.Class+import qualified Control.Monad.Parallel as P+import Control.Monad.Trans.Reader+import Data.Aeson+import Data.List hiding (cycle)+import qualified Data.Map.Strict as M+import qualified Data.Vector as VB+import qualified Data.Vector.Unboxed as VU+import Mcmc.Algorithm.MHG+import Mcmc.Chain.Chain+import Mcmc.Chain.Link+import Mcmc.Chain.Trace+import Mcmc.Environment+import Mcmc.Internal.Random+import Mcmc.Likelihood+import Mcmc.Logger+import Mcmc.Mcmc+import Mcmc.Monitor+import Mcmc.Proposal+import Mcmc.Settings+import Numeric.Log hiding (sum)+import System.Directory+import System.Random.MWC+import Text.Printf+import Text.Show.Pretty+import Prelude hiding (cycle)++-- | Marginal likelihood values are stored in log domain.+type MarginalLikelihood = Log Double++-- Reciprocal temperature value traversed along the path integral.+type Point = Double++-- | The number of points used to approximate the path integral.+newtype NPoints = NPoints {fromNPoints :: Int}+  deriving (Eq, Read, Show)++-- | Algorithms to calculate the marginal likelihood.+data MLAlgorithm+  = -- | Use a classical path integral. Also known as thermodynamic integration.+    -- In particular, /Annealing-Melting Integration/ is used.+    --+    -- See Lartillot, N., & Philippe, H., Computing Bayes Factors Using+    -- Thermodynamic Integration, Systematic Biology, 55(2), 195–207 (2006).+    -- http://dx.doi.org/10.1080/10635150500433722+    ThermodynamicIntegration+  | -- | Use stepping stone sampling.+    --+    -- See Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M., Improving+    -- marginal likelihood estimation for Bayesian phylogenetic model selection,+    -- Systematic Biology, 60(2), 150–160 (2010).+    -- http://dx.doi.org/10.1093/sysbio/syq085+    --+    -- Or Fan, Y., Wu, R., Chen, M., Kuo, L., & Lewis, P. O., Choosing among+    -- partition models in bayesian phylogenetics, Molecular Biology and+    -- Evolution, 28(1), 523–532 (2010). http://dx.doi.org/10.1093/molbev/msq224+    SteppingStoneSampling+  deriving (Eq, Read, Show)++-- | Settings of the marginal likelihood estimation.+data MLSettings = MLSettings+  { mlAnalysisName :: AnalysisName,+    mlAlgorithm :: MLAlgorithm,+    mlNPoints :: NPoints,+    -- | Initial burn in at the starting point of the path.+    mlInitialBurnIn :: BurnInSpecification,+    -- | Repetitive burn in at each point on the path.+    mlPointBurnIn :: BurnInSpecification,+    -- | The number of iterations performed at each point.+    mlIterations :: Iterations,+    mlExecutionMode :: ExecutionMode,+    mlLogMode :: LogMode,+    mlVerbosity :: Verbosity+  }+  deriving (Eq, Read, Show)++instance HasAnalysisName MLSettings where+  getAnalysisName = mlAnalysisName++instance HasExecutionMode MLSettings where+  getExecutionMode = mlExecutionMode++instance HasLogMode MLSettings where+  getLogMode = mlLogMode++instance HasVerbosity MLSettings where+  getVerbosity = mlVerbosity++type ML a = ReaderT (Environment MLSettings) IO a++-- See 'getPoints'. Alpha=0.3 is the standard choice.+alpha :: Double+alpha = 0.3++-- Distribute the points according to a skewed beta distribution with given+-- 'alpha' value. If alpha is below 1.0, more points at lower values, which is+-- desired. It is inconvenient that the reciprocal temperatures are denoted as+-- beta, and we also use the beta distribution :). Don't mix them up!+--+-- See discussion in Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M.,+-- Improving marginal likelihood estimation for bayesian phylogenetic model+-- selection, Systematic Biology, 60(2), 150–160 (2010).+-- http://dx.doi.org/10.1093/sysbio/syq085+--+-- Or Figure 1 in Höhna, S., Landis, M. J., & Huelsenbeck, J. P., Parallel power+-- posterior analyses for fast computation of marginal likelihoods in+-- phylogenetics (2017). http://dx.doi.org/10.1101/104422+getPoints :: NPoints -> [Point]+getPoints x = [f i ** (1.0 / alpha) | i <- [0 .. k1]]+  where+    k = fromNPoints x+    k1 = pred k+    f j = fromIntegral j / fromIntegral k1++sampleAtPoint ::+  ToJSON a =>+  Point ->+  Settings ->+  LikelihoodFunction a ->+  MHG a ->+  ML (MHG a)+sampleAtPoint x ss lhf a = do+  a'' <- liftIO $ mcmc ss' a'+  let ch'' = fromMHG a''+      ac = acceptance ch''+      mAr = sequence $ acceptanceRates ac+  logDebugB "sampleAtPoint: Summarize cycle."+  logDebugB $ summarizeCycle ac $ cycle ch''+  case mAr of+    Nothing -> logWarnB "Some acceptance rates are unavailable. The tuning period may be too small."+    Just ar -> do+      unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1."+      unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9."+  return a''+  where+    -- For debugging set a proper analysis name.+    nm = sAnalysisName ss+    getName :: Point -> AnalysisName+    getName y = nm <> AnalysisName ("/" <> printf "point%.8f" y)+    ss' = ss {sAnalysisName = getName x}+    -- Amend the likelihood function. Don't calculate the likelihood when the+    -- point is 0.0.+    lhf' = if x == 0.0 then const 1.0 else (** Exp (log x)) . lhf+    -- Amend the MHG algorithm.+    ch = fromMHG a+    l = link ch+    ch' =+      ch+        { -- Important: Update the likelihood using the new likelihood function.+          link = l {likelihood = lhf' $ state l},+          iteration = 0,+          start = 0,+          likelihoodFunction = lhf'+        }+    a' = MHG ch'++traversePoints ::+  ToJSON a =>+  -- Current point.+  Int ->+  NPoints ->+  [Point] ->+  Settings ->+  LikelihoodFunction a ->+  MHG a ->+  -- For each point a vector of obtained likelihoods stored in the log domain.+  ML [VU.Vector Likelihood]+traversePoints _ _ [] _ _ _ = return []+traversePoints i k (b : bs) ss lhf a = do+  logInfoS $ "Point " <> show i <> " of " <> show k' <> ": " <> show b <> "."+  a' <- sampleAtPoint b ss lhf a+  -- Get the links samples at this point.+  ls <- liftIO $ takeT n $ trace $ fromMHG a'+  -- Extract the likelihoods.+  --+  -- NOTE: This could be sped up by mapping (** -b) on the power likelihoods.+  --+  -- NOTE: This bang is an important one, because if the lhs are not strictly+  -- calculated here, the complete MCMC runs are dragged along before doing so+  -- resulting in a severe memory leak.+  let !lhs = VU.convert $ VB.map (lhf . state) ls+  -- Sample the other points.+  lhss <- traversePoints (i + 1) k bs ss lhf a'+  return $ lhs : lhss+  where+    n = fromIterations $ sIterations ss+    (NPoints k') = k++mlRun ::+  ToJSON a =>+  NPoints ->+  [Point] ->+  ExecutionMode ->+  Verbosity ->+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  a ->+  GenIO ->+  -- For each point a vector of likelihoods stored in log domain.+  ML [VU.Vector Likelihood]+mlRun k xs em vb prf lhf cc mn i0 g = do+  logDebugB "mlRun: Begin."+  s <- reader settings+  let nm = mlAnalysisName s+      is = mlIterations s+      biI = mlInitialBurnIn s+      biP = mlPointBurnIn s+      -- Only log sub MCMC samplers when debugging.+      vb' = case vb of+        Debug -> Debug+        _ -> Quiet+      ssI = Settings nm biI (Iterations 0) em Sequential NoSave LogFileOnly vb'+      ssP = Settings nm biP is em Sequential NoSave LogFileOnly vb'+      trLen = TraceMinimum $ fromIterations is+  logDebugB "mlRun: Initialize MHG algorithm."+  a0 <- liftIO $ mhg prf lhf cc mn trLen i0 g+  logDebugS $ "mlRun: Perform initial burn in at first point " <> show x0 <> "."+  a1 <- sampleAtPoint x0 ssI lhf a0+  logDebugB "mlRun: Traverse points."+  traversePoints 1 k xs ssP lhf a1+  where+    x0 = head xs++-- Use lists since the number of points is expected to be low.+integrateSimpsonTriangle ::+  -- X values.+  [Point] ->+  -- Y values.+  [Double] ->+  -- Integral.+  Double+integrateSimpsonTriangle xs ys = 0.5 * go xs ys+  where+    go (p0 : p1 : ps) (z0 : z1 : zs) = (z0 + z1) * (p1 - p0) + go (p1 : ps) (z1 : zs)+    go _ _ = 0++tiWrapper ::+  ToJSON a =>+  MLSettings ->+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  a ->+  GenIO ->+  ML MarginalLikelihood+tiWrapper s prf lhf cc mn i0 g = do+  logInfoB "Path integral (thermodynamic integration)."+  [g0, g1] <- splitGen 2 g++  -- Parallel execution of both path integrals.+  [lhssForward, lhssBackward] <-+    P.sequence+      [ mlRun k bsForward em vb prf lhf cc mn i0 g0,+        mlRun k bsBackward em vb prf lhf cc mn i0 g1+      ]+  logInfoEndTime++  logDebugB "tiWrapper: Calculate mean log likelihoods."+  -- It is important to average across the log likelihoods here (and not the+  -- likelihoods). I am not exactly sure why this is.+  let getMeanLogLhs = map (\x -> VU.sum (VU.map ln x) / fromIntegral (VU.length x))+      mlForward = integrateSimpsonTriangle bsForward (getMeanLogLhs lhssForward)+      mlBackward = negate $ integrateSimpsonTriangle bsBackward (getMeanLogLhs lhssBackward)+  logDebugS $ "tiWrapper: Marginal log likelihood of forward integral: " ++ show mlForward+  logDebugS $ "tiWrapper: Marginal log likelihood of backward integral: " ++ show mlBackward+  let mean = 0.5 * (mlForward + mlBackward)+  logDebugS $ "tiWrapper: The mean is: " ++ show mean+  return $ Exp mean+  where+    k = mlNPoints s+    bsForward = getPoints k+    bsBackward = reverse bsForward+    em = mlExecutionMode s+    vb = mlVerbosity s++-- Helper function to exponentiate log domain values with a double value.+pow' :: Log Double -> Double -> Log Double+pow' x p = Exp $ ln x * p++-- See Xie2010 p. 153, bottom left.+sssCalculateMarginalLikelihood :: [Point] -> [VU.Vector Likelihood] -> MarginalLikelihood+sssCalculateMarginalLikelihood xs lhss = product $ zipWith3 f xs (tail xs) lhss+  where+    f :: Point -> Point -> VU.Vector Likelihood -> MarginalLikelihood+    -- f beta_{k-1} beta_k lhs_{k-1}+    f bkm1 bk lhs = n1 * VU.sum lhsPowered+      where+        n1 = recip $ fromIntegral $ VU.length lhs+        dbeta = bk - bkm1+        lhsPowered = VU.map (`pow'` dbeta) lhs++-- -- Numerical stability by factoring out lhMax. But no observed+-- -- improvement towards the standard version.+--+-- f bkm1 bk lhs = n1 * pow' lhMax dbeta * VU.sum lhsNormedPowered+--   where n1 = recip $ fromIntegral $ VU.length lhs+--         lhMax = VU.maximum lhs+--         dbeta = bk - bkm1+--         lhsNormed = VU.map (/lhMax) lhs+--         lhsNormedPowered = VU.map (`pow'` dbeta) lhsNormed++-- -- Computation of the log of the marginal likelihood. According to the paper,+-- -- this estimator is biased and I did not observe any improvements compared+-- -- to the direct estimator implemented above.+--+-- -- See Xie2010 p. 153, top right.+-- sssCalculateMarginalLikelihood' :: [Point] -> [VU.Vector Likelihood] -> MarginalLikelihood+-- sssCalculateMarginalLikelihood' xs lhss = Exp $ sum $ zipWith3 f xs (tail xs) lhss+--   where f :: Point -> Point -> VU.Vector Likelihood -> Double+--         -- f beta_{k-1} beta_k lhs_{k-1}+--         f bkm1 bk lhs = dbeta * llhMax + log (n1 * VU.sum lhsNormedPowered)+--           where dbeta = bk - bkm1+--                 llhMax = ln $ VU.maximum lhs+--                 n1 = recip $ fromIntegral $ VU.length lhs+--                 llhs = VU.map ln lhs+--                 llhsNormed = VU.map (\x -> x - llhMax) llhs+--                 lhsNormedPowered = VU.map (\x -> exp $ dbeta * x) llhsNormed+sssWrapper ::+  ToJSON a =>+  MLSettings ->+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  a ->+  GenIO ->+  ML MarginalLikelihood+sssWrapper s prf lhf cc mn i0 g = do+  logInfoB "Stepping stone sampling."+  logLhss <- mlRun k bsForward' em vb prf lhf cc mn i0 g+  logInfoB "The last point does not need to be sampled with stepping stone sampling."+  logDebugB "sssWrapper: Calculate marginal likelihood."+  return $ sssCalculateMarginalLikelihood bsForward logLhss+  where+    k = mlNPoints s+    bsForward = getPoints k+    bsForward' = init bsForward+    em = mlExecutionMode s+    vb = mlVerbosity s++-- | Estimate the marginal likelihood.+marginalLikelihood ::+  ToJSON a =>+  MLSettings ->+  PriorFunction a ->+  LikelihoodFunction a ->+  Cycle a ->+  Monitor a ->+  InitialState a ->+  -- | A source of randomness. For reproducible runs, make sure to use+  -- generators with the same seed.+  GenIO ->+  IO MarginalLikelihood+marginalLikelihood s prf lhf cc mn i0 g = do+  -- Initialize.+  e <- initializeEnvironment s++  when (mlVerbosity s == Debug) $ do+    let n = fromAnalysisName $ mlAnalysisName s+    createDirectoryIfMissing True n++  -- Run.+  runReaderT+    ( do+        logInfoStartingTime+        logInfoB "Estimate marginal likelihood."+        logDebugB "marginalLikelihood: The marginal likelihood settings are:"+        logDebugS $ ppShow s+        val <- case mlAlgorithm s of+          ThermodynamicIntegration -> tiWrapper s prf lhf cc mn i0 g+          SteppingStoneSampling -> sssWrapper s prf lhf cc mn i0 g+        logInfoS $ "Marginal log likelihood: " ++ show (ln val)+        -- TODO: Simulation variance.+        logInfoS "The simulation variance is not yet available."+        return val+    )+    e

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Mcmc -- Description :  Framework for running Markov chain Monte Carlo samplers--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -23,14 +23,10 @@  import Control.Monad import Control.Monad.IO.Class--- import Control.Monad.Trans.RWS.CPS import Control.Monad.Trans.Reader-import qualified Data.ByteString.Lazy.Char8 as BL-import Data.Maybe-import Data.Time.Clock import Mcmc.Algorithm import Mcmc.Environment-import Mcmc.Monitor.Time+import Mcmc.Logger import Mcmc.Settings import System.IO import Text.Show.Pretty@@ -38,79 +34,24 @@  -- The MCMC algorithm has read access to an environment and uses an algorithm -- transforming the state @a@.-type MCMC a = ReaderT Environment IO a--msgPrepare :: BL.ByteString -> BL.ByteString -> BL.ByteString-msgPrepare pref msg = BL.intercalate "\n" $ map (BL.append pref) $ BL.lines msg---- Write to standard output and log file.-mcmcOutB :: BL.ByteString -> BL.ByteString -> MCMC ()-mcmcOutB pref msg = do-  h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> reader logHandle-  liftIO $ BL.putStrLn msg' >> BL.hPutStrLn h msg'-  where-    msg' = msgPrepare pref msg---- -- Perform warning action.--- mcmcWarnA :: MCMC a () -> MCMC a ()--- mcmcWarnA a = reader (verbosity . settings) >>= \v -> when (v >= Warn) a---- -- Print warning message.--- mcmcWarnB :: BL.ByteString -> MCMC a ()--- mcmcWarnB = mcmcWarnA . mcmcOutB . msgPrepare 'W'---- -- Print warning message.--- mcmcWarnS :: String -> MCMC a ()--- mcmcWarnS = mcmcWarnB . BL.pack---- Perform info action.-mcmcInfoA :: MCMC () -> MCMC ()-mcmcInfoA a = reader (sVerbosity . settings) >>= \v -> when (v >= Info) a---- Print info message.-mcmcInfoB :: BL.ByteString -> MCMC ()-mcmcInfoB = mcmcInfoA . mcmcOutB "I: "---- Print info message.-mcmcInfoS :: String -> MCMC ()-mcmcInfoS = mcmcInfoB . BL.pack---- Perform debug action.-mcmcDebugA :: MCMC () -> MCMC ()-mcmcDebugA a = reader (sVerbosity . settings) >>= \v -> when (v == Debug) a---- Print debug message.-mcmcDebugB :: BL.ByteString -> MCMC ()-mcmcDebugB = mcmcDebugA . mcmcOutB "D: "---- Print debug message.-mcmcDebugS :: String -> MCMC ()-mcmcDebugS = mcmcDebugB . BL.pack--mcmcReportTime :: MCMC ()-mcmcReportTime = do-  mcmcDebugB "Report time."-  ti <- reader startingTime-  mcmcInfoS $ "Starting time of MCMC sampler: " <> renderTime ti+type MCMC a = ReaderT (Environment Settings) IO a  mcmcExecute :: Algorithm a => a -> MCMC a mcmcExecute a = do-  mcmcDebugB "Executing MCMC run."+  logDebugB "Executing MCMC run."   s <- reader settings   a' <- case sExecutionMode s of     Fail -> mcmcNewRun a     Overwrite -> mcmcNewRun a     Continue -> mcmcContinueRun a-  mcmcDebugB "Executed MCMC run."+  logDebugB "Executed MCMC run."   return a' --- Reset acceptance counts. mcmcResetAcceptance :: Algorithm a => a -> MCMC a mcmcResetAcceptance a = do-  mcmcDebugB "Reset acceptance rates."+  logDebugB "Reset acceptance rates."   return $ aResetAcceptance a --- Execute the monitors of the chain. mcmcExecuteMonitors :: Algorithm a => a -> MCMC () mcmcExecuteMonitors a = do   e <- ask@@ -119,7 +60,7 @@       t0 = startingTime e       iTotal = burnInIterations (sBurnIn s) + fromIterations (sIterations s)   mStdLog <- liftIO (aExecuteMonitors vb t0 iTotal a)-  forM_ mStdLog (mcmcOutB "   ")+  forM_ mStdLog (logOutB "   ")  mcmcIterate :: Algorithm a => Int -> a -> MCMC a mcmcIterate n a@@ -134,30 +75,30 @@ mcmcNewRun :: Algorithm a => a -> MCMC a mcmcNewRun a = do   s <- reader settings-  mcmcInfoB "Start new MCMC sampler."-  mcmcInfoB "Initial state."-  mcmcInfoB $ aStdMonitorHeader a+  logInfoB "Start new MCMC sampler."+  logInfoB "Initial state."+  logInfoB $ aStdMonitorHeader a   mcmcExecuteMonitors a-  mcmcInfoB $ aSummarizeCycle a+  when (aIsInValidState a) (logWarnB "The initial state is invalid!")+  logInfoB $ aSummarizeCycle a   a' <- mcmcBurnIn a-  a'' <- mcmcResetAcceptance a'   let i = fromIterations $ sIterations s-  mcmcInfoS $ "Run chain for " ++ show i ++ " iterations."-  mcmcInfoB $ aStdMonitorHeader a''-  mcmcIterate i a''+  logInfoS $ "Run chain for " ++ show i ++ " iterations."+  logInfoB $ aStdMonitorHeader a'+  mcmcIterate i a'  mcmcContinueRun :: Algorithm a => a -> MCMC a mcmcContinueRun a = do   s <- reader settings   let iTotal = fromIterations (sIterations s) + burnInIterations (sBurnIn s)-  mcmcInfoB "Continuation of MCMC sampler."+  logInfoB "Continuation of MCMC sampler."   let iCurrent = aIteration a-  mcmcInfoS $ "Current iteration: " ++ show iCurrent ++ "."-  mcmcInfoS $ "Total iterations: " ++ show iTotal ++ "."+  logInfoS $ "Current iteration: " ++ show iCurrent ++ "."+  logInfoS $ "Total iterations: " ++ show iTotal ++ "."   let di = iTotal - iCurrent-  mcmcInfoB $ aSummarizeCycle a-  mcmcInfoS $ "Run chain for " ++ show di ++ " iterations."-  mcmcInfoB $ aStdMonitorHeader a+  logInfoB $ aSummarizeCycle a+  logInfoS $ "Run chain for " ++ show di ++ " iterations."+  logInfoB $ aStdMonitorHeader a   mcmcIterate di a  mcmcBurnIn :: Algorithm a => a -> MCMC a@@ -165,53 +106,67 @@   s <- reader settings   case sBurnIn s of     NoBurnIn -> do-      mcmcInfoS "No burn in."+      logInfoS "No burn in."       return a     BurnInWithoutAutoTuning n -> do-      mcmcInfoS $ "Burn in for " <> show n <> " iterations."-      mcmcInfoS "Auto tuning is disabled."-      mcmcInfoB $ aStdMonitorHeader a+      logInfoS $ "Burn in for " <> show n <> " iterations."+      logInfoS "Auto tuning is disabled."+      logInfoB $ aStdMonitorHeader a       a' <- mcmcIterate n a-      mcmcInfoB $ aSummarizeCycle a'-      mcmcInfoB "Burn in finished."-      return a'+      logInfoB $ aSummarizeCycle a'+      a'' <- mcmcResetAcceptance a'+      logInfoB "Burn in finished."+      return a''     BurnInWithAutoTuning n t -> do-      mcmcInfoS $ "Burn in for " ++ show n ++ " iterations."-      mcmcInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "."-      mcmcInfoB $ aStdMonitorHeader a-      a' <- mcmcBurnInWithAutoTuning n t a-      mcmcInfoB "Burn in finished."+      logInfoS $ "Burn in for " ++ show n ++ " iterations."+      logInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "."+      logInfoB $ aStdMonitorHeader a+      let (m, r) = n `divMod` t+          -- Don't add if 0. Because then we auto tune without acceptance counts+          -- and get NaNs.+          xs = replicate m t <> [r | r > 0]+      a' <- mcmcBurnInWithAutoTuning xs a+      logInfoB "Burn in finished."       return a'+    BurnInWithCustomAutoTuning xs -> do+      logInfoS $ "Burn in for " ++ show (sum xs) ++ " iterations."+      logInfoS $ "Custom auto tuning is enabled with periods " ++ show xs ++ "."+      logInfoB $ aStdMonitorHeader a+      a' <- mcmcBurnInWithAutoTuning xs a+      logInfoB "Burn in finished."+      return a'  -- Auto tune the proposals. mcmcAutotune :: Algorithm a => a -> MCMC a mcmcAutotune a = do-  mcmcDebugB "Auto tune."+  logDebugB "Auto tune."   return $ aAutoTune a -mcmcBurnInWithAutoTuning :: Algorithm a => Int -> Int -> a -> MCMC a-mcmcBurnInWithAutoTuning b t a-  | b > t = do-    a' <- mcmcResetAcceptance a-    a'' <- mcmcIterate t a'-    mcmcDebugB $ aSummarizeCycle a''-    a''' <- mcmcAutotune a''-    mcmcDebugB $ aStdMonitorHeader a''-    mcmcBurnInWithAutoTuning (b - t) t a'''-  | otherwise = do-    a' <- mcmcResetAcceptance a-    a'' <- mcmcIterate b a'-    mcmcInfoB $ aSummarizeCycle a''-    mcmcInfoS $ "Acceptance rates calculated over the last " <> show b <> " iterations."-    return a''+mcmcBurnInWithAutoTuning :: Algorithm a => [Int] -> a -> MCMC a+mcmcBurnInWithAutoTuning [] _ = error "mcmcBurnInWithAutoTuning: Empty lisst."+mcmcBurnInWithAutoTuning [x] a = do+  -- Last round.+  a' <- mcmcIterate x a+  a'' <- mcmcAutotune a'+  logInfoB $ aSummarizeCycle a''+  logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations."+  mcmcResetAcceptance a''+mcmcBurnInWithAutoTuning (x:xs) a = do+  a' <- mcmcIterate x a+  a'' <- mcmcAutotune a'+  logDebugB $ aSummarizeCycle a''+  logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations."+  logDebugB $ aStdMonitorHeader a''+  a''' <- mcmcResetAcceptance a''+  mcmcBurnInWithAutoTuning xs a'''  mcmcInitialize :: Algorithm a => a -> MCMC a mcmcInitialize a = do-  mcmcInfoS $ aName a ++ " algorithm."+  logInfoS $ aName a ++ " algorithm."   s <- settings <$> ask-  mcmcDebugB "Opening monitors."+  logDebugB "Opening monitors."   a' <- liftIO $ aOpenMonitors (sAnalysisName s) (sExecutionMode s) a-  mcmcDebugB "Monitors opened."+  logDebugB "Monitors opened."   return a'  -- Save the MCMC run.@@ -219,42 +174,41 @@ mcmcSave a = do   s <- reader settings   case sSaveMode s of-    NoSave -> mcmcInfoB "Do not save the MCMC analysis."+    NoSave -> logInfoB "Do not save the MCMC analysis."     Save -> do-      mcmcInfoB "Save settings."+      logInfoB "Save settings."       liftIO $ settingsSave s       let nm = sAnalysisName s-      mcmcInfoB "Save compressed MCMC analysis."-      mcmcInfoB "For long traces, or complex objects, this may take a while."+      logInfoB "Save compressed MCMC analysis."+      logInfoB "For long traces, or complex objects, this may take a while."       liftIO $ aSave nm a-      mcmcInfoB "Markov chain saved."+      logInfoB "Markov chain saved."  -- Report and finish up. mcmcClose :: Algorithm a => a -> MCMC a mcmcClose a = do-  mcmcDebugB "Closing MCMC run."-  mcmcInfoB $ aSummarizeCycle a-  mcmcInfoS $ aName a ++ " algorithm finished."+  logDebugB "Closing MCMC run."+  logInfoB $ aSummarizeCycle a+  logInfoS $ aName a ++ " algorithm finished."   mcmcSave a-  ti <- reader startingTime-  te <- liftIO getCurrentTime-  let dt = te `diffUTCTime` ti-  mcmcInfoB $ "Wall clock run time: " <> renderDuration dt <> "."-  mcmcInfoS $ "End time: " <> renderTime te+  logInfoEndTime   a' <- liftIO $ aCloseMonitors a-  h <- reader logHandle-  liftIO $ forM_ h hClose+  e <- ask+  liftIO $ closeEnvironment e   return a'  -- Initialize the run, execute the run, and close the run. mcmcRun :: Algorithm a => a -> MCMC a mcmcRun a = do-  mcmcDebugB "The settings are:"-  reader settings >>= mcmcDebugS . ppShow+  logInfoHeader +  -- Debug settings.+  logDebugB "The MCMC settings are:"+  reader settings >>= logDebugS . ppShow+   -- Initialize.   a' <- mcmcInitialize a-  mcmcReportTime+  logInfoStartingTime    -- Execute.   a'' <- mcmcExecute a'@@ -274,15 +228,15 @@ -- Currently, it is only possible to continue MCMC algorithms that have -- completed successfully. This restriction is necessary, because for parallel -- chains, it is hardly possible to ensure all chains are synchronized when the--- process is killed.+-- process is killed or fails. -- -- See: ----- - 'Mcmc.Algorithm.Metropolis.mhgLoad'+-- - 'Mcmc.Algorithm.MHG.mhgLoad' -- -- - 'Mcmc.Algorithm.MC3.mc3Load'-mcmcContinue :: Algorithm a => Int -> Settings -> a -> IO a+mcmcContinue :: Algorithm a => Iterations -> Settings -> a -> IO a mcmcContinue dn s = mcmc s'   where-    n' = Iterations $ fromIterations (sIterations s) + dn+    n' = Iterations $ fromIterations (sIterations s) + fromIterations dn     s' = s {sIterations = n', sExecutionMode = Continue}

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor -- Description :  Monitor a Markov chain--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -14,11 +14,14 @@ module Mcmc.Monitor   ( -- * Create monitors     Monitor (..),+    Period,+    simpleMonitor,     MonitorStdOut,     monitorStdOut,     msHeader,     MonitorFile,     monitorFile,+    BatchSize,     MonitorBatch,     monitorBatch,     getMonitorBatchSize,@@ -64,21 +67,31 @@     mBatches :: [MonitorBatch a]   } +-- | Monitor period.+type Period = Int++-- | Do not monitor parameters.+--+-- Monitor prior and likelihood with given period.+simpleMonitor :: Period -> Monitor a+simpleMonitor p+  | p < 1 = error "simpleMonitor: Monitor period must be 1 or larger."+  | otherwise =+    Monitor (MonitorStdOut [] p) [] []+ -- | Monitor to standard output; constructed with 'monitorStdOut'. data MonitorStdOut a = MonitorStdOut   { msParams :: [MonitorParameter a],-    msPeriod :: Int+    msPeriod :: Period   }  -- | Monitor to standard output. monitorStdOut ::-  -- | Instructions about which parameters to log.   [MonitorParameter a] ->-  -- | Logging period.-  Int ->+  Period ->   MonitorStdOut a monitorStdOut ps p-  | p < 1 = error "monitorStdOut: Monitor period has to be 1 or larger."+  | p < 1 = error "monitorStdOut: Monitor period must be 1 or larger."   | otherwise = MonitorStdOut ps p  msIWidth :: Int@@ -151,20 +164,18 @@   { mfName :: String,     mfHandle :: Maybe Handle,     mfParams :: [MonitorParameter a],-    mfPeriod :: Int+    mfPeriod :: Period   }  -- | Monitor parameters to a file. monitorFile ::   -- | Name; used as part of the file name.   String ->-  -- | Instructions about which parameters to log.   [MonitorParameter a] ->-  -- | Logging period.-  Int ->+  Period ->   MonitorFile a monitorFile n ps p-  | p < 1 = error "monitorFile: Monitor period has to be 1 or larger."+  | p < 1 = error "monitorFile: Monitor period must be 1 or larger."   | otherwise = MonitorFile n Nothing ps p  mfRenderRow :: [BL.ByteString] -> BL.ByteString@@ -216,6 +227,9 @@   Just h -> hClose h   Nothing -> error $ "mfClose: File was not opened for monitor " <> mfName m <> "." +-- | Batch size.+type BatchSize = Int+ -- | Batch monitor to a file. -- -- Calculate summary statistics over the last given number of iterations (batch@@ -224,26 +238,24 @@   { mbName :: String,     mbHandle :: Maybe Handle,     mbParams :: [MonitorParameterBatch a],-    mbSize :: Int+    mbSize :: BatchSize   }  -- | Batch monitor parameters to a file, see 'MonitorBatch'. monitorBatch ::   -- | Name; used as part of the file name.   String ->-  -- | Instructions about how to calculate the summary statistics.   [MonitorParameterBatch a] ->-  -- | Batch size.-  Int ->+  BatchSize ->   MonitorBatch a-monitorBatch n ps p-  | p < 2 = error "monitorBatch: Batch size has to be 2 or larger."-  | otherwise = MonitorBatch n Nothing ps p+monitorBatch n ps b+  | b < 2 = error "monitorBatch: Batch size must be 2 or larger."+  | otherwise = MonitorBatch n Nothing ps b  -- | Batch monitor size. -- -- Useful to determine the trace length.-getMonitorBatchSize :: MonitorBatch a -> Int+getMonitorBatchSize :: MonitorBatch a -> BatchSize getMonitorBatchSize = mbSize  mbOpen :: String -> String -> ExecutionMode -> MonitorBatch a -> IO (MonitorBatch a)@@ -322,7 +334,6 @@ -- | Execute monitors; print status information to files and return text to be -- printed to standard output and log file. mExec ::-  -- | Verbosity   Verbosity ->   -- | Iteration.   Int ->@@ -330,11 +341,9 @@   Int ->   -- | Starting time.   UTCTime ->-  -- | Trace of Markov chain.   Trace a ->   -- | Total number of iterations; to calculate ETA.   Int ->-  -- | The monitor.   Monitor a ->   IO (Maybe BL.ByteString) mExec v i ss st xs j (Monitor s fs bs) = do

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Monitor.Log -- Description :  Monitor logarithmic values--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor.Parameter -- Description :  Monitor parameters--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor.ParameterBatch -- Description :  Batch monitor parameters--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor.Time -- Description :  Print time related values--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -0,0 +1,20 @@+-- |+-- Module      :  Mcmc.Posterior+-- Description :  Types for posterior values and functions+-- Copyright   :  (c) 2021 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Fri May 28 12:26:35 2021.+module Mcmc.Posterior+  ( Posterior,+  )+where++import Numeric.Log++-- | Posterior values are stored in log domain.+type Posterior = Log Double

@@ -2,8 +2,8 @@  -- | -- Module      :  Prior--- Description :  Convenience functions for computing priors--- Copyright   :  (c) Dominik Schrempf, 2020+-- Description :  Types and convenience functions for computing priors+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -12,7 +12,11 @@ -- -- Creation date: Thu Jul 23 13:26:14 2020. module Mcmc.Prior-  ( -- * Improper priors+  ( Prior,+    PriorFunction,+    noPrior,++    -- * Improper priors     largerThan,     positive,     lowerThan,@@ -21,6 +25,10 @@     -- * Continuous priors     exponential,     gamma,+    gammaMeanVariance,+    gammaMeanOne,+    gammaShapeScaleToMeanVariance,+    gammaMeanVarianceToShapeScale,     normal,     uniform, @@ -34,6 +42,7 @@  import Control.Monad import Data.Maybe (fromMaybe)+import Mcmc.Statistics.Types import Numeric.Log import qualified Statistics.Distribution as S import qualified Statistics.Distribution.Exponential as S@@ -41,80 +50,98 @@ import qualified Statistics.Distribution.Normal as S import qualified Statistics.Distribution.Poisson as S +-- | Prior values are stored in log domain.+type Prior = Log Double++-- | Prior function.+type PriorFunction a = a -> Prior++-- | Flat prior function. Useful for testing and debugging.+noPrior :: PriorFunction a+noPrior = const 1.0+ -- | Improper uniform prior; strictly larger than a given value.-largerThan :: Double -> Double -> Log Double+largerThan :: LowerBoundary -> PriorFunction Double largerThan a x   | x <= a = 0   | otherwise = 1  -- | Improper uniform prior; strictly larger than zero.-positive :: Double -> Log Double+positive :: PriorFunction Double positive = largerThan 0  -- | Improper uniform prior; strictly lower than a given value.-lowerThan :: Double -> Double -> Log Double+lowerThan :: UpperBoundary -> PriorFunction Double lowerThan b x   | x >= b = 0   | otherwise = 1  -- | Improper uniform prior; strictly lower than zero.-negative :: Double -> Log Double+negative :: PriorFunction Double negative = lowerThan 0  -- | Exponential distributed prior.-exponential ::-  -- | Rate.-  Double ->-  Double ->-  Log Double-exponential l x = Exp $ S.logDensity d x+exponential :: Rate -> PriorFunction Double+exponential l = Exp . S.logDensity d   where     d = S.exponential l  -- | Gamma distributed prior.-gamma ::-  -- | Shape.-  Double ->-  -- | Scale.-  Double ->-  Double ->-  Log Double-gamma k t x = Exp $ S.logDensity d x+gamma :: Shape -> Scale -> PriorFunction Double+gamma k t = Exp . S.logDensity d   where     d = S.gammaDistr k t +-- | See 'gamma' but parametrized using mean and variance.+gammaMeanVariance :: Mean -> Variance -> PriorFunction Double+gammaMeanVariance m v = Exp . S.logDensity d+  where (k, th) = gammaMeanVarianceToShapeScale m v+        d = S.gammaDistr k th++-- | Gamma disstributed prior with given shape and mean 1.0.+gammaMeanOne :: Shape -> PriorFunction Double+gammaMeanOne k = Exp . S.logDensity d+  where+    d = S.gammaDistr k (recip k)++-- The mean and variance of the gamma distribution are+--+-- m = k*t+--+-- v = k*t*t+--+-- Hence, the shape and scale are+--+-- k = m^2/v+--+-- t = v/m++-- | Calculate mean and variance of the gamma distribution given the shape and+-- the scale.+gammaShapeScaleToMeanVariance :: Shape -> Scale -> (Mean, Variance)+gammaShapeScaleToMeanVariance k t = let m = k * t in (m, m * t)++-- | Calculate shape and scale of the gamma distribution given the mean and+-- the variance.+gammaMeanVarianceToShapeScale :: Mean -> Variance -> (Shape, Scale)+gammaMeanVarianceToShapeScale m v = (m * m / v, v / m)+ -- | Normal distributed prior.-normal ::-  -- | Mean.-  Double ->-  -- | Standard deviation.-  Double ->-  Double ->-  Log Double-normal m s x = Exp $ S.logDensity d x+normal :: Mean -> StandardDeviation -> PriorFunction Double+normal m s = Exp . S.logDensity d   where     d = S.normalDistr m s  -- | Uniform prior on [a, b].-uniform ::-  -- | Lower bound a.-  Double ->-  -- | Upper bound b.-  Double ->-  Double ->-  Log Double+uniform :: LowerBoundary -> UpperBoundary -> PriorFunction Double uniform a b x   | x <= a = 0   | x >= b = 0   | otherwise = Exp 0  -- | Poisson distributed prior.-poisson ::-  -- | Rate.-  Double ->-  Int ->-  Log Double-poisson l x = Exp $ S.logProbability d x+poisson :: Rate -> PriorFunction Int+poisson l = Exp . S.logProbability d   where     d = S.poisson l 

@@ -6,7 +6,7 @@ -- | -- Module      :  Mcmc.Proposal -- Description :  Proposals are instruction to move around the state space--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -15,34 +15,39 @@ -- -- Creation date: Wed May 20 13:42:53 2020. module Mcmc.Proposal-  ( -- * Proposal+  ( -- * Proposals and types     PName (..),     PDescription (..),-    PWeight (..),+    PWeight (fromPWeight),+    pWeight,     PDimension (..),     Proposal (..),+    KernelRatio,+    Jacobian,+    JacobianFunction,     (@~),+    liftProposal,+    liftProposalWith,     ProposalSimple,-    Tuner (tParam, tFunc),+    Tuner (..),     Tune (..),     createProposal,+    TuningParameter,+    tuningParameterMin,+    tuningParameterMax,     tune,     getOptimalRate,-    proposalHeader,-    proposalHLine,-    summarizeProposal, -    -- * Cycle+    -- * Cycles     Order (..),     Cycle (ccProposals),     cycleFromList,     setOrder,-    orderProposals,+    prepareProposals,     tuneCycle,     autoTuneCycle,-    summarizeCycle, -    -- * Acceptance+    -- * Acceptance rates     Acceptance (fromAcceptance),     emptyA,     pushA,@@ -50,6 +55,12 @@     transformKeysA,     acceptanceRate,     acceptanceRates,++    -- * Output+    proposalHeader,+    proposalHLine,+    summarizeProposal,+    summarizeCycle,   ) where @@ -79,11 +90,16 @@ newtype PDescription = PDescription {fromPDescription :: String}   deriving (Show, Eq, Ord) --- | The weight determines how often a 'Proposal' is executed per iteration of--- the Markov chain.+-- | The positive weight determines how often a 'Proposal' is executed per+-- iteration of the Markov chain. newtype PWeight = PWeight {fromPWeight :: Int}   deriving (Show, Eq, Ord) +-- | Check if the weight is positive.+pWeight :: Int -> PWeight+pWeight n | n <= 0 = error "pWeight: Proposal weight is zero or negative."+          | otherwise = PWeight n+ -- | Proposal dimension. -- -- The number of affected, independent parameters.@@ -127,48 +143,90 @@ -- -- A 'Proposal' may be tuneable in that it contains information about how to enlarge -- or shrink the step size to tune the acceptance rate.+--+-- Predefined proposals are provided. To create custom proposals, one may use+-- the convenience function 'createProposal'. data Proposal a = Proposal   { -- | Name of the affected variable.-    pName :: PName,+    prName :: PName,     -- | Description of the proposal type and parameters.-    pDescription :: PDescription,+    prDescription :: PDescription,     -- | Dimension of the proposal. The dimension is used to calculate the     -- optimal acceptance rate, and does not have to be exact.-    pDimension :: PDimension,+    prDimension :: PDimension,     -- | The weight determines how often a 'Proposal' is executed per iteration of     -- the Markov chain.-    pWeight :: PWeight,+    prWeight :: PWeight,     -- | Simple proposal without name, weight, and tuning information.-    pSimple :: ProposalSimple a,+    prSimple :: ProposalSimple a,     -- | Tuning is disabled if set to 'Nothing'.-    pTuner :: Maybe (Tuner a)+    prTuner :: Maybe (Tuner a)   }  instance Eq (Proposal a) where-  m == n = pName m == pName n && pDescription m == pDescription n+  m == n = prName m == prName n && prDescription m == prDescription n  instance Ord (Proposal a) where-  compare = compare `on` (\p -> (pDescription p, pName p, pWeight p))+  compare = compare `on` (\p -> (prDescription p, prName p, prWeight p)) --- | Convert a proposal from one data type to another using a lens.+-- | Ratio of the proposal kernels. --+-- Part of the MHG acceptance ratio.+--+-- 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?+type KernelRatio = Log Double++-- | Absolute value of the determinant of the Jacobian matrix.+--+-- Part of the MHG acceptance ratio.+--+-- See also 'Jacobian'.+type Jacobian = Log Double++-- | Function calculating the 'Jacobian' of a proposal.+type JacobianFunction a = a -> Jacobian++-- | Lift a proposal from one data type to another.+--+-- Assume the Jacobian is 1.0 (see also 'liftProposal' and 'liftProposalWith').+-- -- For example: -- -- @ -- scaleFirstEntryOfTuple = _1 @~ scale -- @+infixl 7 @~ (@~) :: Lens' b a -> Proposal a -> Proposal b-(@~) l (Proposal n r d w s t) = Proposal n r d w (convertProposalSimple l s) (convertTuner l <$> t)+(@~) = liftProposal +-- | Lift a proposal from one data type to another.+--+-- Assume the Jacobian is 1.0 (see also '(@~)' and 'liftProposalWith').+liftProposal :: Lens' b a -> Proposal a -> Proposal b+liftProposal = liftProposalWith (const 1.0)++-- | Lift a proposal from one data type to another.+--+-- A function to calculate the Jacobian has to be provided (see also+-- 'liftProposal').+--+-- For further reference, please see the [example+-- @Pair@](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs).+liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b+liftProposalWith jf l (Proposal n r d w s t) = Proposal n r d w (convertProposalSimple jf l s) (convertTuner jf l <$> t)+ -- | Simple proposal without tuning information. -- -- Instruction about randomly moving from the current state to a new state, -- given some source of randomness. -- -- In order to calculate the Metropolis-Hastings-Green ratio, we need to know--- the ratio of the backward to forward kernels (i.e., the probability masses or--- probability densities) and the absolute value of the determinant of the--- Jacobian matrix.+-- the ratio of the backward to forward kernels (the 'KernelRatio' or the+-- probability masses or probability densities) and the 'Jacobian'. -- -- For unbiased proposals, these values are 1.0 such that --@@ -179,38 +237,45 @@ -- For biased proposals, the kernel ratio is qYX / qXY, where qXY is the -- probability density to move from X to Y, and the absolute value of the -- determinant of the Jacobian matrix differs from 1.0.-type ProposalSimple a = a -> GenIO -> IO (a, Log Double, Log Double)+type ProposalSimple a = a -> GenIO -> IO (a, KernelRatio, Jacobian) -convertProposalSimple :: Lens' b a -> ProposalSimple a -> ProposalSimple b-convertProposalSimple l s = s'+convertProposalSimple :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b+convertProposalSimple jf l s = s'   where-    s' v g = do-      (x', r, j) <- s (v ^. l) g-      return (set l x' v, r, j)+    s' y g = do+      (x', r, j) <- s (y ^. l) g+      let y' = set l x' y+          jxy = jf y+          jyx = jf y'+          j' = j * jyx / jxy+      return (y', r, j')  -- | Tune the acceptance rate of a 'Proposal'; see 'tune', or 'autoTuneCycle'. data Tuner a = Tuner-  { tParam :: Double,-    tFunc :: Double -> ProposalSimple a+  { tParam :: TuningParameter,+    tFunc :: TuningParameter -> ProposalSimple a   } -convertTuner :: Lens' b a -> Tuner a -> Tuner b-convertTuner l (Tuner p f) = Tuner p f'+convertTuner :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b+convertTuner jf l (Tuner p f) = Tuner p f'   where-    f' x = convertProposalSimple l $ f x+    f' x = convertProposalSimple jf l $ f x  -- | Tune the proposal? data Tune = Tune | NoTune   deriving (Show, Eq) +-- | Tuning parameter.+type TuningParameter = Double+ -- | Create a tuneable proposal. createProposal ::   -- | Description of the proposal type and parameters.   PDescription ->   -- | Function creating a simple proposal for a given tuning parameter. The-  -- larger the tuning parameter, the larger the proposal (and the lower the-  -- expected acceptance rate), and vice versa.-  (Double -> ProposalSimple a) ->+  -- larger the tuning parameter, the larger the proposal and the lower the+  -- expected acceptance rate; and vice versa.+  (TuningParameter -> ProposalSimple a) ->   -- | Dimension.   PDimension ->   -- | Name.@@ -223,19 +288,34 @@ createProposal r f d n w Tune = Proposal n r d w (f 1.0) (Just $ Tuner 1.0 f) createProposal r f d n w NoTune = Proposal n r d w (f 1.0) Nothing --- Minimal tuning parameter; subject to change.-tuningParamMin :: Double-tuningParamMin = 1e-12+-- | Minimal tuning parameter; @1e-12@, subject to change.+--+-- >>> tuningParameterMin+-- 1e-5+tuningParameterMin :: TuningParameter+tuningParameterMin = 1e-5 --- | Tune a 'Proposal'. Return 'Nothing' if 'Proposal' is not tuneable. The size---   of the proposal is proportional to the tuning parameter. Negative tuning---   parameters are not allowed.-tune :: (Double -> Double) -> Proposal a -> Maybe (Proposal a)+-- | Maximal tuning parameter; @1e12@, subject to change.+-- >>> tuningParameterMax+-- 1e3+tuningParameterMax :: TuningParameter+tuningParameterMax = 1e3++-- | Tune a 'Proposal'.+--+-- The size of the proposal is proportional to the tuning parameter which has a+-- positive lower bound of 'tuningParameterMin'.+--+-- The tuning function maps the current tuning parameter to a new one.+--+-- Return 'Nothing' if 'Proposal' is not tuneable.+tune :: (TuningParameter -> TuningParameter) -> Proposal a -> Maybe (Proposal a) tune f m = do-  (Tuner t g) <- pTuner m-  -- Ensure that the tuning parameter is strictly positive.-  let t' = max tuningParamMin (f t)-  return $ m {pSimple = g t', pTuner = Just $ Tuner t' g}+  (Tuner t g) <- prTuner m+  -- Ensure that the tuning parameter is strictly positive and well bounded.+  let t' = max tuningParameterMin (f t)+      t'' = min tuningParameterMax t'+  return $ m {prSimple = g t'', prTuner = Just $ Tuner t'' g}  -- | See 'PDimension'. getOptimalRate :: PDimension -> Double@@ -329,8 +409,8 @@ setOrder o c = c {ccOrder = o}  -- | Replicate 'Proposal's according to their weights and possibly shuffle them.-orderProposals :: Cycle a -> GenIO -> IO [Proposal a]-orderProposals (Cycle xs o) g = case o of+prepareProposals :: Cycle a -> GenIO -> IO [Proposal a]+prepareProposals (Cycle xs o) g = case o of   RandomO -> shuffle ps g   SequentialO -> return ps   RandomReversibleO -> do@@ -338,7 +418,7 @@     return $ psR ++ reverse psR   SequentialReversibleO -> return $ ps ++ reverse ps   where-    !ps = concat [replicate (fromPWeight $ pWeight p) p | p <- xs]+    !ps = concat [replicate (fromPWeight $ prWeight p) p | p <- xs]  -- The number of proposals depends on the order. getNProposalsPerCycle :: Cycle a -> Int@@ -348,10 +428,10 @@   RandomReversibleO -> 2 * once   SequentialReversibleO -> 2 * once   where-    once = sum $ map (fromPWeight . pWeight) xs+    once = sum $ map (fromPWeight . prWeight) xs  -- | Tune 'Proposal's in the 'Cycle'. See 'tune'.-tuneCycle :: M.Map (Proposal a) (Double -> Double) -> Cycle a -> Cycle a+tuneCycle :: M.Map (Proposal a) (TuningParameter -> TuningParameter) -> Cycle a -> Cycle a tuneCycle m c =   if sort (M.keys m) == sort ps     then c {ccProposals = map tuneF ps}@@ -368,9 +448,11 @@ autoTuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a autoTuneCycle a = tuneCycle (M.mapWithKey tuningF $ acceptanceRates a)   where-    tuningF proposal currentRate currentTuningParam =-      let optimalRate = getOptimalRate (pDimension proposal)-       in exp (currentRate - optimalRate) * currentTuningParam+    tuningF proposal mCurrentRate currentTuningParam = case mCurrentRate of+      Nothing -> currentTuningParam+      Just currentRate ->+        let optimalRate = getOptimalRate (prDimension proposal)+         in exp (2 * (currentRate - optimalRate)) * currentTuningParam  renderRow ::   BL.ByteString ->@@ -395,6 +477,67 @@     tp = alignRight 20 tuneParam     mt = alignRight 30 manualAdjustment +-- | For each key @k@, store the number of accepted and rejected proposals.+newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)}+  deriving (Eq, Read, Show)++instance ToJSONKey k => ToJSON (Acceptance k) where+  toJSON (Acceptance m) = toJSON m+  toEncoding (Acceptance m) = toEncoding m++instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where+  parseJSON v = Acceptance <$> 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, (0, 0)) | k <- ks]++-- | For key @k@, prepend an accepted (True) or rejected (False) proposal.+pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k+pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance+pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance+{-# INLINEABLE pushA #-}++-- | Reset acceptance storage.+resetA :: Ord k => Acceptance k -> Acceptance k+resetA = emptyA . M.keys . fromAcceptance++transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> M.Map k1 v -> M.Map k2 v+transformKeys ks1 ks2 m = foldl' insrt M.empty $ zip ks1 ks2+  where+    insrt m' (k1, k2) = M.insert k2 (m M.! k1) m'++-- | 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 ks1 ks2 = Acceptance . transformKeys ks1 ks2 . fromAcceptance++-- | Acceptance counts and rate for a specific proposal.+--+-- Return 'Nothing' if no proposals have been accepted or rejected (division by+-- zero).+acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, Double)+acceptanceRate k a = case fromAcceptance a M.!? k of+  Just (0, 0) -> Nothing+  Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))+  Nothing -> error "acceptanceRate: Key not found in map."++-- | 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 Double)+acceptanceRates =+  M.map+    ( \(as, rs) ->+        if as + rs == 0+          then Nothing+          else Just $ fromIntegral as / fromIntegral (as + rs)+    )+    . fromAcceptance+ -- | Header of proposal summaries. proposalHeader :: BL.ByteString proposalHeader =@@ -418,12 +561,11 @@   PName ->   PDescription ->   PWeight ->-  -- Tuning parameter.-  Maybe Double ->+  Maybe TuningParameter ->   PDimension ->   Maybe (Int, Int, Double) ->   BL.ByteString-summarizeProposal name description weight tuningParam dimension r =+summarizeProposal name description weight tuningParameter dimension ar =   renderRow     (BL.pack $ fromPName name)     (BL.pack $ fromPDescription description)@@ -436,16 +578,27 @@     manualAdjustmentStr   where     weightStr = BB.toLazyByteString $ BB.intDec $ fromPWeight weight-    nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) r-    nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) r-    acceptRate = BL.fromStrict $ maybe "" (BC.toFixed 2 . (^. _3)) r+    nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) ar+    nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) ar+    acceptRate = BL.fromStrict $ maybe "" (BC.toFixed 2 . (^. _3)) ar     optimalRate = BL.fromStrict $ BC.toFixed 2 $ getOptimalRate dimension-    tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) tuningParam-    check v-      | v < rateMin = "rate too low"-      | v > rateMax = "rate too high"-      | otherwise = ""-    manualAdjustmentStr = BL.fromStrict $ maybe "" (check . (^. _3)) r+    tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) tuningParameter+    checkRate rate+      | rate < rateMin = Just "rate too low"+      | rate > rateMax = Just "rate too high"+      | otherwise = Nothing+    checkTuningParam tp+      | tp <= (1.1 * tuningParameterMin) = Just "tuning parameter too low"+      | tp >= (0.9 * tuningParameterMax) = Just "tuning parameter too high"+      | otherwise = Nothing+    tps = checkTuningParam =<< tuningParameter+    ars = (checkRate . (^. _3)) =<< ar+    manualAdjustmentStr =+      let+       in case (ars, tps) of+            (Nothing, Nothing) -> ""+            (Just s, _) -> s+            (_, Just s) -> s  -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates. summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString@@ -458,11 +611,11 @@       proposalHLine     ]       ++ [ summarizeProposal-             (pName p)-             (pDescription p)-             (pWeight p)-             (tParam <$> pTuner p)-             (pDimension p)+             (prName p)+             (prDescription p)+             (prWeight p)+             (tParam <$> prTuner p)+             (prDimension p)              (ar p)            | p <- ps          ]@@ -475,51 +628,3 @@       1 -> nProposalsStr <> " proposal is performed per iteration."       _ -> nProposalsStr <> " proposals are performed per iterations."     ar m = acceptanceRate m a---- | For each key @k@, store the number of accepted and rejected proposals.-newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)}-  deriving (Eq, Read, Show)--instance ToJSONKey k => ToJSON (Acceptance k) where-  toJSON (Acceptance m) = toJSON m-  toEncoding (Acceptance m) = toEncoding m--instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where-  parseJSON v = Acceptance <$> 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, (0, 0)) | k <- ks]---- | For key @k@, prepend an accepted (True) or rejected (False) proposal.-pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k-pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance-pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance-{-# INLINEABLE pushA #-}---- | Reset acceptance storage.-resetA :: Ord k => Acceptance k -> Acceptance k-resetA = emptyA . M.keys . fromAcceptance--transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> M.Map k1 v -> M.Map k2 v-transformKeys ks1 ks2 m = foldl' insrt M.empty $ zip ks1 ks2-  where-    insrt m' (k1, k2) = M.insert k2 (m M.! k1) m'---- | 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 ks1 ks2 = Acceptance . transformKeys ks1 ks2 . fromAcceptance---- | Acceptance counts and rate for a specific proposal.-acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, Double)-acceptanceRate k a = case fromAcceptance a M.!? k of-  Just (0, 0) -> Nothing-  Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))-  Nothing -> error "acceptanceRate: Key not found in map."---- | Acceptance rates for all proposals.-acceptanceRates :: Acceptance k -> M.Map k Double-acceptanceRates = M.map (\(as, rs) -> fromIntegral as / fromIntegral (as + rs)) . fromAcceptance

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Proposal.Bactrian -- Description :  Bactrian proposals--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -12,21 +12,26 @@ -- -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. module Mcmc.Proposal.Bactrian-  ( slideBactrian,+  ( SpikeParameter,+    slideBactrian,     scaleBactrian,   ) where  import Mcmc.Proposal+import Mcmc.Statistics.Types import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Normal import System.Random.MWC import System.Random.MWC.Distributions +-- | Type synonym indicating the spike parameter.+type SpikeParameter = Double+ genBactrian ::-  Double ->-  Double ->+  SpikeParameter ->+  StandardDeviation ->   GenIO ->   IO Double genBactrian m s g = do@@ -37,7 +42,7 @@   b <- bernoulli 0.5 g   return $ if b then x else - x -logDensityBactrian :: Double -> Double -> Double -> Log Double+logDensityBactrian :: SpikeParameter -> StandardDeviation  -> Double -> Log Double logDensityBactrian m s x = Exp $ log $ kernel1 + kernel2   where     mn = m * s@@ -48,20 +53,18 @@     kernel2 = density dist2 x  bactrianAdditive ::-  Double ->-  Double ->-  Double ->-  GenIO ->-  IO (Double, Log Double, Log Double)+  SpikeParameter ->+  StandardDeviation  ->+  ProposalSimple Double bactrianAdditive m s x g = do   dx <- genBactrian m s g   return (x + dx, 1.0, 1.0)  -- bactrianSimple lens spike stdDev tune forwardOp backwardOp bactrianAdditiveSimple ::-  Double ->-  Double ->-  Double ->+  SpikeParameter ->+  StandardDeviation  ->+  TuningParameter  ->   ProposalSimple Double bactrianAdditiveSimple m s t   | m < 0 = error "bactrianAdditiveSimple: Spike parameter negative."@@ -81,15 +84,10 @@ -- -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. slideBactrian ::-  -- | Spike parameter \(m\).-  Double ->-  -- | Standard deviation \(s\).-  Double ->-  -- | Name.+  SpikeParameter ->+  StandardDeviation ->   PName ->-  -- | Weight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) (PDimension 1)@@ -106,11 +104,9 @@ fInv dx = recip (1 - dx) - 1  bactrianMult ::-  Double ->-  Double ->-  Double ->-  GenIO ->-  IO (Double, Log Double, Log Double)+  SpikeParameter ->+  StandardDeviation  ->+  ProposalSimple Double bactrianMult m s x g = do   du <- genBactrian m s g   let qXY = logDensityBactrian m s du@@ -119,7 +115,7 @@       jac = Exp $ log $ recip u   return (x * u, qYX / qXY, jac) -bactrianMultSimple :: Double -> Double -> Double -> ProposalSimple Double+bactrianMultSimple :: SpikeParameter -> StandardDeviation  -> TuningParameter  -> ProposalSimple Double bactrianMultSimple m s t   | m < 0 = error "bactrianMultSimple: Spike parameter negative."   | m >= 1 = error "bactrianMultSimple: Spike parameter 1.0 or larger."@@ -129,15 +125,10 @@ -- | Multiplicative proposal with kernel similar to the silhouette of a Bactrian -- camel. See 'slideBactrian'. scaleBactrian ::-  -- | Spike parameter.-  Double ->-  -- | Standard deviation.-  Double ->-  -- | Name.+  SpikeParameter ->+  StandardDeviation ->   PName ->-  -- | Weight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double scaleBactrian m s = createProposal description (bactrianMultSimple m s) (PDimension 1)

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Proposal.Generic -- Description :  Generic interface for creating proposals--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -51,7 +51,7 @@   --   -- Conversion to log domain is necessary, because some determinants of   -- Jacobians are very small (or large).-  Maybe (a -> Double -> Log Double) ->+  Maybe (a -> Double -> Jacobian) ->   ProposalSimple a genericContinuous d f mInv mJac x g = do   u <- genContVar d g

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Proposal.Scale -- Description :  Multiplicative proposals--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -20,12 +20,13 @@  import Mcmc.Proposal import Mcmc.Proposal.Generic+import Mcmc.Statistics.Types import Numeric.Log import Statistics.Distribution.Gamma  -- The actual proposal with tuning parameter. The tuning parameter does not -- change the mean.-scaleSimple :: Double -> Double -> Double -> ProposalSimple Double+scaleSimple :: Shape -> Scale -> TuningParameter -> ProposalSimple Double scaleSimple k th t =   genericContinuous     (gammaDistr (k / t) (th * t))@@ -37,10 +38,8 @@  -- | Multiplicative proposal with gamma distributed kernel. scale ::-  -- | Shape.-  Double ->-  -- | Scale.-  Double ->+  Shape ->+  Scale ->   PName ->   PWeight ->   Tune ->@@ -54,8 +53,7 @@ -- The scale of the gamma distribution is set to (shape)^{-1}, so that the mean -- of the gamma distribution is 1.0. scaleUnbiased ::-  -- | Shape.-  Double ->+  Shape ->   PName ->   PWeight ->   Tune ->@@ -64,7 +62,7 @@   where     description = PDescription $ "Scale unbiased; shape: " ++ show k -scaleContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)+scaleContrarilySimple :: Shape -> Scale -> TuningParameter -> ProposalSimple (Double, Double) scaleContrarilySimple k th t =   genericContinuous     (gammaDistr (k / t) (th * t))@@ -80,10 +78,8 @@ -- The two values are scaled contrarily so that their product stays constant. -- Contrary proposals are useful when parameters are confounded. scaleContrarily ::-  -- | Shape.-  Double ->-  -- | Scale.-  Double ->+  Shape ->+  Scale ->   PName ->   PWeight ->   Tune ->

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Proposal.Simplex -- Description :  Proposals on simplices--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -27,6 +27,7 @@ import Data.Aeson.TH import qualified Data.Vector.Unboxed as V import Mcmc.Proposal+import Mcmc.Statistics.Types import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Beta@@ -77,11 +78,12 @@ -- | Create the uniform element of the K-dimensional simplex. -- -- Set all values to \(1/D\).-simplexUniform :: Int -> Simplex+simplexUniform :: Dimension -> Simplex+--                 Should never fail. simplexUniform k = either error id $ simplexFromVector $ V.replicate k (1.0 / fromIntegral k)  -- Tuning function is inverted (high alpha means small steps).-getTuningFunction :: Double -> (Double -> Double)+getTuningFunction :: TuningParameter -> (TuningParameter -> TuningParameter) getTuningFunction t = (/ t'')   where     -- Start with small steps.@@ -98,7 +100,7 @@ -- The values determining the proposal size have been set using an example -- analysis. They are good values for this analysis, but may fail for other -- analyses.-dirichletSimple :: Double -> ProposalSimple Simplex+dirichletSimple :: TuningParameter -> ProposalSimple Simplex dirichletSimple t (SimplexUnsafe xs) g = do   -- If @t@ is high and above 1.0, the parameter vector will be low, and the   -- variance will be high. If @t@ is low and below 1.0, the parameter vector@@ -135,7 +137,7 @@ -- The proposal dimension, which is the dimension of the simplex, is used to -- determine the optimal acceptance rate. ----- For high dimensional simplices, this proposal may have low acceptance ratios.+-- For high dimensional simplices, this proposal may have low acceptance rates. -- In this case, please see the coordinate wise 'beta' proposal. dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple@@ -147,7 +149,7 @@ -- analyses. -- -- See also the 'dirichlet' proposal.-betaSimple :: Int -> Double -> ProposalSimple Simplex+betaSimple :: Dimension -> TuningParameter -> ProposalSimple Simplex betaSimple i t (SimplexUnsafe xs) g = do   -- Shape parameters of beta distribution. Do not assume that the sum of the   -- elements of 'xs' is 1.0, because then repeated proposals let the sum of the@@ -203,7 +205,7 @@ -- -- This proposal has been assigned a dimension of 2. See the discussion at -- 'PDimension'.-beta :: Int -> PName -> PWeight -> Tune -> Proposal Simplex+beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex beta i = createProposal description (betaSimple i) (PDimension 2)   where     description = PDescription $ "Beta; coordinate: " ++ show i

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Proposal.Slide -- Description :  Additive proposals--- Copyright   :  (c) Dominik Schrempf 2020+-- Copyright   :  (c) Dominik Schrempf 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -21,25 +21,21 @@  import Mcmc.Proposal import Mcmc.Proposal.Generic+import Mcmc.Statistics.Types import Statistics.Distribution.Normal import Statistics.Distribution.Uniform  -- The actual proposal with tuning parameter.-slideSimple :: Double -> Double -> Double -> ProposalSimple Double+slideSimple :: Mean -> StandardDeviation -> TuningParameter -> ProposalSimple Double slideSimple m s t =   genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing  -- | Additive proposal with normally distributed kernel. slide ::-  -- | Mean.-  Double ->-  -- | Standard deviation.-  Double ->-  -- | Name.+  Mean ->+  StandardDeviation ->   PName ->-  -- | Weight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double slide m s = createProposal description (slideSimple m s) (PDimension 1)@@ -47,7 +43,7 @@     description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s  -- The actual proposal with tuning parameter.-slideSymmetricSimple :: Double -> Double -> ProposalSimple Double+slideSymmetricSimple :: StandardDeviation -> TuningParameter -> ProposalSimple Double slideSymmetricSimple s t =   genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing @@ -55,13 +51,9 @@ -- proposal is very fast, because the Metropolis-Hastings-Green ratio does not -- include calculation of the forwards and backwards kernels. slideSymmetric ::-  -- | Standard deviation.-  Double ->-  -- | Name.+  StandardDeviation ->   PName ->-  -- | Weight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double slideSymmetric s = createProposal description (slideSymmetricSimple s) (PDimension 1)@@ -69,7 +61,7 @@     description = PDescription $ "Slide symmetric; sd: " ++ show s  -- The actual proposal with tuning parameter.-slideUniformSimple :: Double -> Double -> ProposalSimple Double+slideUniformSimple :: Size -> TuningParameter -> ProposalSimple Double slideUniformSimple d t =   genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing @@ -77,13 +69,9 @@ -- proposal is very fast, because the Metropolis-Hastings-Green ratio does not -- include calculation of the forwards and backwards kernels. slideUniformSymmetric ::-  -- | Delta.-  Double ->-  -- | Name.+  Size ->   PName ->-  -- | Weight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal Double slideUniformSymmetric d = createProposal description (slideUniformSimple d) (PDimension 1)@@ -93,7 +81,11 @@ contra :: (Double, Double) -> Double -> (Double, Double) contra (x, y) u = (x + u, y - u) -slideContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)+slideContrarilySimple ::+  Mean ->+  StandardDeviation ->+  TuningParameter ->+  ProposalSimple (Double, Double) slideContrarilySimple m s t =   genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing @@ -102,15 +94,10 @@ -- The two values are slid contrarily so that their sum stays constant. Contrary -- proposals are useful when parameters are confounded. slideContrarily ::-  -- | Mean.-  Double ->-  -- | Standard deviation.-  Double ->-  -- | Name.+  Mean ->+  StandardDeviation ->   PName ->-  -- | Weight.   PWeight ->-  -- | Enable tuning.   Tune ->   Proposal (Double, Double) slideContrarily m s = createProposal description (slideContrarilySimple m s) (PDimension 2)

@@ -1,10 +1,11 @@+{-# LANGUAGE DerivingVia #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-}  -- | -- Module      :  Mcmc.Settings -- Description :  Settings of Markov chain Monte Carlo samplers--- Copyright   :  (c) Dominik Schrempf, 2020+-- Copyright   :  (c) Dominik Schrempf, 2021 -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -15,13 +16,17 @@ module Mcmc.Settings   ( -- * Data types     AnalysisName (..),+    HasAnalysisName (..),     BurnInSpecification (..),     burnInIterations,     Iterations (..),+    TraceLength (..),     ExecutionMode (..),+    HasExecutionMode (..),     openWithExecutionMode,     ParallelizationMode (..),     SaveMode (..),+    LogMode (..),     Verbosity (..),      -- * Settings@@ -35,15 +40,21 @@ import Data.Aeson import Data.Aeson.TH import qualified Data.ByteString.Lazy.Char8 as BL+import Mcmc.Logger import System.Directory import System.IO  -- | Analysis name of the MCMC sampler. newtype AnalysisName = AnalysisName {fromAnalysisName :: String}   deriving (Eq, Read, Show)+  deriving (Monoid, Semigroup) via String  $(deriveJSON defaultOptions ''AnalysisName) +-- | Types with analysis names.+class HasAnalysisName s where+  getAnalysisName :: s -> AnalysisName+ -- | Burn in specification. data BurnInSpecification   = -- | No burn in.@@ -53,15 +64,32 @@   | -- | Burn in for a given number of iterations. Enable auto tuning with a     -- given period.     BurnInWithAutoTuning Int Int+  | -- | Burn in with the given list of auto tuning periods.+    --+    -- For example, @BurnInWithCustomAutoTuning [100,200]@ performs 300+    -- iterations with two auto tuning steps. One after 100 iterations, the+    -- second one after 200 more iterations.+    --+    -- Usually it is useful to auto tune more frequently in the beginning of the+    -- MCMC run.+    BurnInWithCustomAutoTuning [Int]   deriving (Eq, Read, Show)  $(deriveJSON defaultOptions ''BurnInSpecification) +-- Check if the burn in specification is valid.+burnInValid :: BurnInSpecification -> Bool+burnInValid NoBurnIn = True+burnInValid (BurnInWithoutAutoTuning n) = n > 0+burnInValid (BurnInWithAutoTuning n t) = n > 0 && t > 0+burnInValid (BurnInWithCustomAutoTuning xs) = not (null xs) && all (> 0) xs+ -- | Get the number of burn in iterations. burnInIterations :: BurnInSpecification -> Int burnInIterations NoBurnIn = 0 burnInIterations (BurnInWithoutAutoTuning n) = n burnInIterations (BurnInWithAutoTuning n _) = n+burnInIterations (BurnInWithCustomAutoTuning xs) = sum xs  -- | Number of normal iterations after burn in. --@@ -71,6 +99,17 @@  $(deriveJSON defaultOptions ''Iterations) +-- | The length of the stored "Mcmc.Chain.Trace".+--+-- Be careful, this setting determines the memory requirement of the MCMC chain.+data TraceLength+  = -- | Automatically determine the length of the trace. The value is+    -- determined by the 'Mcmc.Monitor.MonitorBatch' with largest batch size.+    TraceAuto+  | -- | Store a given minimum number of iterations of the chain. Store more+    --  iterations if required (see 'TraceAuto').+    TraceMinimum Int+ -- | Execution mode. data ExecutionMode   = -- | Perform new run.@@ -89,6 +128,10 @@  $(deriveJSON defaultOptions ''ExecutionMode) +-- | Types with execution modes.+class HasExecutionMode s where+  getExecutionMode :: s -> ExecutionMode+ -- | Open a file honoring the execution mode. -- -- Call 'error' if execution mode is@@ -121,17 +164,16 @@ -- | Parallelization mode. -- -- Parallel execution of the chains is only beneficial when the algorithm allows--- for parallelization, and if computation of the next iteration takes a long+-- for parallelization, and if computation of the next iteration takes some -- time. If the calculation of the next state is fast, sequential execution is--- usually beneficial, even for algorithms involving parallel chains. If the--- calculation of the next state is slow, parallel execution may be beneficial.+-- usually beneficial, even for algorithms involving parallel chains. ----- - The "Mcmc.Algorithm.Metropolis" algorithm is inherently sequential.+-- - The "Mcmc.Algorithm.MHG" algorithm is inherently sequential. -- -- - The "Mcmc.Algorithm.MC3" algorithm works well with parallelization. -- -- Of course, also the prior or likelihood functions can be computed in--- parallel. However, this library is not aware of how these functions are+-- parallel. However, this library is unaware about how these functions are -- computed. data ParallelizationMode   = Sequential@@ -140,20 +182,18 @@  $(deriveJSON defaultOptions ''ParallelizationMode) --- | Should the MCMC run be saved at the end of the run?-data SaveMode = NoSave | Save+-- | Define information stored on disk.+data SaveMode+  = -- | Do not save the MCMC analysis. The analysis can not be continued.+    NoSave+  | -- | Save the MCMC analysis so that it can be continued. This can be slow,+    -- if the trace is long, or if the states are large objects. See+    -- 'TraceLength'.+    Save   deriving (Eq, Read, Show)  $(deriveJSON defaultOptions ''SaveMode) --- $(deriveJSON defaultOptions ''CleaningMode)---- | Not much to say here.-data Verbosity = Quiet | Warn | Info | Debug-  deriving (Eq, Ord, Read, Show)--$(deriveJSON defaultOptions ''Verbosity)- -- | Settings of an MCMC sampler. data Settings = Settings   { sAnalysisName :: AnalysisName,@@ -162,14 +202,27 @@     sExecutionMode :: ExecutionMode,     sParallelizationMode :: ParallelizationMode,     sSaveMode :: SaveMode,+    sLogMode :: LogMode,     sVerbosity :: Verbosity   }   deriving (Eq, Show) +instance HasAnalysisName Settings where+  getAnalysisName = sAnalysisName++instance HasExecutionMode Settings where+  getExecutionMode = sExecutionMode++instance HasLogMode Settings where+  getLogMode = sLogMode++instance HasVerbosity Settings where+  getVerbosity = sVerbosity+ $(deriveJSON defaultOptions ''Settings)  settingsFn :: String -> FilePath-settingsFn n = n ++ ".settings"+settingsFn n = n ++ ".mcmc.settings"  -- | Save settings to a file determined by the analysis name. settingsSave :: Settings -> IO ()@@ -217,10 +270,10 @@   -- | Current iteration.   Int ->   IO ()-settingsCheck s@(Settings nm bi i em _ _ _) iCurrent+settingsCheck s@(Settings nm bi i em _ _ _ _) iCurrent   | null (fromAnalysisName nm) = serr "Analysis name is the empty string."   | burnInIterations bi < 0 = serr "Number of burn in iterations is negative."-  | not $ burnInAutoTuningPeriodValid bi = serr "Auto tuning period is zero or negative."+  | 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."@@ -231,6 +284,3 @@   | otherwise = return ()   where     serr = settingsError s iCurrent-    burnInAutoTuningPeriodValid :: BurnInSpecification -> Bool-    burnInAutoTuningPeriodValid (BurnInWithAutoTuning _ t) = t > 0-    burnInAutoTuningPeriodValid _ = True

@@ -0,0 +1,56 @@+-- |+-- Module      :  Mcmc.Statistics.Types+-- Description :  Types indicating properties of distributions+-- Copyright   :  (c) 2021 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Tue Feb 23 14:51:06 2021.+module Mcmc.Statistics.Types+  ( Mean,+    StandardDeviation,+    Variance,+    Shape,+    Scale,+    Rate,+    Dimension,+    Size,+    LowerBoundary,+    UpperBoundary,+  )+where++-- | Mean of a distribution.+type Mean = Double++-- | Standard deviation of a distribution.+type StandardDeviation = Double++-- | Variance of a distribution.+type Variance = Double++-- | Shape of a distribution.+type Shape = Double++-- | Scale of a distribution.+type Scale = Double++-- | Rate of a distribution.+type Rate = Double++-- | Dimension of a distribution.+type Dimension = Int++-- | Size of a distribution.+--+-- For example, the size of the interval of the uniform distribution.+type Size = Double++-- | Lower boundary of a distribution.+type LowerBoundary = Double++-- | Upper boundary of a distribution.+type UpperBoundary = Double

@@ -1,7 +1,7 @@ -- | --   Module      :  Mcmc.ProposalSpec --   Description :  Unit tests for Mcmc.Proposal---   Copyright   :  (c) Dominik Schrempf, 2020+--   Copyright   :  (c) Dominik Schrempf, 2021 --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -20,24 +20,24 @@ import Test.Hspec  p1 :: Proposal Double-p1 = slideSymmetric 1.0 (PName "Test 1") (PWeight 1) Tune+p1 = slideSymmetric 1.0 (PName "Test 1") (pWeight 1) Tune  p2 :: Proposal Double-p2 = slideSymmetric 1.0 (PName "Test 2") (PWeight 3) Tune+p2 = slideSymmetric 1.0 (PName "Test 2") (pWeight 3) Tune  c :: Cycle Double c = cycleFromList [p1, p2]  spec :: Spec spec =-  describe "orderProposals" $+  describe "prepareProposals" $     it "returns the correct number of proposals in a cycle" $       do         g <- create-        l1 <- length <$> orderProposals c g+        l1 <- length <$> prepareProposals c g         l1 `shouldBe` 4-        l2 <- length <$> orderProposals (setOrder RandomReversibleO c) g+        l2 <- length <$> prepareProposals (setOrder RandomReversibleO c) g         l2 `shouldBe` 8-        o3 <- orderProposals (setOrder SequentialReversibleO c) g+        o3 <- prepareProposals (setOrder SequentialReversibleO c) g         length o3 `shouldBe` 8         o3 == [p1, p2, p2, p2, p2, p2, p2, p1] `shouldBe` True

@@ -1,7 +1,7 @@ -- | --   Module      :  Mcmc.SaveSpec --   Description :  Unit tests for Mcmc.Save---   Copyright   :  (c) Dominik Schrempf, 2020+--   Copyright   :  (c) Dominik Schrempf, 2021 --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -18,7 +18,6 @@ import Mcmc.Chain.Chain import Mcmc.Chain.Save import Mcmc.Chain.Trace-import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Normal import qualified System.Random.MWC as R@@ -36,10 +35,10 @@ proposals :: Cycle Double proposals =   cycleFromList-    [ slideSymmetric 0.1 (PName "Small") (PWeight 5) Tune,-      slideSymmetric 1.0 (PName "Medium") (PWeight 2) Tune,-      slideSymmetric 5.0 (PName "Large") (PWeight 2) Tune,-      slide 1.0 4.0 (PName "Skewed") (PWeight 1) Tune+    [ slideSymmetric 0.1 (PName "Small") (pWeight 5) Tune,+      slideSymmetric 1.0 (PName "Medium") (pWeight 2) Tune,+      slideSymmetric 5.0 (PName "Large") (pWeight 2) Tune,+      slide 1.0 4.0 (PName "Skewed") (pWeight 1) Tune     ]  monStd :: MonitorStdOut Double@@ -62,8 +61,9 @@                 Overwrite                 Sequential                 NoSave+                LogStdOutOnly                 Quiet-        c <- fromMHG <$> mhg noPrior lh proposals mon 0 gen+        c <- fromMHG <$> mhg noPrior lh proposals mon TraceAuto 0 gen         savedChain <- toSavedChain c         c' <- fromSavedChain noPrior lh proposals mon savedChain         putStrLn "@load . save@ should be @id@."@@ -87,7 +87,7 @@         g2' <- R.save $ generator r'         g2 `shouldBe` g2' --- -- TODO.+-- -- TODO: 'mhContinue'. -- describe "mhContinue" --   $ it "mh 200 + mhContinue 200 == mh 400" --   $ do