mcmc 0.1.3 → 0.2.0

raw patch · 29 files changed

+1405/−1043 lines, 29 filesdep +data-default

Dependencies added: data-default

Files

• README.md +3/−2
• bench/Bench.hs +1/−0
• bench/Normal.hs +22/−9
• bench/Poisson.hs +9/−9
• mcmc.cabal +10/−6
• src/Mcmc.hs +55/−42
• src/Mcmc/Mcmc.hs +148/−50
• src/Mcmc/Metropolis.hs +51/−65
• src/Mcmc/Monitor.hs +70/−44
• src/Mcmc/Monitor/Parameter.hs +25/−25
• src/Mcmc/Monitor/ParameterBatch.hs +31/−45
• src/Mcmc/Monitor/Time.hs +8/−0
• src/Mcmc/Move.hs +0/−285
• src/Mcmc/Move/Generic.hs +0/−136
• src/Mcmc/Move/Scale.hs +0/−69
• src/Mcmc/Move/Slide.hs +0/−95
• src/Mcmc/Proposal.hs +344/−0
• src/Mcmc/Proposal/Bactrian.hs +137/−0
• src/Mcmc/Proposal/Generic.hs +90/−0
• src/Mcmc/Proposal/Scale.hs +57/−0
• src/Mcmc/Proposal/Slide.hs +82/−0
• src/Mcmc/Save.hs +35/−31
• src/Mcmc/Status.hs +101/−65
• src/Mcmc/Tools/Shuffle.hs +9/−0
• src/Mcmc/Trace.hs +0/−2
• src/Mcmc/Verbosity.hs +40/−0
• test/Mcmc/Move/SlideSpec.hs +0/−34
• test/Mcmc/ProposalSpec.hs +43/−0
• test/Mcmc/SaveSpec.hs +34/−29

@@ -3,7 +3,7 @@  <p align="center"><img src="https://travis-ci.org/dschrempf/mcmc.svg?branch=master"/></p> -Sample from a posterior using Markov chain Monte Carlo.+Sample from a posterior using Markov chain Monte Carlo methods.  At the moment, the library is tailored to the Metropolis-Hastings algorithm since it covers most use cases. However, implementation of more algorithms is@@ -12,7 +12,7 @@  ## Documentation -The source code contains detailed documentation about general concepts as well+The [source code](https://hackage.haskell.org/package/mcmc) contains detailed documentation about general concepts as well as specific functions.  @@ -22,6 +22,7 @@ attached to this repository.      git clone https://github.com/dschrempf/mcmc.git+    cd mcmc     stack build  For example, estimate the [accuracy of an archer](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Archery/Archery.hs) with

@@ -23,5 +23,6 @@   g <- create   defaultMain     [ bench "Normal" $ nfIO (normalBench g),+      bench "NormalBactrian" $ nfIO (normalBactrianBench g),       bench "Poisson" $ nfIO (poissonBench g)     ]

@@ -11,6 +11,7 @@ -- Creation date: Wed May  6 00:10:11 2020. module Normal   ( normalBench,+    normalBactrianBench,   ) where @@ -33,18 +34,20 @@ lh :: Double -> Log Double lh = Exp . logDensity (normalDistr trueMean trueStdDev) -moveCycle :: Cycle Double-moveCycle =+proposals :: Cycle Double+proposals =   fromList-    [ slideSymmetric "small" 5 id 0.1 True,-      slideSymmetric "medium" 2 id 1.0 True,-      slideSymmetric "large" 2 id 5.0 True,-      slide "skewed" 1 id 1.0 4.0 True-    ]+    [slideSymmetric "medium" 1 1.0 True] +mons :: [MonitorParameter Double]+mons = [monitorRealFloat "mu"]+ monStd :: MonitorStdOut Double-monStd = monitorStdOut [monitorRealFloat "mu" id] 200+monStd = monitorStdOut mons 200 +-- monFile :: MonitorFile Double+-- monFile = monitorFile "Mu" mons 200+ mon :: Monitor Double mon = Monitor monStd [] [] @@ -59,5 +62,15 @@  normalBench :: GenIO -> IO () normalBench g = do-  let s = noSave $ status "Normal" (const 1) lh moveCycle mon 0 nBurn nAutoTune nIter g+  let s = quiet $ noSave $ status "Normal" (const 1) lh proposals mon 0 nBurn nAutoTune nIter g+  void $ mh s++proposalsBactrian :: Cycle Double+proposalsBactrian =+  fromList+    [slideBactrian "bactrian" 1 0.5 1.0 True]++normalBactrianBench :: GenIO -> IO ()+normalBactrianBench g = do+  let s = quiet $ noSave $ status "NormalBactrian" (const 1) lh proposalsBactrian mon 0 nBurn nAutoTune nIter g   void $ mh s

@@ -47,23 +47,23 @@ lh x =   product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears] -moveAlpha :: Move I-moveAlpha = slideSymmetric "alpha" 2 _1 0.2 False+proposalAlpha :: Proposal I+proposalAlpha = _1 @~ slideSymmetric "alpha" 2 0.2 False -moveBeta :: Move I-moveBeta = slideSymmetric "beta" 1 _2 0.2 False+proposalBeta :: Proposal I+proposalBeta = _2 @~ slideSymmetric "beta" 1 0.2 False -moveCycle :: Cycle I-moveCycle = fromList [moveAlpha, moveBeta]+proposals :: Cycle I+proposals = fromList [proposalAlpha, proposalBeta]  initial :: I initial = (0, 0)  monAlpha :: MonitorParameter I-monAlpha = monitorRealFloat "alpha" _1+monAlpha = _1 @. monitorRealFloat "alpha"  monBeta :: MonitorParameter I-monBeta = monitorRealFloat "beta" _2+monBeta = _2 @. monitorRealFloat "beta"  monStd :: MonitorStdOut I monStd = monitorStdOut [monAlpha, monBeta] 150@@ -82,5 +82,5 @@  poissonBench :: GenIO -> IO () poissonBench g = do-  let s = noSave $ status "Poisson" (const 1) lh moveCycle mon initial nBurn nAutoTune nIter g+  let s = quiet $ noSave $ status "Poisson" (const 1) lh proposals mon initial nBurn nAutoTune nIter g   void $ mh s

@@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack  name:           mcmc-version:        0.1.3+version:        0.2.0 synopsis:       Sample from a posterior using Markov chain Monte Carlo description:    Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme> category:       Math, Statistics@@ -35,14 +35,16 @@       Mcmc.Monitor.Parameter       Mcmc.Monitor.ParameterBatch       Mcmc.Monitor.Time-      Mcmc.Move-      Mcmc.Move.Generic-      Mcmc.Move.Scale-      Mcmc.Move.Slide+      Mcmc.Proposal+      Mcmc.Proposal.Bactrian+      Mcmc.Proposal.Generic+      Mcmc.Proposal.Scale+      Mcmc.Proposal.Slide       Mcmc.Save       Mcmc.Status       Mcmc.Tools.Shuffle       Mcmc.Trace+      Mcmc.Verbosity   other-modules:       Paths_mcmc   autogen-modules:@@ -55,6 +57,7 @@     , base >=4.7 && <5     , bytestring     , containers+    , data-default     , directory     , log-domain     , microlens@@ -71,7 +74,7 @@   type: exitcode-stdio-1.0   main-is: Spec.hs   other-modules:-      Mcmc.Move.SlideSpec+      Mcmc.ProposalSpec       Mcmc.SaveSpec       Paths_mcmc   hs-source-dirs:@@ -80,6 +83,7 @@   build-depends:       QuickCheck     , base >=4.7 && <5+    , directory     , hspec     , hspec-discover     , log-domain

@@ -19,71 +19,81 @@ -- For examples, please see -- [mcmc-examples](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples). ----- The import of this module alone should cover most use cases.+-- __The import of this module alone should cover most use cases.__ module Mcmc-  ( -- * Moves+  ( -- * Proposals -    -- | A 'Move' is an instruction about how to advance a given Markov chain so-    -- that it possibly reaches a new state. That is, 'Move's specify how the+    -- | A 'Proposal' is an instruction about how to advance a given Markov chain so+    -- that it possibly reaches a new state. That is, 'Proposal's specify how the     -- chain traverses the state space. As far as this MCMC library is-    -- concerned, 'Move's are /elementary updates/ in that they cannot be+    -- concerned, 'Proposal's are /elementary updates/ in that they cannot be     -- decomposed into smaller updates.     ---    -- 'Move's can be combined to form composite updates, a technique often+    -- 'Proposal's can be combined to form composite updates, a technique often     -- referred to as /composition/. On the other hand, /mixing/ (used in the-    -- sense of mixture models) is the random choice of a 'Move' (or a-    -- composition of 'Move's) from a given set.+    -- sense of mixture models) is the random choice of a 'Proposal' (or a+    -- composition of 'Proposal's) from a given set.     ---    -- The __composition__ and __mixture__ of 'Move's allows specification of+    -- The __composition__ and __mixture__ of 'Proposal's allows specification of     -- nearly all MCMC algorithms involving a single chain (i.e., population     -- methods such as particle filters are excluded). In particular, Gibbs     -- samplers of all sorts can be specified using this procedure.     ---    -- This library enables composition and mixture of 'Move's via the 'Cycle'-    -- data type. Essentially, a 'Cycle' is a set of 'Move's. The chain advances+    -- This library enables composition and mixture of 'Proposal's via the 'Cycle'+    -- data type. Essentially, a 'Cycle' is a set of 'Proposal's. The chain advances     -- after the completion of each 'Cycle', which is called an __iteration__,     -- and the iteration counter is increased by one.     ---    -- The 'Move's in a 'Cycle' can be executed in the given order or in a+    -- The 'Proposal's in a 'Cycle' can be executed in the given order or in a     -- random sequence which allows, for example, specification of a fixed scan-    -- Gibbs sampler, or a random sequence scan Gibbs sampler, respectively.+    -- 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.     ---    -- Moves are named according to what they do, i.e., how they change the+    -- 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     -- probability distributions. For example, 'slideSymmetric' is a sliding-    -- move. Under the hood, it uses the normal distribution with mean zero and+    -- proposal. Under the hood, it uses the normal distribution with mean zero and     -- given variance. The sampled variate is added to the current value of the     -- variable (hence, the name slide). The same nomenclature is used by     -- RevBayes [1]. The probability distributions and intrinsic properties of a-    -- specific move are specified in detail in the documentation.+    -- specific proposal are specified in detail in the documentation.     --     -- The other method, which is used intrinsically, is more systematic, but     -- also a little bit more complicated: we separate between the proposal     -- distribution and how the state is affected. And here, I am not only     -- referring to the accessor (i.e., the lens), but also to the operator     -- (addition, multiplication, any other binary operator). For example, the-    -- sliding move (without tuning information) is implemented as+    -- sliding proposal (without tuning information) is implemented as     --     -- @-    -- slideSimple :: Lens' a Double -> Double -> Double -> Double -> MoveSimple a-    -- slideSimple l m s t = moveGenericContinuous l (normalDistr m (s * t)) (+) (-)+    -- slideSimple :: Lens' a Double -> Double -> Double -> Double -> ProposalSimple a+    -- slideSimple l m s t = proposalGenericContinuous l (normalDistr m (s * t)) (+) (-)     -- @     --     -- This specification is more involved. Especially since we need to know the     -- probability of jumping back, and so we need to know the inverse operator.-    -- However, it also allows specification of new moves with great ease.+    -- However, it also allows specification of new proposals with great ease.     --     -- [1] Höhna, S., Landis, M. J., Heath, T. A., Boussau, B., Lartillot, N., Moore,     -- B. R., Huelsenbeck, J. P., …, Revbayes: bayesian phylogenetic inference using     -- graphical models and an interactive model-specification language, Systematic     -- Biology, 65(4), 726–736 (2016). http://dx.doi.org/10.1093/sysbio/syw021-    module Mcmc.Move,-    module Mcmc.Move.Slide,-    module Mcmc.Move.Scale,+    Proposal,+    (@~),+    scale,+    scaleUnbiased,+    slide,+    slideBactrian,+    slideSymmetric,+    slideUniform,+    Cycle,+    fromList,+    Order (..),+    setOrder,      -- * Initialization @@ -91,16 +101,19 @@     -- constructed using the function 'status'.     --     -- The 'Status' of a Markov chain includes information about current state-    -- ('Item') and iteration, the history of the chain ('Trace'), the 'Acceptance'-    -- ratios, and the random number generator.+    -- ('Mcmc.Item.Item') and iteration, the history of the chain+    -- ('Mcmc.Trace.Trace'), the 'Acceptance' ratios, and the random number+    -- generator.     --     -- Further, the 'Status' includes auxiliary variables and functions such as     -- the prior and likelihood functions, instructions to move around the state     -- space (see above) and to monitor the MCMC run, as well as some auxiliary     -- information.-    module Mcmc.Status,-    module Mcmc.Item,-    module Mcmc.Trace,+    status,+    noSave,+    force,+    quiet,+    debug,      -- * Monitor @@ -110,7 +123,13 @@     -- - 'MonitorFile': Log to a file.     -- - 'MonitorBatch': Log summary statistics such as the mean of the last     -- - states to a file.-    module Mcmc.Monitor,+    Monitor (Monitor),+    MonitorStdOut,+    monitorStdOut,+    MonitorFile,+    monitorFile,+    MonitorBatch,+    monitorBatch,     module Mcmc.Monitor.Parameter,     module Mcmc.Monitor.ParameterBatch, @@ -119,27 +138,21 @@     -- | At the moment, the library is tailored to the Metropolis-Hastings     -- algorithm ('mh') since it covers most use cases. However, implementation     -- of more algorithms is planned in the future.-    module Mcmc.Metropolis,+    mh,+    mhContinue,      -- * Misc-    pzero,-    module Mcmc.Save,+    loadStatus,   ) where -import Mcmc.Item import Mcmc.Metropolis import Mcmc.Monitor import Mcmc.Monitor.Parameter import Mcmc.Monitor.ParameterBatch-import Mcmc.Move-import Mcmc.Move.Scale-import Mcmc.Move.Slide+import Mcmc.Proposal+import Mcmc.Proposal.Bactrian+import Mcmc.Proposal.Scale+import Mcmc.Proposal.Slide import Mcmc.Save import Mcmc.Status-import Mcmc.Trace-import Numeric.Log---- | Because we need a probability of zero for likelihoods of really bad moves.-pzero :: Fractional a => Log a-pzero = Exp $ - (1 / 0)

@@ -13,70 +13,160 @@ -- Creation date: Fri May 29 10:19:45 2020. module Mcmc.Mcmc   ( Mcmc,+    mcmcOutT,+    mcmcOutS,+    mcmcWarnT,+    mcmcWarnS,+    mcmcInfoT,+    mcmcInfoS,+    mcmcDebugT,+    mcmcDebugS,     mcmcAutotune,+    mcmcResetA,     mcmcSummarizeCycle,-    mcmcInit,     mcmcReport,-    mcmcMonitorHeader,+    mcmcMonitorStdOutHeader,     mcmcMonitorExec,-    mcmcClose,+    mcmcRun,   ) where  import Control.Monad import Control.Monad.IO.Class-import Control.Monad.Trans.State.Strict+import Control.Monad.Trans.State import Data.Aeson import Data.Maybe+import qualified Data.Text.Lazy as T+import Data.Text.Lazy (Text) import qualified Data.Text.Lazy.IO as T import Data.Time.Clock import Data.Time.Format import Mcmc.Monitor import Mcmc.Monitor.Time-import Mcmc.Move+import Mcmc.Proposal import Mcmc.Save-import Mcmc.Status+import Mcmc.Status hiding (debug)+import Mcmc.Verbosity+import System.Directory+import System.IO import Prelude hiding (cycle)  -- | An Mcmc state transformer; usually fiddling around with this type is not -- required, but it is used by the different inference algorithms. type Mcmc a = StateT (Status a) IO --- | Auto tune the 'Move's in the 'Cycle' of the chain. See 'autotune'.-mcmcAutotune :: Int -> Mcmc a ()-mcmcAutotune t = do+msgPrepare :: Char -> Text -> Text+msgPrepare c t = T.cons c $ ": " <> t++-- | Write to standard output and log file.+mcmcOutT :: Text -> Mcmc a ()+mcmcOutT msg = do+  h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> gets logHandle+  liftIO $ T.putStrLn msg >> T.hPutStrLn h msg++-- | Write to standard output and log file.+mcmcOutS :: String -> Mcmc a ()+mcmcOutS = mcmcOutT . T.pack++-- Perform warning action.+mcmcWarnA :: Mcmc a () -> Mcmc a ()+mcmcWarnA a = gets verbosity >>= \v -> info v a++-- | Print warning message.+mcmcWarnT :: Text -> Mcmc a ()+mcmcWarnT = mcmcWarnA . mcmcOutT . msgPrepare 'W'++-- | Print warning message.+mcmcWarnS :: String -> Mcmc a ()+mcmcWarnS = mcmcWarnT . T.pack++-- Perform info action.+mcmcInfoA :: Mcmc a () -> Mcmc a ()+mcmcInfoA a = gets verbosity >>= \v -> info v a++-- | Print info message.+mcmcInfoT :: Text -> Mcmc a ()+mcmcInfoT = mcmcInfoA . mcmcOutT . msgPrepare 'I'++-- | Print info message.+mcmcInfoS :: String -> Mcmc a ()+mcmcInfoS = mcmcInfoT . T.pack++-- Perform debug action.+mcmcDebugA :: Mcmc a () -> Mcmc a ()+mcmcDebugA a = gets verbosity >>= \v -> debug v a++-- | Print debug message.+mcmcDebugT :: Text -> Mcmc a ()+mcmcDebugT = mcmcDebugA . mcmcOutT . msgPrepare 'D'++-- | Print debug message.+mcmcDebugS :: String -> Mcmc a ()+mcmcDebugS = mcmcDebugT . T.pack++-- | Auto tune the 'Proposal's in the 'Cycle' of the chain. Reset acceptance counts.+-- See 'autotuneCycle'.+mcmcAutotune :: Mcmc a ()+mcmcAutotune = do+  mcmcDebugT "Auto tune."   s <- get   let a = acceptance s       c = cycle s-      c' = autotuneC t a c+      c' = autotuneCycle a c   put $ s {cycle = c'} --- | Print short summary of 'Move's in 'Cycle'. See 'summarizeCycle'.-mcmcSummarizeCycle :: Maybe Int -> Mcmc a ()-mcmcSummarizeCycle Nothing = do-  c <- gets cycle-  liftIO $ T.putStr $ summarizeCycle Nothing c-mcmcSummarizeCycle (Just n) = do+-- | Reset acceptance counts.+mcmcResetA :: Mcmc a ()+mcmcResetA = do+  mcmcDebugT "Reset acceptance ratios."+  s <- get+  let a = acceptance s+  put $ s {acceptance = resetA a}++-- | Print short summary of 'Proposal's in 'Cycle'. See 'summarizeCycle'.+mcmcSummarizeCycle :: Mcmc a Text+mcmcSummarizeCycle = do   a <- gets acceptance   c <- gets cycle-  liftIO $ T.putStr $ summarizeCycle (Just (n, a)) c+  return $ summarizeCycle a c  fTime :: FormatTime t => t -> String fTime = formatTime defaultTimeLocale "%B %-e, %Y, at %H:%M %P, %Z." --- | Set the total number of iterations, the current time and open the--- 'Monitor's of the chain. See 'mOpen'.+-- Open log file.+mcmcOpenLog :: Mcmc a ()+mcmcOpenLog = do+  s <- get+  let lfn = name s ++ ".log"+      n = iteration s+      frc = forceOverwrite s+  fe <- liftIO $ doesFileExist lfn+  mh <- liftIO $ case verbosity s of+    Quiet -> return Nothing+    _ -> Just <$> case (fe, n, frc) of+      (False, _, _) -> openFile lfn WriteMode+      (True, 0, True) -> openFile lfn WriteMode+      (True, 0, False) -> error "mcmcInit: Log file exists; use 'force' to overwrite output files."+      (True, _, _) -> openFile lfn AppendMode+  put s {logHandle = mh}+  mcmcDebugS $ "Log file name: " ++ lfn ++ "."+  mcmcDebugT "Log file opened."++-- Set the total number of iterations, the current time and open the 'Monitor's+-- of the chain. See 'mOpen'. mcmcInit :: Mcmc a () mcmcInit = do-  t <- liftIO getCurrentTime-  liftIO $ putStrLn $ "-- Start time: " <> fTime t+  mcmcOpenLog   s <- get+  -- Start time.+  t <- liftIO getCurrentTime+  mcmcInfoS $ "Start time: " <> fTime t+  -- Monitor.   let m = monitor s       n = iteration s       nm = name s-  m' <- if n == 0-        then liftIO $ mOpen nm m-        else liftIO $ mAppend nm m+      frc = forceOverwrite s+  m' <- if n == 0 then liftIO $ mOpen nm frc m else liftIO $ mAppend nm m   put $ s {monitor = m', start = Just (n, t)}  -- | Report what is going to be done.@@ -87,35 +177,32 @@       t = autoTuningPeriod s       n = iterations s   case b of-    Just b' -> liftIO $ putStrLn $ "-- Burn in for " <> show b' <> " iterations."+    Just b' -> mcmcInfoS $ "Burn in for " <> show b' <> " iterations."     Nothing -> return ()   case t of-    Just t' ->-      liftIO $ putStrLn $-        "-- Auto tune every "-          <> show t'-          <> " iterations (during burn in only)."+    Just t' -> mcmcInfoS $ "Auto tune every " <> show t' <> " iterations (during burn in only)."     Nothing -> return ()-  liftIO $ putStrLn $ "-- Run chain for " <> show n <> " iterations."-  liftIO $ putStrLn "-- Initial state."-  mcmcMonitorHeader+  mcmcInfoS $ "Run chain for " <> show n <> " iterations."+  mcmcInfoT "Initial state."+  mcmcMonitorStdOutHeader   mcmcMonitorExec --- | Print header line of 'Monitor' (only standard output).-mcmcMonitorHeader :: Mcmc a ()-mcmcMonitorHeader = do+-- | Print header line of standard output monitor.+mcmcMonitorStdOutHeader :: Mcmc a ()+mcmcMonitorStdOutHeader = do   m <- gets monitor-  liftIO $ mHeader m+  mcmcInfoA $ mcmcOutT $ mHeader m  -- Save the status of an MCMC run. See 'saveStatus'. mcmcSave :: ToJSON a => Mcmc a () mcmcSave = do   s <- get-  liftIO $ if save s-    then do putStrLn "-- Save Markov chain. For long chains, this may take a while."-            saveStatus (name s <> ".mcmc") s-            putStrLn "-- Done saving Markov chain."-    else putStrLn "-- Do not save the Markov chain."+  if save s+    then do+      mcmcInfoT "Save Markov chain. For long chains, this may take a while."+      liftIO $ saveStatus (name s <> ".mcmc") s+      mcmcInfoT "Done saving Markov chain."+    else mcmcInfoT "Do not save the Markov chain."  -- | Execute the 'Monitor's of the chain. See 'mExec'. mcmcMonitorExec :: ToJSON a => Mcmc a ()@@ -124,17 +211,18 @@   let i = iteration s       j = iterations s + fromMaybe 0 (burnInIterations s)       m = monitor s-      st = fromMaybe (error "mcmcMonitorExec: Starting state and time not set.") (start s)+      (ss, st) = fromMaybe (error "mcmcMonitorExec: Starting state and time not set.") (start s)       tr = trace s-  liftIO $ mExec i st tr j m+      vb = verbosity s+  mt <- liftIO $ mExec vb i ss st tr j m+  forM_ mt mcmcOutT --- | Close the 'Monitor's of the chain. See 'mClose'.+-- Close the 'Monitor's of the chain. See 'mClose'. mcmcClose :: ToJSON a => Mcmc a () mcmcClose = do   s <- get-  let n = iterations s-  mcmcSummarizeCycle (Just n)-  liftIO $ putStrLn "-- Metropolis-Hastings sampler finished."+  mcmcSummarizeCycle >>= mcmcInfoT+  mcmcInfoT "Metropolis-Hastings sampler finished."   let m = monitor s   m' <- liftIO $ mClose m   put $ s {monitor = m'}@@ -143,5 +231,15 @@   let rt = case start s of         Nothing -> error "mcmcClose: Start time not set."         Just (_, st) -> t `diffUTCTime` st-  liftIO $ T.putStrLn $ "-- Wall clock run time: " <> renderDuration rt <> "."-  liftIO $ putStrLn $ "-- End time: " <> fTime t+  mcmcInfoT $ "Wall clock run time: " <> renderDuration rt <> "."+  mcmcInfoS $ "End time: " <> fTime t+  case logHandle s of+    Just h -> liftIO $ hClose h+    Nothing -> return ()++-- | Run an MCMC algorithm.+mcmcRun :: ToJSON a => Mcmc a () -> Status a -> IO (Status a)+mcmcRun algorithm = execStateT $ do+  mcmcInit+  algorithm+  mcmcClose

@@ -1,4 +1,5 @@ {-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}  -- | -- Module      :  Mcmc.Metropolis@@ -19,33 +20,28 @@  import Control.Monad import Control.Monad.IO.Class-import Control.Monad.Trans.State.Strict+import Control.Monad.Trans.State import Data.Aeson import Data.Maybe import Mcmc.Item import Mcmc.Mcmc-import Mcmc.Move+import Mcmc.Proposal import Mcmc.Status-import Mcmc.Tools.Shuffle import Mcmc.Trace import Numeric.Log import System.Random.MWC import Prelude hiding (cycle) --- For non-symmetric moves.-mhRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double-mhRatio lX lY qXY qYX = lY * qYX / lX / qXY+-- For non-symmetric proposals.+--+-- q = qYX / qXY+mhRatio :: Log Double -> Log Double -> Log Double -> Log Double+mhRatio fX fY q = fY * q / fX {-# INLINE mhRatio #-} --- For symmetric moves.-mhRatioSymmetric :: Log Double -> Log Double -> Log Double-mhRatioSymmetric lX lY = lY / lX-{-# INLINE mhRatioSymmetric #-}--mhMove :: Move a -> Mcmc a ()-mhMove m = do-  let p = mvSample $ mvSimple m-      mq = mvDensity $ mvSimple m+mhPropose :: Proposal a -> Mcmc a ()+mhPropose m = do+  let p = pSample $ pSimple m   s <- get   let (Item x pX lX) = item s       pF = priorF s@@ -53,13 +49,11 @@       a = acceptance s       g = generator s   -- 1. Sample new state.-  !y <- liftIO $ p x g+  (!y, !q) <- liftIO $ p x g   -- 2. Calculate Metropolis-Hastings ratio.   let !pY = pF y       !lY = lF y-      !r = case mq of-        Nothing -> mhRatioSymmetric (pX * lX) (pY * lY)-        Just q -> mhRatio (pX * lX) (pY * lY) (q x y) (q y x)+      !r = mhRatio (pX * lX) (pY * lY) q   -- 3. Accept or reject.   if ln r >= 0.0     then put $ s {item = Item y pY lY, acceptance = pushA m True a}@@ -69,18 +63,12 @@         then put $ s {item = Item y pY lY, acceptance = pushA m True a}         else put $ s {acceptance = pushA m False a} --- TODO: Split the generator here. See SaveSpec -> mhContinue.---- Replicate 'Move's according to their weights and shuffle them.-getNCycles :: Cycle a -> Int -> GenIO -> IO [[Move a]]-getNCycles c = shuffleN mvs-  where-    !mvs = concat [replicate (mvWeight m) m | m <- fromCycle c]+-- TODO: Splitmix. Split the generator here. See SaveSpec -> mhContinue. --- Run one iterations; perform all moves in a Cycle.-mhIter :: ToJSON a => [Move a] -> Mcmc a ()-mhIter mvs = do-  mapM_ mhMove mvs+-- Run one iterations; perform all proposals in a Cycle.+mhIter :: ToJSON a => [Proposal a] -> Mcmc a ()+mhIter ps = do+  mapM_ mhPropose ps   s <- get   let i = item s       t = trace s@@ -91,6 +79,7 @@ -- Run N iterations. mhNIter :: ToJSON a => Int -> Mcmc a () mhNIter n = do+  mcmcDebugS $ "Run " <> show n <> " iterations."   c <- gets cycle   g <- gets generator   cycles <- liftIO $ getNCycles c n g@@ -100,8 +89,17 @@ mhBurnInN :: ToJSON a => Int -> Maybe Int -> Mcmc a () mhBurnInN b (Just t)   | t <= 0 = error "mhBurnInN: Auto tuning period smaller equal 0."-  | b > t = mhNIter t >> mcmcAutotune t >> mhBurnInN (b - t) (Just t)-  | otherwise = mhNIter b >> mcmcAutotune b+  | b > t = do+    mcmcResetA+    mhNIter t+    mcmcSummarizeCycle >>= mcmcDebugT+    mcmcAutotune+    mhBurnInN (b - t) (Just t)+  | otherwise = do+    mcmcResetA+    mhNIter b+    mcmcSummarizeCycle >>= mcmcInfoT+    mcmcInfoS $ "Acceptance ratios calculated over the last " <> show b <> " iterations." mhBurnInN b Nothing = mhNIter b  -- Initialize burn in for given number of iterations.@@ -110,48 +108,35 @@   | b < 0 = error "mhBurnIn: Negative number of burn in iterations."   | b == 0 = return ()   | otherwise = do-    liftIO $ putStrLn $ "-- Burn in for " <> show b <> " cycles."-    mcmcMonitorHeader+    mcmcInfoS $ "Burn in for " <> show b <> " cycles."+    mcmcDebugS $ "Auto tuning period is " <> show t <> "."+    mcmcMonitorStdOutHeader     mhBurnInN b t-    liftIO $ putStrLn "-- Burn in finished."-    case t of-      Nothing -> return ()-      Just _ -> mcmcSummarizeCycle t-    s <- get-    let a = acceptance s-    put $ s {acceptance = resetA a}+    mcmcInfoT "Burn in finished."  -- Run for given number of iterations. mhRun :: ToJSON a => Int -> Mcmc a () mhRun n = do-  liftIO $ putStrLn $ "-- Run chain for " <> show n <> " iterations."-  mcmcMonitorHeader+  mcmcInfoS $ "Run chain for " <> show n <> " iterations."+  mcmcMonitorStdOutHeader   mhNIter n  mhT :: ToJSON a => Mcmc a () mhT = do-  s <- get-  liftIO $ putStrLn "-- Start of Metropolis-Hastings sampler."-  mcmcInit+  mcmcInfoT "Metropolis-Hastings sampler."   mcmcReport-  mcmcSummarizeCycle Nothing+  mcmcSummarizeCycle >>= mcmcInfoT+  s <- get   let b = fromMaybe 0 (burnInIterations s)   mhBurnIn b (autoTuningPeriod s)-  let n = iterations s-  mhRun n-  mcmcClose+  mhRun $ iterations s  mhContinueT :: ToJSON a => Int -> Mcmc a () mhContinueT dn = do-  liftIO $ putStrLn "-- Continue Metropolis-Hastings sampler."-  liftIO $ putStrLn $ "-- Run chain for " <> show dn <> " additional iterations."-  s <- get-  let n = iterations s-  put s {iterations = n + dn}-  mcmcInit-  mcmcSummarizeCycle Nothing+  mcmcInfoT "Continuation of Metropolis-Hastings sampler."+  mcmcInfoS $ "Run chain for " <> show dn <> " additional iterations."+  mcmcSummarizeCycle >>= mcmcInfoT   mhRun dn-  mcmcClose  -- | Continue a Markov chain for a given number of Metropolis-Hastings steps. mhContinue ::@@ -161,9 +146,11 @@   -- | Loaded status of the Markov chain.   Status a ->   IO (Status a)-mhContinue dn+mhContinue dn s   | dn <= 0 = error "mhContinue: The number of iterations is zero or negative."-  | otherwise = execStateT $ mhContinueT dn+  | otherwise = mcmcRun (mhContinueT dn) s'+    where n' = iterations s + dn+          s' = s {iterations = n'}  -- | Run a Markov chain for a given number of Metropolis-Hastings steps. mh ::@@ -171,10 +158,9 @@   -- | Initial (or last) status of the Markov chain.   Status a ->   IO (Status a)-mh s = do-  let m = iteration s-  if m == 0-    then execStateT mhT s+mh s =+  if iteration s == 0+    then mcmcRun mhT s     else do       putStrLn "To continue a Markov chain run, please use 'mhContinue'."-      error $ "mh: Current iteration " ++ show m ++ " is non-zero."+      error $ "mh: Current iteration " ++ show (iteration s) ++ " is non-zero."

@@ -43,6 +43,7 @@ import Mcmc.Monitor.ParameterBatch import Mcmc.Monitor.Time import Mcmc.Trace+import Mcmc.Verbosity import Numeric.Log import System.Directory import System.IO@@ -62,7 +63,7 @@         mBatches :: [MonitorBatch a]       } --- | Monitor to standard output.+-- | Monitor to standard output; constructed with 'monitorStdOut'. data MonitorStdOut a   = MonitorStdOut       { msParams :: [MonitorParameter a],@@ -91,8 +92,8 @@   where     vals = map (T.justifyRight msWidth ' ') (tail xs) -msHeader :: MonitorStdOut a -> IO ()-msHeader m = T.hPutStr stdout $ T.unlines [row, sep]+msHeader :: MonitorStdOut a -> Text+msHeader m = T.intercalate "\n" [row, sep]   where     row =       msRenderRow $@@ -105,27 +106,32 @@ msExec ::   Int ->   Item a ->-  (Int, UTCTime) ->   Int ->+  UTCTime ->+  Int ->   MonitorStdOut a ->-  IO ()-msExec i (Item x p l) (ss, st) j m-  | i `mod` msPeriod m /= 0 =-    return ()+  IO (Maybe Text)+msExec i (Item x p l) ss st j m+  | i `mod` msPeriod m /= 0 = return Nothing   | otherwise = do       ct <- getCurrentTime       let dt = ct `diffUTCTime` st+          -- Careful, don't evaluate this when i == ss.           timePerIter = dt / fromIntegral (i - ss)-          eta = if i < 10-                then ""-                else renderDuration $ timePerIter * fromIntegral (j - i)-      T.hPutStrLn stdout-        $ msRenderRow-        $ [T.pack (show i), renderLog p, renderLog l, renderLog (p * l)]-          ++ [T.toLazyText $ mpFunc mp x | mp <- msParams m]-          ++ [renderDuration dt , eta]+          -- -- Always 0; doesn't make much sense.+          -- tpi = if (i - ss) < 10+          --       then ""+          --       else renderDurationS timePerIter+          eta =+            if (i - ss) < 10+            then ""+            else renderDuration $ timePerIter * fromIntegral (j - i)+      return $ Just $ msRenderRow $+        [T.pack (show i), renderLog p, renderLog l, renderLog (p * l)]+        ++ [T.toLazyText $ mpFunc mp x | mp <- msParams m]+        ++ [renderDuration dt, eta] --- | Monitor to a file.+-- | Monitor to a file; constructed with 'monitorFile'. data MonitorFile a   = MonitorFile       { mfName :: String,@@ -154,9 +160,18 @@ mfRenderRow :: [Text] -> Text mfRenderRow = T.intercalate "\t" -mfOpen :: String -> MonitorFile a -> IO (MonitorFile a)-mfOpen n m = do-  h <- openFile (n <> mfName m <> ".monitor") WriteMode+open' :: String -> Bool -> IO Handle+open' n frc = do+  fe <- doesFileExist n+  case (fe, frc) of+    (False, _) -> openFile n WriteMode+    (True, True) -> openFile n WriteMode+    (True, False) -> error $ "open': File \"" <> n <> "\" exists; probably use 'force'?"++mfOpen :: String -> Bool -> MonitorFile a -> IO (MonitorFile a)+mfOpen n frc m = do+  let mfn = n <> mfName m <> ".monitor"+  h <- open' mfn frc   hSetBuffering h LineBuffering   return $ m {mfHandle = Just h} @@ -165,9 +180,10 @@   let fn = n <> mfName m <> ".monitor"   fe <- doesFileExist fn   if fe-    then do h <- openFile fn AppendMode-            hSetBuffering h LineBuffering-            return $ m {mfHandle = Just h}+    then do+      h <- openFile fn AppendMode+      hSetBuffering h LineBuffering+      return $ m {mfHandle = Just h}     else error $ "mfAppend: Monitor file does not exist: " ++ fn ++ "."  mfHeader :: MonitorFile a -> IO ()@@ -210,10 +226,11 @@   Just h -> hClose h   Nothing -> error $ "mfClose: File was not opened for monitor " <> mfName m <> "." --- | Monitor to a file, but calculate batch means for the given batch size.+-- | Monitor to a file, but calculate batch means for the given batch size;+-- constructed with 'monitorBatch'. ----- XXX: Batch monitors are slow at the moment because the monitored parameter--- has to be extracted from the state for each iteration.+-- Batch monitors are slow at the moment because the monitored parameter has to+-- be extracted from the state for each iteration. data MonitorBatch a   = MonitorBatch       { mbName :: String,@@ -240,9 +257,10 @@   | p < 2 = error "monitorBatch: Batch size has to be 2 or larger."   | otherwise = MonitorBatch n Nothing ps p -mbOpen :: String -> MonitorBatch a -> IO (MonitorBatch a)-mbOpen n m = do-  h <- openFile (n <> mbName m <> ".batch") WriteMode+mbOpen :: String -> Bool -> MonitorBatch a -> IO (MonitorBatch a)+mbOpen n frc m = do+  let mfn = n <> mbName m <> ".batch"+  h <- open' mfn frc   hSetBuffering h LineBuffering   return $ m {mbHandle = Just h} @@ -251,9 +269,10 @@   let fn = n <> mbName m <> ".batch"   fe <- doesFileExist fn   if fe-    then do h <- openFile fn AppendMode-            hSetBuffering h LineBuffering-            return $ m {mbHandle = Just h}+    then do+      h <- openFile fn AppendMode+      hSetBuffering h LineBuffering+      return $ m {mbHandle = Just h}     else error $ "mbAppend: Monitor file does not exist: " ++ fn ++ "."  mbHeader :: MonitorBatch a -> IO ()@@ -308,11 +327,11 @@   Nothing -> error $ "mfClose: File was not opened for batch monitor: " <> mbName m <> "."  -- | Open the files associated with the 'Monitor'.-mOpen :: String -> Monitor a -> IO (Monitor a)-mOpen n (Monitor s fs bs) = do-  fs' <- mapM (mfOpen n) fs+mOpen :: String -> Bool -> Monitor a -> IO (Monitor a)+mOpen n frc (Monitor s fs bs) = do+  fs' <- mapM (mfOpen n frc) fs   mapM_ mfHeader fs'-  bs' <- mapM (mbOpen n) bs+  bs' <- mapM (mbOpen n frc) bs   mapM_ mbHeader bs'   return $ Monitor s fs' bs' @@ -323,27 +342,34 @@   bs' <- mapM (mbAppend n) bs   return $ Monitor s fs' bs' --- | Print header line of 'Monitor' (standard output only).-mHeader :: Monitor a -> IO ()+-- | Get header line of 'MonitorStdOut'.+mHeader :: Monitor a -> Text mHeader (Monitor s _ _) = msHeader s --- | Execute monitors; print status information to standard output and files.+-- | Execute monitors; print status information to files and return text to be+-- printed to standard output and log file. mExec ::+  -- | Verbosity+  Verbosity ->   -- | Iteration.   Int ->-  -- | Starting state and time.-  (Int, UTCTime) ->+  -- | Starting state.+  Int ->+  -- | Starting time.+  UTCTime ->   -- | Trace of Markov chain.   Trace a ->   -- | Total number of iterations; to calculate ETA.   Int ->   -- | The monitor.   Monitor a ->-  IO ()-mExec i t xs j (Monitor s fs bs) = do-  msExec i (headT xs) t j s+  IO (Maybe Text)+mExec v i ss st xs j (Monitor s fs bs) = do   mapM_ (mfExec i $ headT xs) fs   mapM_ (mbExec i xs) bs+  if v == Quiet+    then return Nothing+    else msExec i (headT xs) ss st j s  -- | Close the files associated with the 'Monitor'. mClose :: Monitor a -> IO (Monitor a)

@@ -12,7 +12,9 @@ -- -- Creation date: Fri May 29 11:11:49 2020. module Mcmc.Monitor.Parameter-  ( MonitorParameter (..),+  ( -- * Parameter monitors+    MonitorParameter (..),+    (@.),     monitorInt,     monitorRealFloat,     monitorRealFloatF,@@ -34,52 +36,50 @@         mpFunc :: a -> Builder       } +-- | Convert a parameter monitor from one data type to another using a lens.+--+-- For example, to monitor a real float value being the first entry of a tuple:+--+-- @+-- mon = _1 @@ monitorRealFloat+-- @+(@.) :: Lens' b a -> MonitorParameter a -> MonitorParameter b+(@.) l (MonitorParameter n f) = MonitorParameter n (\x -> f $ x^.l)+ -- | Monitor integral parameters such as 'Int'. monitorInt ::-  Integral b =>+  Integral a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameter a-monitorInt n l = MonitorParameter n (\x -> T.decimal $ x ^. l)-{-# SPECIALIZE monitorInt :: String -> Lens' a Int -> MonitorParameter a #-}+monitorInt n = MonitorParameter n T.decimal+{-# SPECIALIZE monitorInt :: String -> MonitorParameter Int #-}  -- | Monitor real float parameters such as 'Double' with eight decimal places -- (half precision). monitorRealFloat ::-  RealFloat b =>+  RealFloat a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameter a-monitorRealFloat n l =-  MonitorParameter n (\x -> T.formatRealFloat T.Fixed (Just 8) $ x ^. l)-{-# SPECIALIZE monitorRealFloat :: String -> Lens' a Double -> MonitorParameter a #-}+monitorRealFloat n = MonitorParameter n (T.formatRealFloat T.Fixed (Just 8))+{-# SPECIALIZE monitorRealFloat :: String -> MonitorParameter Double #-}  -- | Monitor real float parameters such as 'Double' with full precision. monitorRealFloatF ::-  RealFloat b =>+  RealFloat a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameter a-monitorRealFloatF n l = MonitorParameter n (\x -> T.realFloat $ x ^. l)-{-# SPECIALIZE monitorRealFloatF :: String -> Lens' a Double -> MonitorParameter a #-}+monitorRealFloatF n = MonitorParameter n T.realFloat+{-# SPECIALIZE monitorRealFloatF :: String -> MonitorParameter Double #-}  -- | Monitor real float parameters such as 'Double' with scientific notation and -- eight decimal places. monitorRealFloatS ::-  RealFloat b =>+  RealFloat a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameter a-monitorRealFloatS n l =-  MonitorParameter-    n-    (\x -> T.formatRealFloat T.Exponent (Just 8) $ x ^. l)-{-# SPECIALIZE monitorRealFloatS :: String -> Lens' a Double -> MonitorParameter a #-}+monitorRealFloatS n = MonitorParameter n (T.formatRealFloat T.Exponent (Just 8))+{-# SPECIALIZE monitorRealFloatS :: String -> MonitorParameter Double #-}

@@ -14,7 +14,9 @@ -- -- Batch mean monitors. module Mcmc.Monitor.ParameterBatch-  ( MonitorParameterBatch (..),+  ( -- * Batch parameter monitors+    MonitorParameterBatch (..),+    (@#),     monitorBatchMeanInt,     monitorBatchMeanIntF,     monitorBatchMeanRealFloat,@@ -43,8 +45,16 @@         mbpFunc :: [a] -> Builder       } -mapL :: Lens' a b -> [a] -> [b]-mapL l = map (^. l)+-- | Convert a batch parameter monitor from one data type to another using a+-- lens.+--+-- For example, to batch monitor a real float value being the first entry of a tuple:+--+-- @+-- mon = _1 @# monitorBatchMeanRealFloat+-- @+(@#) :: Lens' b a -> MonitorParameterBatch a -> MonitorParameterBatch b+(@#) l (MonitorParameterBatch n f) = MonitorParameterBatch n (f . map (^. l))  mean :: Real a => [a] -> Double mean xs = realToFrac (sum xs) / fromIntegral (length xs)@@ -54,83 +64,59 @@ -- | Batch monitor integral parameters such as 'Int'. Print the mean with eight -- decimal places (half precision). monitorBatchMeanInt ::-  Integral b =>+  Integral a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameterBatch a-monitorBatchMeanInt n l =-  MonitorParameterBatch-    n-    (T.formatRealFloat T.Fixed (Just 8) . mean . mapL l)-{-# SPECIALIZE monitorBatchMeanInt :: String -> Lens' a Int -> MonitorParameterBatch a #-}+monitorBatchMeanInt n = MonitorParameterBatch n (T.formatRealFloat T.Fixed (Just 8) . mean)+{-# SPECIALIZE monitorBatchMeanInt :: String -> MonitorParameterBatch Int #-}  -- | Batch monitor integral parameters such as 'Int'. Print the mean with full -- precision. monitorBatchMeanIntF ::-  Integral b =>+  Integral a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameterBatch a-monitorBatchMeanIntF n l =-  MonitorParameterBatch n (T.realFloat . mean . mapL l)-{-# SPECIALIZE monitorBatchMeanIntF :: String -> Lens' a Int -> MonitorParameterBatch a #-}+monitorBatchMeanIntF n = MonitorParameterBatch n (T.realFloat . mean)+{-# SPECIALIZE monitorBatchMeanIntF :: String -> MonitorParameterBatch Int #-}  -- | Batch monitor real float parameters such as 'Double' with eight decimal -- places (half precision). monitorBatchMeanRealFloat ::-  RealFloat b =>+  RealFloat a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameterBatch a-monitorBatchMeanRealFloat n l =-  MonitorParameterBatch-    n-    (T.formatRealFloat T.Fixed (Just 8) . mean . mapL l)-{-# SPECIALIZE monitorBatchMeanRealFloat :: String -> Lens' a Double -> MonitorParameterBatch a #-}+monitorBatchMeanRealFloat n = MonitorParameterBatch n (T.formatRealFloat T.Fixed (Just 8) . mean)+{-# SPECIALIZE monitorBatchMeanRealFloat :: String -> MonitorParameterBatch Double #-}  -- | Batch monitor real float parameters such as 'Double' with full precision. monitorBatchMeanRealFloatF ::-  RealFloat b =>+  RealFloat a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameterBatch a-monitorBatchMeanRealFloatF n l =-  MonitorParameterBatch n (T.realFloat . mean . mapL l)-{-# SPECIALIZE monitorBatchMeanRealFloatF :: String -> Lens' a Double -> MonitorParameterBatch a #-}+monitorBatchMeanRealFloatF n = MonitorParameterBatch n (T.realFloat . mean)+{-# SPECIALIZE monitorBatchMeanRealFloatF :: String -> MonitorParameterBatch Double #-}  -- | Batch monitor real float parameters such as 'Double' with scientific -- notation and eight decimal places. monitorBatchMeanRealFloatS ::-  RealFloat b =>+  RealFloat a =>   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   MonitorParameterBatch a-monitorBatchMeanRealFloatS n l =-  MonitorParameterBatch-    n-    (T.formatRealFloat T.Exponent (Just 8) . mean . mapL l)-{-# SPECIALIZE monitorBatchMeanRealFloatS :: String -> Lens' a Double -> MonitorParameterBatch a #-}+monitorBatchMeanRealFloatS n = MonitorParameterBatch n (T.formatRealFloat T.Exponent (Just 8) . mean)+{-# SPECIALIZE monitorBatchMeanRealFloatS :: String -> MonitorParameterBatch Double #-}  -- | Batch monitor parameters with custom lens and builder. monitorBatchCustom ::   -- | Name of monitor.   String ->-  -- | Instruction about which parameter to monitor.-  Lens' a b ->   -- | Function to calculate the batch mean.-  ([b] -> b) ->+  ([a] -> a) ->   -- | Custom builder.-  (b -> Builder) ->+  (a -> Builder) ->   MonitorParameterBatch a-monitorBatchCustom n l f b =-  MonitorParameterBatch n (b . f . mapL l)+monitorBatchCustom n f b = MonitorParameterBatch n (b . f)

@@ -13,6 +13,7 @@ -- Creation date: Fri May 29 12:36:43 2020. module Mcmc.Monitor.Time   ( renderDuration,+    renderDurationS,   ) where @@ -36,3 +37,10 @@     ts :: Int     ts = round dt     renderDecimal n = T.justifyRight 2 '0' $ T.toLazyText $ T.decimal n++-- | Render duration in seconds.+renderDurationS :: NominalDiffTime -> Text+renderDurationS dt = T.toLazyText $ T.decimal ts+  where+    ts :: Int+    ts = round dt

@@ -1,285 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE RankNTypes #-}---- TODO: Allow execution of moves in order of appearance in the cycle.---- TODO: Moves and monitors both use lenses and names for what they move and--- monitor. Should a data structure be used combining the lens and the name, so--- that things are cohesive?---- TODO: Moves on simplices: SimplexElementScale (?).---- TODO: Moves on tree branch lengths.--- - Slide a node on the tree.--- - Scale a tree.---- TODO: Moves on tree topologies.--- - NNI--- - Narrow (what is this, see RevBayes)--- - FNPR (dito)---- TODO: Bactrian moves; https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/.------ slideBactrian------ scaleBactrian---- |--- Module      :  Mcmc.Move--- Description :  Moves and cycles--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 20 13:42:53 2020.-module Mcmc.Move-  ( -- * Move-    Move (..),-    MoveSimple (..),-    Tuner (tParam, tFunc),-    tuner,-    tune,-    autotune,--    -- * Cycle-    Cycle (fromCycle),-    fromList,-    autotuneC,-    summarizeCycle,--    -- * Acceptance-    Acceptance (..),-    emptyA,-    pushA,-    resetA,-    acceptanceRatios,-  )-where--import Data.Aeson-import Data.Function-import Data.List-import qualified Data.Map.Strict as M-import Data.Map.Strict (Map)-import Data.Maybe-import qualified Data.Text.Lazy as T-import Data.Text.Lazy (Text)-import qualified Data.Text.Lazy.Builder as B-import qualified Data.Text.Lazy.Builder.Int as B-import qualified Data.Text.Lazy.Builder.RealFloat as B-import Numeric.Log hiding (sum)-import System.Random.MWC---- | A 'Move' is an instruction about how the Markov chain will traverse the--- state space @a@. Essentially, it is a probability density conditioned on the--- current state.------ A 'Move' may be tuneable in that it contains information about how to enlarge--- or shrink the step size to tune the acceptance ratio.-data Move a = Move-  { -- | Name (no moves with the same name are allowed in a 'Cycle').-    mvName :: String,-    -- | The weight determines how often a 'Move' is executed per iteration of-    -- the Markov chain.-    mvWeight :: Int,-    -- | Simple move without tuning information.-    mvSimple :: MoveSimple a,-    -- | Tuning is disabled if set to 'Nothing'.-    mvTune :: Maybe (Tuner a)-  }--instance Show (Move a) where-  show m = show $ mvName m--instance Eq (Move a) where-  m == n = mvName m == mvName n--instance Ord (Move a) where-  compare = compare `on` mvName---- XXX: One could also use a different type for 'mvSample', so that--- 'mvDensity' can be avoided. In detail,------ @---   mvSample :: a -> GenIO -> IO (a, Log Double, Log, Double)--- @------ where the densities describe the probability of going there and back.--- However, we may need more information about the move for other MCMC samplers--- different from Metropolis-Hastings.---- | Simple move without tuning information.------ In order to calculate the Metropolis-Hastings ratio, we need to know the--- probability (density) of jumping forth, and the probability (density) of--- jumping back.-data MoveSimple a = MoveSimple-  { -- | Instruction about randomly moving from the current state to a new-    -- state, given some source of randomness.-    mvSample :: a -> GenIO -> IO a,-    -- | The density of going from one state to another. Set to 'Nothing' for-    -- symmetric moves.-    mvDensity :: Maybe (a -> a -> Log Double)-  }---- | Tune the acceptance ratio of a 'Move'; see 'tune', or 'autotune'.-data Tuner a = Tuner-  { tParam :: Double,-    tFunc :: Double -> MoveSimple a-  }---- | Create a 'Tuner'. The tuning function accepts a tuning parameter, and--- returns a corresponding 'MoveSimple'. The larger the tuning parameter, the--- larger the 'Move', and vice versa.-tuner :: (Double -> MoveSimple a) -> Tuner a-tuner = Tuner 1.0---- | Tune a 'Move'. Return 'Nothing' if 'Move' is not tuneable. If the parameter---   @dt@ is larger than 1.0, the 'Move' is enlarged, if @0<dt<1.0@, it is---   shrunk. Negative tuning parameters are not allowed.-tune :: Double -> Move a -> Maybe (Move a)-tune dt tm-  | dt <= 0 = error $ "tune: Tuning parameter not positive: " <> show dt <> "."-  | otherwise = do-    (Tuner t f) <- mvTune tm-    let t' = t * dt-    return $ tm {mvSimple = f t', mvTune = Just $ Tuner t' f}--ratioOpt :: Double-ratioOpt = 0.44---- | For a given acceptance ratio, auto tune the 'Move'. For now, a 'Move' is--- enlarged when the acceptance ratio is above 0.44, and shrunk otherwise.--- Return 'Nothing' if 'Move' is not tuneable.------ XXX: The desired acceptance ratio 0.44 is optimal for one-dimensional--- 'Move's; one could also store the affected number of dimensions with the--- 'Move' and tune towards an acceptance ratio accounting for the number of--- dimensions.-autotune :: Double -> Move a -> Maybe (Move a)-autotune a = tune (a / ratioOpt)---- | In brief, a 'Cycle' is a list of moves. The state of the Markov chain will--- be logged only after each 'Cycle', and the iteration counter will be--- increased by one. __Moves must have unique names__, so that they can be--- identified.------ 'Move's are replicated according to their weights and executed in random--- order. That is, they are not executed in the order they appear in the--- 'Cycle'. However, if a 'Move' has weight @w@, it is executed exactly @w@--- times per iteration.-newtype Cycle a = Cycle {fromCycle :: [Move a]}---- | Create a 'Cycle' from a list of 'Move's.-fromList :: [Move a] -> Cycle a-fromList [] =-  error "fromList: Received an empty list but cannot create an empty Cycle."-fromList xs =-  if length (nub nms) == length nms-    then Cycle xs-    else error "fromList: Moves don't have unique names."-  where-    nms = map mvName xs---- | Tune the 'Move's in the 'Cycle'. Tuning has no effect on 'Move's that--- cannot be tuned. See 'autotune'.-autotuneC :: Int -> Acceptance (Move a) -> Cycle a -> Cycle a-autotuneC n a = Cycle . map tuneF . fromCycle-  where-    ars = acceptanceRatios n a-    tuneF m = fromMaybe m (autotune (ars M.! m) m)--renderRow :: Text -> Text -> Text -> Text -> Text-renderRow name weight acceptRatio tuneParam = "   " <> nB <> wB <> rB <> tB-  where-    nB = T.justifyLeft 30 ' ' name-    wB = T.justifyRight 8 ' ' weight-    rB = T.justifyRight 20 ' ' acceptRatio-    tB = T.justifyRight 20 ' ' tuneParam--moveHeader :: Text-moveHeader =-  renderRow "Move name" "Weight" "Acceptance ratio" "Tuning parameter"--summarizeMove :: Move a -> Maybe Double -> Text-summarizeMove m r = renderRow (T.pack name) weight acceptRatio tuneParamStr-  where-    name = mvName m-    weight = B.toLazyText $ B.decimal $ mvWeight m-    acceptRatio = B.toLazyText $ maybe "" (B.formatRealFloat B.Fixed (Just 3)) r-    tuneParamStr = B.toLazyText $ maybe "" (B.formatRealFloat B.Fixed (Just 3)) (tParam <$> mvTune m)---- | Summarize the 'Move's in the 'Cycle'. Also report acceptance ratios for the--- given number of last iterations.-summarizeCycle :: Maybe (Int, Acceptance (Move a)) -> Cycle a -> Text-summarizeCycle acc c =-  T.unlines $-    [ "-- Summary of move(s) in cycle.",-      -- T.replicate (T.length moveHeader) "─",-      moveHeader,-      "   " <> T.replicate (T.length moveHeader - 3) "─"-    ]-      ++ [summarizeMove m (ar m) | m <- mvs]-      ++ [ "   " <> T.replicate (T.length moveHeader - 3) "─",-           B.toLazyText $-             B.fromLazyText "-- "-               <> B.decimal mpi-               <> B.fromString " move(s) per iteration."-               <> arStr-         ]-  where-    mvs = fromCycle c-    mpi = sum $ map mvWeight mvs-    arStr = case acc of-      Nothing -> ""-      Just (n, _) ->-        " Acceptance ratio(s) calculated over " <> B.decimal n <> " iterations."-    ars = case acc of-      Nothing -> M.empty-      Just (n, a) -> acceptanceRatios n a-    ar m = ars M.!? m---- | For each key @k@, store the list of accepted (True) and rejected (False)--- proposals. For reasons of efficiency, the lists are stored in reverse order;--- latest first.-newtype Acceptance k = Acceptance {fromAcceptance :: Map k [Bool]}--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, []) | k <- ks]---- | For key @k@, prepend an accepted (True) or rejected (False) proposal.-pushA :: (Ord k, Show k) => k -> Bool -> Acceptance k -> Acceptance k--- Unsafe; faster.-pushA k v (Acceptance m) = Acceptance $ M.adjust (v :) k m--- -- Safe; slower.--- prependA k v (Acceptance m) | k `M.member` m = Acceptance $ M.adjust (v:) k m---                             | otherwise = error msg---   where msg = "prependA: Can not add acceptance value for key: " <> show k <> "."-{-# INLINEABLE pushA #-}---- | Reset acceptance storage.-resetA :: Acceptance k -> Acceptance k-resetA = Acceptance . M.map (const []) . fromAcceptance--ratio :: Int -> [Bool] -> Double-ratio n xs = fromIntegral (length ts) / fromIntegral n-  where-    ts = filter (== True) xs---- | Acceptance ratios averaged over the given number of last iterations. If--- less than @n@ iterations are available, only those are used.-acceptanceRatios :: Int -> Acceptance k -> Map k Double-acceptanceRatios n (Acceptance m) = M.map (ratio n . take n) m

@@ -1,136 +0,0 @@-{-# LANGUAGE RankNTypes #-}---- Technically, only a Getter is needed when calculating the density of the move--- ('densCont', and similar functions). I tried splitting the lens into a getter--- and a setter. However, speed improvements were marginal, and some times not--- even measurable. Using a 'Lens'' is just easier, and has no real drawbacks.---- |--- Module      :  Mcmc.Move.Generic--- Description :  Generic interface to create moves--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu May 14 20:26:27 2020.-module Mcmc.Move.Generic-  ( moveGenericContinuous,-    moveSymmetricGenericContinuous,-    moveGenericDiscrete,-    moveSymmetricGenericDiscrete,-  )-where--import Lens.Micro-import Mcmc.Move-import Numeric.Log-import Statistics.Distribution-import System.Random.MWC--jumpCont ::-  (ContDistr d, ContGen d) =>-  Lens' a Double ->-  d ->-  (Double -> Double -> Double) ->-  a ->-  GenIO ->-  IO a-jumpCont l d f x g = do-  dx <- genContVar d g-  return $ set l ((x ^. l) `f` dx) x-{-# INLINEABLE jumpCont #-}--densCont ::-  (ContDistr d, ContGen d) =>-  Lens' a Double ->-  d ->-  (Double -> Double -> Double) ->-  a ->-  a ->-  Log Double-densCont l d fInv x y = Exp $ logDensity d ((y ^. l) `fInv` (x ^. l))-{-# INLINEABLE densCont #-}---- | Generic function to create moves for continuous parameters ('Double').-moveGenericContinuous ::-  (ContDistr d, ContGen d) =>-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Probability distribution-  d ->-  -- | Forward operator, e.g. (+), so that x + dx = y.-  (Double -> Double -> Double) ->-  -- | Inverse operator, e.g.,(-), so that y - dx = x.-  (Double -> Double -> Double) ->-  MoveSimple a-moveGenericContinuous l d f fInv =-  MoveSimple (jumpCont l d f) (Just $ densCont l d fInv)---- | Generic function to create symmetric moves for continuous parameters ('Double').-moveSymmetricGenericContinuous ::-  (ContDistr d, ContGen d) =>-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Probability distribution-  d ->-  -- | Forward operator, e.g. (+), so that x + dx = y.-  (Double -> Double -> Double) ->-  MoveSimple a-moveSymmetricGenericContinuous l d f =-  MoveSimple (jumpCont l d f) Nothing--jumpDiscrete ::-  (DiscreteDistr d, DiscreteGen d) =>-  Lens' a Int ->-  d ->-  (Int -> Int -> Int) ->-  a ->-  GenIO ->-  IO a-jumpDiscrete l d f x g = do-  dx <- genDiscreteVar d g-  return $ set l ((x ^. l) `f` dx) x-{-# INLINEABLE jumpDiscrete #-}--densDiscrete ::-  (DiscreteDistr d, DiscreteGen d) =>-  Lens' a Int ->-  d ->-  (Int -> Int -> Int) ->-  a ->-  a ->-  Log Double-densDiscrete l d fInv x y =-  Exp $ logProbability d ((y ^. l) `fInv` (x ^. l))-{-# INLINEABLE densDiscrete #-}---- | Generic function to create moves for discrete parameters ('Int').-moveGenericDiscrete ::-  (DiscreteDistr d, DiscreteGen d) =>-  -- | Instruction about which parameter to change.-  Lens' a Int ->-  -- | Probability distribution.-  d ->-  -- | Forward operator, e.g. (+), so that x + dx = y.-  (Int -> Int -> Int) ->-  -- | Inverse operator, e.g.,(-), so that y - dx = x.-  (Int -> Int -> Int) ->-  MoveSimple a-moveGenericDiscrete l fd f fInv =-  MoveSimple (jumpDiscrete l fd f) (Just $ densDiscrete l fd fInv)---- | Generic function to create symmetric moves for discrete parameters ('Int').-moveSymmetricGenericDiscrete ::-  (DiscreteDistr d, DiscreteGen d) =>-  -- | Instruction about which parameter to change.-  Lens' a Int ->-  -- | Probability distribution.-  d ->-  -- | Forward operator, e.g. (+), so that x + dx = y.-  (Int -> Int -> Int) ->-  MoveSimple a-moveSymmetricGenericDiscrete l fd f =-  MoveSimple (jumpDiscrete l fd f) Nothing

@@ -1,69 +0,0 @@-{-# LANGUAGE RankNTypes #-}---- |--- Module      :  Mcmc.Move.Scale--- Description :  Scaling move with Gamma distribution--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu May 14 21:49:23 2020.-module Mcmc.Move.Scale-  ( scale,-    scaleUnbiased,-  )-where--import Lens.Micro-import Mcmc.Move-import Mcmc.Move.Generic-import Statistics.Distribution.Gamma---- The actual move with tuning parameter.-scaleSimple :: Lens' a Double -> Double -> Double -> Double -> MoveSimple a-scaleSimple l k th t = moveGenericContinuous l (gammaDistr k (t * th)) (*) (/)---- | Multiplicative move with Gamma distributed density.-scale ::-  -- | Name.-  String ->-  -- | Weight.-  Int ->-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Shape.-  Double ->-  -- | Scale.-  Double ->-  -- | Enable tuning.-  Bool ->-  Move a-scale n w l k th True =-  Move n w (scaleSimple l k th 1.0) (Just $ tuner $ scaleSimple l k th)-scale n w l k th False = Move n w (scaleSimple l k th 1.0) Nothing---- | Multiplicative move with Gamma distributed density. The scale of the Gamma--- distributions is set to (shape)^{-1}, so that the mean of the Gamma--- distribution is 1.0.-scaleUnbiased ::-  -- | Name.-  String ->-  -- | Weight.-  Int ->-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Shape.-  Double ->-  -- | Enable tuning.-  Bool ->-  Move a-scaleUnbiased n w l k True =-  Move-    n-    w-    (scaleSimple l k (1 / k) 1.0)-    (Just $ tuner $ scaleSimple l k (1 / k))-scaleUnbiased n w l k False = Move n w (scaleSimple l k (1 / k) 1.0) Nothing

@@ -1,95 +0,0 @@-{-# LANGUAGE RankNTypes #-}---- |--- Module      :  Mcmc.Move.Slide--- Description :  Normally distributed move--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May  6 10:59:13 2020.-module Mcmc.Move.Slide-  ( slide,-    slideSymmetric,-    slideUniform,-  )-where--import Lens.Micro-import Mcmc.Move-import Mcmc.Move.Generic-import Statistics.Distribution.Normal-import Statistics.Distribution.Uniform---- The actual move with tuning parameter.-slideSimple :: Lens' a Double -> Double -> Double -> Double -> MoveSimple a-slideSimple l m s t = moveGenericContinuous l (normalDistr m (s * t)) (+) (-)---- | Additive move with normally distributed density.-slide ::-  -- | Name.-  String ->-  -- | Weight.-  Int ->-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Mean.-  Double ->-  -- | Standard deviation.-  Double ->-  -- | Enable tuning.-  Bool ->-  Move a-slide n w l m s True =-  Move n w (slideSimple l m s 1.0) (Just $ tuner (slideSimple l m s))-slide n w l m s False = Move n w (slideSimple l m s 1.0) Nothing---- The actual move with tuning parameter.-slideSymmetricSimple :: Lens' a Double -> Double -> Double -> MoveSimple a-slideSymmetricSimple l s t = moveSymmetricGenericContinuous l (normalDistr 0.0 (s * t)) (+)---- | Additive move with normally distributed density with mean zero. This move--- is very fast, because the Metropolis-Hastings ratio does not include--- calculation of the forwards and backwards densities.-slideSymmetric ::-  -- | Name.-  String ->-  -- | Weight.-  Int ->-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Standard deviation.-  Double ->-  -- | Enable tuning.-  Bool ->-  Move a-slideSymmetric n w l s True =-  Move n w (slideSymmetricSimple l s 1.0) (Just $ tuner (slideSymmetricSimple l s))-slideSymmetric n w l s False = Move n w (slideSymmetricSimple l s 1.0) Nothing---- The actual move with tuning parameter.-slideUniformSimple :: Lens' a Double -> Double -> Double -> MoveSimple a-slideUniformSimple l d t =-  moveSymmetricGenericContinuous l (uniformDistr (- t * d) (t * d)) (+)---- | Additive move with uniformly distributed density. This move is very fast,--- because the Metropolis-Hastings ratio does not include calculation of the--- forwards and backwards densities.-slideUniform ::-  -- | Name.-  String ->-  -- | Weight.-  Int ->-  -- | Instruction about which parameter to change.-  Lens' a Double ->-  -- | Delta.-  Double ->-  -- | Enable tuning.-  Bool ->-  Move a-slideUniform n w l d True =-  Move n w (slideUniformSimple l d 1.0) (Just $ tuner (slideUniformSimple l d))-slideUniform n w l d False = Move n w (slideUniformSimple l d 1.0) Nothing

@@ -0,0 +1,344 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE RankNTypes #-}++-- TODO: Proposals on simplices: SimplexElementScale (?).++-- TODO: Proposals on tree topologies.+-- - NNI+-- - Narrow (what is this, see RevBayes)+-- - FNPR (dito)++-- |+-- Module      :  Mcmc.Proposal+-- Description :  Proposals and cycles+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 20 13:42:53 2020.+module Mcmc.Proposal+  ( -- * Proposal+    Proposal (..),+    (@~),+    ProposalSimple (..),+    Tuner (tParam, tFunc),+    createProposal,+    tune,++    -- * Cycle+    Order (..),+    Cycle (ccProposals),+    fromList,+    setOrder,+    getNCycles,+    tuneCycle,+    autotuneCycle,+    summarizeCycle,++    -- * Acceptance+    Acceptance (fromAcceptance),+    emptyA,+    pushA,+    resetA,+    transformKeysA,+    acceptanceRatios,+  )+where++import Data.Aeson+import Data.Default+import Data.Function+import Data.List+import qualified Data.Map.Strict as M+import Data.Map.Strict (Map)+import Data.Maybe+import qualified Data.Text.Lazy as T+import Data.Text.Lazy (Text)+import qualified Data.Text.Lazy.Builder as B+import qualified Data.Text.Lazy.Builder.Int as B+import qualified Data.Text.Lazy.Builder.RealFloat as B+import Lens.Micro+import Mcmc.Tools.Shuffle+import Numeric.Log hiding (sum)+import System.Random.MWC++-- | A 'Proposal' is an instruction about how the Markov chain will traverse the+-- state space @a@. Essentially, it is a probability mass or probability density+-- conditioned on the current state (i.e., a kernel).+--+-- A 'Proposal' may be tuneable in that it contains information about how to enlarge+-- or shrink the step size to tune the acceptance ratio.+data Proposal a = Proposal+  { -- | Name (no proposals with the same name are allowed in a 'Cycle').+    pName :: String,+    -- | The weight determines how often a 'Proposal' is executed per iteration of+    -- the Markov chain.+    pWeight :: Int,+    -- | Simple proposal without tuning information.+    pSimple :: ProposalSimple a,+    -- | Tuning is disabled if set to 'Nothing'.+    pTuner :: Maybe (Tuner a)+  }++instance Show (Proposal a) where+  show m = show $ pName m++instance Eq (Proposal a) where+  m == n = pName m == pName n++instance Ord (Proposal a) where+  compare = compare `on` pName++convertP :: Lens' b a -> Proposal a -> Proposal b+convertP l (Proposal n w s t) = Proposal n w (convertS l s) (convertT l <$> t)++-- | Convert a proposal from one data type to another using a lens.+--+-- For example:+--+-- @+-- scaleFirstEntryOfTuple = scale >>> _1+-- @+(@~) :: Lens' b a -> Proposal a -> Proposal b+(@~) = convertP++-- | 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 ratio, we need to know the+-- ratio of the backward to forward kernels (i.e., the probability masses or+-- probability densities). For unbiased proposals, this ratio is 1.0.+newtype ProposalSimple a = ProposalSimple+  {     pSample :: a -> GenIO -> IO (a, Log Double)+  }++convertS :: Lens' b a -> ProposalSimple a -> ProposalSimple b+convertS l (ProposalSimple s) = ProposalSimple s'+  where+    s' v g = do+      (x', r) <- s (v ^. l) g+      return (set l x' v, r)++-- | Tune the acceptance ratio of a 'Proposal'; see 'tune', or 'autotuneCycle'.+data Tuner a = Tuner+  { tParam :: Double,+    tFunc :: Double -> ProposalSimple a+  }++convertT :: Lens' b a -> Tuner a -> Tuner b+convertT l (Tuner p f) = Tuner p f'+  where+    f' x = convertS l $ f x++-- | Create a possibly tuneable proposal.+createProposal ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Function creating a simple proposal for a given tuning parameter. The+  -- larger the tuning parameter, the larger the proposal (and the lower the+  -- expected acceptance ratio), and vice versa.+  (Double -> ProposalSimple a) ->+  -- | Activate tuning?+  Bool ->+  Proposal a+createProposal n w f True = Proposal n w (f 1.0) (Just $ Tuner 1.0 f)+createProposal n w f False = Proposal n w (f 1.0) Nothing++-- Minimal tuning parameter; subject to change.+tuningParamMin :: Double+tuningParamMin = 1e-12++-- | Tune a 'Proposal'. Return 'Nothing' if 'Proposal' is not tuneable. If the parameter+--   @dt@ is larger than 1.0, the 'Proposal' is enlarged, if @0<dt<1.0@, it is+--   shrunk. Negative tuning parameters are not allowed.+tune :: Double -> Proposal a -> Maybe (Proposal a)+tune dt m+  | dt <= 0 = error $ "tune: Tuning parameter not positive: " <> show dt <> "."+  | otherwise = do+    (Tuner t f) <- pTuner m+    -- Ensure that the tuning parameter is not too small.+    let t' = max tuningParamMin (t * dt)+    return $ m {pSimple = f t', pTuner = Just $ Tuner t' f}++-- XXX: The desired acceptance ratio 0.44 is optimal for one-dimensional+-- 'Proposal's; one could also store the affected number of dimensions with the+-- 'Proposal' and tune towards an acceptance ratio accounting for the number of+-- dimensions.+ratioOpt :: Double+ratioOpt = 0.44++-- | Define the order in which 'Proposal's are executed in a 'Cycle'. The total+-- number of 'Proposal's per 'Cycle' may differ between 'Order's (e.g., compare+-- 'RandomO' and 'RandomReversibleO').+data Order+  = -- | Shuffle the 'Proposal's in the 'Cycle'. The 'Proposal's are replicated+    -- according to their weights and executed in random order. If a 'Proposal' has+    -- weight @w@, it is executed exactly @w@ times per iteration.+    RandomO+  | -- | The 'Proposal's are executed sequentially, in the order they appear in the+    -- 'Cycle'. 'Proposal's with weight @w>1@ are repeated immediately @w@ times+    -- (and not appended to the end of the list).+    SequentialO+  | -- | Similar to 'RandomO'. However, a reversed copy of the list of+    --  shuffled 'Proposal's is appended such that the resulting Markov chain is+    --  reversible.+    --  Note: the total number of 'Proposal's executed per cycle is twice the number+    --  of 'RandomO'.+    RandomReversibleO+  | -- | Similar to 'SequentialO'. However, a reversed copy of the list of+    -- sequentially ordered 'Proposal's is appended such that the resulting Markov+    -- chain is reversible.+    SequentialReversibleO+  deriving (Eq, Show)++instance Default Order where def = RandomO++-- | In brief, a 'Cycle' is a list of proposals. The state of the Markov chain will+-- be logged only after all 'Proposal's in the 'Cycle' have been completed, and the+-- iteration counter will be increased by one. The order in which the 'Proposal's+-- are executed is specified by 'Order'. The default is 'RandomO'.+--+-- __Proposals must have unique names__, so that they can be identified.+data Cycle a = Cycle+  { ccProposals :: [Proposal a],+    ccOrder :: Order+  }++-- | Create a 'Cycle' from a list of 'Proposal's.+fromList :: [Proposal a] -> Cycle a+fromList [] =+  error "fromList: Received an empty list but cannot create an empty Cycle."+fromList xs =+  if length (nub nms) == length nms+    then Cycle xs def+    else error "fromList: Proposals don't have unique names."+  where+    nms = map pName xs++-- | Set the order of 'Proposal's in a 'Cycle'.+setOrder :: Order -> Cycle a -> Cycle a+setOrder o c = c {ccOrder = o}++-- | Replicate 'Proposal's according to their weights and possibly shuffle them.+getNCycles :: Cycle a -> Int -> GenIO -> IO [[Proposal a]]+getNCycles (Cycle xs o) n g = case o of+  RandomO -> shuffleN ps n g+  SequentialO -> return $ replicate n ps+  RandomReversibleO -> do+    psRs <- shuffleN ps n g+    return [psR ++ reverse psR | psR <- psRs]+  SequentialReversibleO -> return $ replicate n $ ps ++ reverse ps+  where+    !ps = concat [replicate (pWeight m) m | m <- xs]++-- | Tune 'Proposal's in the 'Cycle'. See 'tune'.+tuneCycle :: Map (Proposal a) Double -> Cycle a -> Cycle a+tuneCycle m c =+  if sort (M.keys m) == sort ps+    then c {ccProposals = map tuneF ps}+    else error "tuneCycle: Map contains proposals that are not in the cycle."+  where+    ps = ccProposals c+    tuneF p = case m M.!? p of+      Nothing -> p+      Just x -> fromMaybe p (tune x p)++-- | Calculate acceptance ratios and auto tune the 'Proposal's in the 'Cycle'. For+-- now, a 'Proposal' is enlarged when the acceptance ratio is above 0.44, and+-- shrunk otherwise. Do not change 'Proposal's that are not tuneable.+autotuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a+autotuneCycle a = tuneCycle (M.map (\x -> exp $ x - ratioOpt) $ acceptanceRatios a)++renderRow :: Text -> Text -> Text -> Text -> Text -> Text -> Text+renderRow name weight nAccept nReject acceptRatio tuneParam = "   " <> nm <> wt <> na <> nr <> ra <> tp+  where+    nm = T.justifyLeft 30 ' ' name+    wt = T.justifyRight 8 ' ' weight+    na = T.justifyRight 15 ' ' nAccept+    nr = T.justifyRight 15 ' ' nReject+    ra = T.justifyRight 15 ' ' acceptRatio+    tp = T.justifyRight 20 ' ' tuneParam++proposalHeader :: Text+proposalHeader =+  renderRow "Proposal" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter"++summarizeProposal :: Proposal a -> Maybe (Int, Int, Double) -> Text+summarizeProposal m r = renderRow (T.pack name) weight nAccept nReject acceptRatio tuneParamStr+  where+    name = pName m+    weight = B.toLazyText $ B.decimal $ pWeight m+    nAccept = B.toLazyText $ maybe "" (B.decimal . (^. _1)) r+    nReject = B.toLazyText $ maybe "" (B.decimal . (^. _2)) r+    acceptRatio = B.toLazyText $ maybe "" (B.formatRealFloat B.Fixed (Just 3) . (^. _3)) r+    tuneParamStr = B.toLazyText $ maybe "" (B.formatRealFloat B.Fixed (Just 3)) (tParam <$> pTuner m)++-- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance ratios.+summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> Text+summarizeCycle a c =+  T.intercalate "\n" $+    [ "Summary of proposal(s) in cycle. " <> mpi <> " proposal(s) per iteration.",+      proposalHeader,+      "   " <> T.replicate (T.length proposalHeader - 3) "─"+    ]+      ++ [summarizeProposal m (ar m) | m <- ps]+      ++ ["   " <> T.replicate (T.length proposalHeader - 3) "─"]+  where+    ps = ccProposals c+    mpi = B.toLazyText $ B.decimal $ sum $ map pWeight ps+    ar m = acceptanceRatio m a++-- | For each key @k@, store the number of accepted and rejected proposals.+newtype Acceptance k = Acceptance {fromAcceptance :: Map k (Int, Int)}++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, Show k) => k -> Bool -> Acceptance k -> Acceptance k+pushA k True = Acceptance . M.adjust (\(a, r) -> (succ a, r)) k . fromAcceptance+pushA k False = Acceptance . M.adjust (\(a, r) -> (a, succ r)) 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] -> Map k1 v -> 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 ratio for a specific proposal.+acceptanceRatio :: (Show k, Ord k) => k -> Acceptance k -> Maybe (Int, Int, Double)+acceptanceRatio 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 $ "acceptanceRatio: Key not found in map: " ++ show k ++ "."++-- | Acceptance ratios for all proposals.+acceptanceRatios :: Acceptance k -> Map k Double+acceptanceRatios = M.map (\(as, rs) -> fromIntegral as / fromIntegral (as + rs)) . fromAcceptance

@@ -0,0 +1,137 @@+-- |+-- Module      :  Mcmc.Proposal.Bactrian+-- Description :  Bactrian proposals+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jun 25 15:49:48 2020.+--+-- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/.+module Mcmc.Proposal.Bactrian+  ( slideBactrian,+    scaleBactrian,+  )+where++import Mcmc.Proposal+import Numeric.Log+import Statistics.Distribution+import Statistics.Distribution.Normal+import System.Random.MWC+import System.Random.MWC.Distributions++genBactrian ::+  Double ->+  Double ->+  GenIO ->+  IO Double+genBactrian m s g = do+  let mn = m * s+      sd = sqrt (1 - m * m) * s+      d = normalDistr mn sd+  x <- genContVar d g+  b <- bernoulli 0.5 g+  return $ if b then x else - x++logDensityBactrian :: Double -> Double -> Double -> Log Double+logDensityBactrian m s x = Exp $ log $ kernel1 + kernel2+  where+    mn = m * s+    sd = sqrt (1 - m * m) * s+    dist1 = normalDistr (- mn) sd+    dist2 = normalDistr mn sd+    kernel1 = density dist1 x+    kernel2 = density dist2 x++bactrianAdditive ::+  Double ->+  Double ->+  Double ->+  GenIO ->+  IO (Double, Log Double)+bactrianAdditive m s x g = do+  dx <- genBactrian m s g+  return (x + dx, 1.0)++-- bactrianSimple lens spike stdDev tune forwardOp backwardOp+bactrianAdditiveSimple ::+  Double ->+  Double ->+  Double ->+  ProposalSimple Double+bactrianAdditiveSimple m s t+  | m < 0 = error "bactrianAdditiveSimple: Spike parameter negative."+  | m >= 1 = error "bactrianAdditiveSimple: Spike parameter 1.0 or larger."+  | s <= 0 = error "bactrianAdditiveSimple: Standard deviation 0.0 or smaller."+  | otherwise = ProposalSimple $ bactrianAdditive m (t * s)++-- | Additive symmetric proposal with kernel similar to the silhouette of a+-- Bactrian camel. The [Bactrian+-- kernel](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/figure/fig01)+-- is a mixture of two symmetrically arranged normal distributions. The spike+-- parameter loosely determines the standard deviations of the individual humps+-- while the other parameter refers to the standard deviation of the complete+-- Bactrian kernel.+--+-- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/.+slideBactrian ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Spike parameter.+  Double ->+  -- | Standard deviation.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+slideBactrian n w m s = createProposal n w (bactrianAdditiveSimple m s)++-- We have:+-- x  (1+dx ) = x'+-- x' (1+dx') = x.+--+-- Hence,+-- dx' = 1/(1-dx) - 1.+fInv :: Double -> Double+fInv dx = recip (1-dx) - 1++bactrianMult ::+  Double ->+  Double ->+  Double ->+  GenIO ->+  IO (Double, Log Double)+bactrianMult m s x g = do+  dx <- genBactrian m s g+  let qXY = logDensityBactrian m s dx+      qYX = logDensityBactrian m s (fInv dx)+  return (x * (1 + dx), qYX / qXY)++bactrianMultSimple :: Double -> Double -> Double -> ProposalSimple Double+bactrianMultSimple m s t+  | m < 0 = error "bactrianMultSimple: Spike parameter negative."+  | m >= 1 = error "bactrianMultSimple: Spike parameter 1.0 or larger."+  | s <= 0 = error "bactrianMultSimple: Standard deviation 0.0 or smaller."+  | otherwise = ProposalSimple $ bactrianMult m (t * s)++-- | Multiplicative proposal with kernel similar to the silhouette of a Bactrian+-- camel. See 'slideBactrian'.+scaleBactrian ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Spike parameter.+  Double ->+  -- | Standard deviation.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+scaleBactrian n w m s = createProposal n w (bactrianMultSimple m s)

@@ -0,0 +1,90 @@+-- Technically, only a Getter is needed when calculating the kernel of the proposal+-- ('kernelCont', and similar functions). I tried splitting the lens into a getter+-- and a setter. However, speed improvements were marginal, and some times not+-- even measurable. Using a 'Lens'' is just easier, and has no real drawbacks.++-- |+-- Module      :  Mcmc.Proposal.Generic+-- Description :  Generic interface to create proposals+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu May 14 20:26:27 2020.+module Mcmc.Proposal.Generic+  ( proposalGenericContinuous,+    proposalGenericDiscrete,+  )+where++import Mcmc.Proposal+import Numeric.Log+import Statistics.Distribution+import System.Random.MWC++sampleCont ::+  (ContDistr d, ContGen d) =>+  d ->+  (a -> Double -> a) ->+  Maybe (Double -> Double) ->+  a ->+  GenIO ->+  IO (a, Log Double)+sampleCont d f mfInv x g = do+  dx <- genContVar d g+  let r = case mfInv of+        Nothing -> 1.0+        Just fInv ->+          let qXY = Exp $ logDensity d dx+              qYX = Exp $ logDensity d (fInv dx)+           in qYX / qXY+  return (x `f` dx, r)+{-# INLINEABLE sampleCont #-}++-- | Generic function to create proposals for continuous parameters ('Double').+proposalGenericContinuous ::+  (ContDistr d, ContGen d) =>+  -- | Probability distribution+  d ->+  -- | Forward operator, e.g. (+), so that x + dx = x'.+  (a -> Double -> a) ->+  -- | Inverse operator, e.g., 'negate', so that x' + (negate dx) = x. Only+  -- required for biased proposals.+  Maybe (Double -> Double) ->+  ProposalSimple a+proposalGenericContinuous d f fInv = ProposalSimple $ sampleCont d f fInv++sampleDiscrete ::+  (DiscreteDistr d, DiscreteGen d) =>+  d ->+  (a -> Int -> a) ->+  Maybe (Int -> Int) ->+  a ->+  GenIO ->+  IO (a, Log Double)+sampleDiscrete d f mfInv x g = do+  dx <- genDiscreteVar d g+  let r = case mfInv of+        Nothing -> 1.0+        Just fInv ->+          let qXY = Exp $ logProbability d dx+              qYX = Exp $ logProbability d (fInv dx)+           in qYX / qXY+  return (x `f` dx, r)+{-# INLINEABLE sampleDiscrete #-}++-- | Generic function to create proposals for discrete parameters ('Int').+proposalGenericDiscrete ::+  (DiscreteDistr d, DiscreteGen d) =>+  -- | Probability distribution.+  d ->+  -- | Forward operator, e.g. (+), so that x + dx = x'.+  (a -> Int -> a) ->+  -- | Inverse operator, e.g., 'negate', so that x' + (negate dx) = x. Only+  -- required for biased proposals.+  Maybe (Int -> Int) ->+  ProposalSimple a+proposalGenericDiscrete fd f fInv = ProposalSimple $ sampleDiscrete fd f fInv

@@ -0,0 +1,57 @@+{-# LANGUAGE RankNTypes #-}++-- |+-- Module      :  Mcmc.Proposal.Scale+-- Description :  Scaling proposal with Gamma distribution+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu May 14 21:49:23 2020.+module Mcmc.Proposal.Scale+  ( scale,+    scaleUnbiased,+  )+where++import Mcmc.Proposal+import Mcmc.Proposal.Generic+import Statistics.Distribution.Gamma++-- The actual proposal with tuning parameter. The tuning parameter does not+-- change the mean.+scaleSimple :: Double -> Double -> Double -> ProposalSimple Double+scaleSimple k th t = proposalGenericContinuous (gammaDistr (k / t) (th * t)) (*) (Just recip)++-- | Multiplicative proposal with Gamma distributed kernel.+scale ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Shape.+  Double ->+  -- | Scale.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+scale n w k th = createProposal n w (scaleSimple k th)++-- | Multiplicative proposal with Gamma distributed kernel. The scale of the Gamma+-- distributions is set to (shape)^{-1}, so that the mean of the Gamma+-- distribution is 1.0.+scaleUnbiased ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Shape.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+scaleUnbiased n w k = createProposal n w (scaleSimple k (1 / k))

@@ -0,0 +1,82 @@+{-# LANGUAGE RankNTypes #-}++-- |+-- Module      :  Mcmc.Proposal.Slide+-- Description :  Normally distributed proposal+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May  6 10:59:13 2020.+module Mcmc.Proposal.Slide+  ( slide,+    slideSymmetric,+    slideUniform,+  )+where++import Mcmc.Proposal+import Mcmc.Proposal.Generic+import Statistics.Distribution.Normal+import Statistics.Distribution.Uniform++-- The actual proposal with tuning parameter.+slideSimple :: Double -> Double -> Double -> ProposalSimple Double+slideSimple m s t = proposalGenericContinuous (normalDistr m (s * t)) (+) (Just negate)++-- | Additive proposal with normally distributed kernel.+slide ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Mean.+  Double ->+  -- | Standard deviation.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+slide n w m s = createProposal n w (slideSimple m s)++-- The actual proposal with tuning parameter.+slideSymmetricSimple :: Double -> Double -> ProposalSimple Double+slideSymmetricSimple s t = proposalGenericContinuous (normalDistr 0.0 (s * t)) (+) Nothing++-- | Additive proposal with normally distributed kernel with mean zero. This+-- proposal is very fast, because the Metropolis-Hastings ratio does not include+-- calculation of the forwards and backwards kernels.+slideSymmetric ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Standard deviation.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+slideSymmetric n w s = createProposal n w (slideSymmetricSimple s)++-- The actual proposal with tuning parameter.+slideUniformSimple :: Double -> Double -> ProposalSimple Double+slideUniformSimple d t =+  proposalGenericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing++-- | Additive proposal with uniformly distributed kernel. This proposal is very fast,+-- because the Metropolis-Hastings ratio does not include calculation of the+-- forwards and backwards kernels.+slideUniform ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Delta.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal Double+slideUniform n w d = createProposal n w (slideUniformSimple d)

@@ -14,7 +14,7 @@ -- Creation date: Tue Jun 16 10:18:54 2020. -- -- Save and load an MCMC run. It is easy to save and restore the current state and--- likelihood (or the trace), but it is not feasible to store all the moves and so+-- 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.Save   ( saveStatus,@@ -28,18 +28,17 @@ import Data.Aeson.TH import qualified Data.ByteString.Lazy as B import Data.List hiding (cycle)-import Data.Map (Map) import qualified Data.Map as M-import Data.Time.Clock import Data.Vector.Unboxed (Vector) import Data.Word--- TODO: Remove as soon as split mix is used and is available with the+-- TODO: Splitmix. Reproposal as soon as split mix is used and is available with the -- statistics package. import Mcmc.Item import Mcmc.Monitor-import Mcmc.Move+import Mcmc.Proposal import Mcmc.Status hiding (save) import Mcmc.Trace+import Mcmc.Verbosity import Numeric.Log import System.IO.Unsafe (unsafePerformIO) import System.Random.MWC@@ -47,27 +46,27 @@  data Save a   = Save-      String+      -- Variables related to the chain.+      String -- Name.       (Item a)-      Int+      Int -- Iteration.       (Trace a)       (Acceptance Int)-      (Maybe Int)-      (Maybe Int)-      Int-      (Maybe (Int, UTCTime))-      Bool-      (Vector Word32)+      (Maybe Int) -- Burn in.+      (Maybe Int) -- Auto tune.+      Int -- Iterations.+      Bool -- Force.+      Bool -- Save.+      Verbosity+      (Vector Word32) -- Current seed. -$(deriveJSON defaultOptions 'Save)+      -- Variables related to the algorithm.+      [Maybe Double] -- Tuning parameters. -mapKeys :: (Ord k1, Ord k2) => [(k1, k2)] -> Map k1 v -> Map k2 v-mapKeys xs m = foldl' insrt M.empty xs-  where-    insrt m' (k1, k2) = M.insert k2 (m M.! k1) m'+$(deriveJSON defaultOptions ''Save)  toSave :: Status a -> Save a-toSave (Status nm it i tr ac br at is st sv g _ _ c _) =+toSave (Status nm it i tr ac br at is f sv vb g _ _ _ _ c _) =   Save     nm     it@@ -77,15 +76,17 @@     br     at     is-    st+    f     sv+    vb     g'+    ts   where-    moveToInt = zip (fromCycle c) [0 ..]-    ac' = Acceptance $ mapKeys moveToInt (fromAcceptance ac)-    -- TODO: Remove as soon as split mix is used and is available with the-    -- statistics package.+    ac' = transformKeysA (ccProposals c) [0 ..] ac+    -- TODO: Splitmix. Remove as soon as split mix is used and is available with+    -- the statistics package.     g' = fromSeed $ unsafePerformIO $ save g+    ts = [fmap tParam mt | mt <- map pTuner $ ccProposals c]  -- | Save a 'Status' to file. --@@ -113,7 +114,7 @@   Monitor a ->   Save a ->   Status a-fromSave p l c m (Save nm it i tr ac' br at is st sv g') =+fromSave p l c m (Save nm it i tr ac' br at is f sv vb g' ts) =   Status     nm     it@@ -123,19 +124,22 @@     br     at     is-    st+    f     sv+    vb     g+    Nothing+    Nothing     p     l-    c+    c'     m   where-    intToMove = zip [0 ..] $ fromCycle c-    ac = Acceptance $ mapKeys intToMove (fromAcceptance ac')-    -- TODO: Remove as soon as split mix is used and is available with the-    -- statistics package.+    ac = transformKeysA [0 ..] (ccProposals c) ac'+    -- TODO: Splitmix. Remove as soon as split mix is used and is available with+    -- the statistics package.     g = unsafePerformIO $ restore $ toSeed g'+    c' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals c) ts) c  -- | Load a 'Status' from file. -- Important information that cannot be saved and has to be provided again when

@@ -1,9 +1,9 @@--- TODO: Add possibility to store supplementary information about the chain.+-- XXX: Add possibility to store supplementary information about the chain. -- -- Maybe something like Trace b; and give a function a -> b to extract -- supplementary info. --- TODO: Status tuned exclusively to the Metropolis-Hastings algorithm. We+-- XXX: Status tuned exclusively to the Metropolis-Hastings algorithm. We -- should abstract the algorithm from the chain. For example, -- -- @@@ -29,68 +29,85 @@   ( Status (..),     status,     noSave,+    force,+    quiet,+    debug,   ) where +import Data.Maybe import Data.Time.Clock import Mcmc.Item import Mcmc.Monitor-import Mcmc.Move+import Mcmc.Proposal import Mcmc.Trace+import Mcmc.Verbosity (Verbosity (..)) import Numeric.Log+import System.IO import System.Random.MWC hiding (save) import Prelude hiding (cycle)  -- | The 'Status' contains all information to run an MCMC chain. It is -- constructed using the function 'status'.-data Status a = Status-  { -- Variables saved to disc.+data Status a+  = Status+      { -- MCMC related variables; saved. -    -- | The name of the MCMC chain; used as file prefix.-    name :: String,-    -- | The current 'Item' of the chain combines the current state and the-    -- current likelihood.-    item :: Item a,-    -- | The iteration is the number of completed cycles.-    iteration :: Int,-    -- | The 'Trace' of the Markov chain in reverse order, the most recent-    -- 'Item' is at the head of the list.-    trace :: Trace a,-    -- | For each 'Move', store the list of accepted (True) and rejected (False)-    -- proposals; for reasons of efficiency, the list is also stored in reverse-    -- order.-    acceptance :: Acceptance (Move a),-    -- | Number of burn in iterations; deactivate burn in with 'Nothing'.-    burnInIterations :: Maybe Int,-    -- | Auto tuning period (only during burn in); deactivate auto tuning with-    -- 'Nothing'.-    autoTuningPeriod :: Maybe Int,-    -- | Number of normal iterations excluding burn in. Note that auto tuning-    -- only happens during burn in.-    iterations :: Int,-    -- | Starting time and starting iteration of chain; used to calculate-    -- run time and ETA.-    start :: Maybe (Int, UTCTime),-    -- | Save the chain? Defaults to 'True'.-    save :: Bool,-    -- | The random number generator.-    generator :: GenIO,-    -- Auxiliary functions.+        -- | The name of the MCMC chain; used as file prefix.+        name :: String,+        -- | The current 'Item' of the chain combines the current state and the+        -- current likelihood.+        item :: Item a,+        -- | The iteration is the number of completed cycles.+        iteration :: Int,+        -- | The 'Trace' of the Markov chain in reverse order, the most recent+        -- 'Item' is at the head of the list.+        trace :: Trace a,+        -- | For each 'Proposal', store the list of accepted (True) and rejected (False)+        -- proposals; for reasons of efficiency, the list is also stored in reverse+        -- order.+        acceptance :: Acceptance (Proposal a),+        -- | Number of burn in iterations; deactivate burn in with 'Nothing'.+        burnInIterations :: Maybe Int,+        -- | Auto tuning period (only during burn in); deactivate auto tuning with+        -- 'Nothing'.+        autoTuningPeriod :: Maybe Int,+        -- | Number of normal iterations excluding burn in. Note that auto tuning+        -- only happens during burn in.+        iterations :: Int,+        -- Auxiliary variables; saved. -    -- | The prior function. The un-normalized posterior is the product of the-    -- prior and the likelihood.-    priorF :: a -> Log Double,-    -- | The likelihood function. The un-normalized posterior is the product of-    -- the prior and the likelihood.-    likelihoodF :: a -> Log Double,-    -- Variables related to the algorithm.+        -- | Overwrite output files? Default is 'False', change with 'force'.+        forceOverwrite :: Bool,+        -- | Save the chain at the end of the run? Default is 'True', change with 'noSave'.+        save :: Bool,+        -- | Verbosity.+        verbosity :: Verbosity,+        -- | The random number generator.+        generator :: GenIO,+        -- Auxiliary variables; not saved. -    -- | A set of 'Move's form a 'Cycle'.-    cycle :: Cycle a,-    -- | A 'Monitor' observing the chain.-    monitor :: Monitor a-  }+        -- | Starting time and starting iteration of chain; used to calculate+        -- run time and ETA.+        start :: Maybe (Int, UTCTime),+        -- | Handle to log file.+        logHandle :: Maybe Handle,+        -- Auxiliary functions; not saved. +        -- | The prior function. The un-normalized posterior is the product of the+        -- prior and the likelihood.+        priorF :: a -> Log Double,+        -- | The likelihood function. The un-normalized posterior is the product of+        -- the prior and the likelihood.+        likelihoodF :: a -> Log Double,+        -- Variables related to the algorithm; not saved.++        -- | A set of 'Proposal's form a 'Cycle'.+        cycle :: Cycle a,+        -- | A 'Monitor' observing the chain.+        monitor :: Monitor a+      }+ -- | Initialize the 'Status' of a Markov chain Monte Carlo run. status ::   -- | Name of the Markov chain; used as file prefix.@@ -99,7 +116,7 @@   (a -> Log Double) ->   -- | The likelihood function.   (a -> Log Double) ->-  -- | A list of 'Move's executed in forward order. The+  -- | A list of 'Proposal's executed in forward order. The   -- chain will be logged after each cycle.   Cycle a ->   -- | A 'Monitor' observing the chain.@@ -118,26 +135,45 @@   -- sure to use a generator with the same seed.   GenIO ->   Status a-status n p l c m x mB mT nI g =-  Status-    n-    i-    0-    (singletonT i)-    (emptyA $ fromCycle c)-    mB-    mT-    nI-    Nothing-    True-    g-    p-    l-    c-    m+status n p l c m x mB mT nI g+  | isJust mT && isNothing mB = error "status: Auto tuning period given, but no burn in."+  | otherwise =+    Status+      n+      i+      0+      (singletonT i)+      (emptyA $ ccProposals c)+      mB+      mT+      nI+      False+      True+      Info+      g+      Nothing+      Nothing+      p+      l+      c+      m   where     i = Item x (p x) (l x)  -- | Do not save the Markov chain at the end. noSave :: Status a -> Status a noSave s = s {save = False}++-- | Overwrite existing files; it is not necessary to use 'force', when a chain+-- is continued.+force :: Status a -> Status a+force s = s {forceOverwrite = True}++-- | Do not print anything to standard output. Do not create log file. File+-- monitors and batch monitors are executed normally.+quiet :: Status a -> Status a+quiet s = s {verbosity = Quiet}++-- | Be verbose.+debug :: Status a -> Status a+debug s = s {verbosity = Debug}

@@ -36,6 +36,15 @@ shuffleN :: [a] -> Int -> GenIO -> IO [[a]] shuffleN xs n = grabble xs n (length xs) +-- -- Using System.Random.Shuffle. Speed is the same, so stay without additional dependency.+-- -- | Shuffle a list @n@ times.+-- shuffleN :: [a] -> Int -> GenIO -> IO [[a]]+-- shuffleN xs n g = replicateM n $ fmap (shuffle xs) (rseqM (length xs - 1) g)+--   where+--     rseqM :: Int -> GenIO -> IO [Int]+--     rseqM 0 _ = return []+--     rseqM i gen = liftM2 (:) (uniformR (0, i) gen) (rseqM (i - 1) gen)+ -- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@ -- elements from @xs@, without replacement, and that @m@ times. grabble :: [a] -> Int -> Int -> GenIO -> IO [[a]]

@@ -40,8 +40,6 @@ instance FromJSON a => FromJSON (Trace a) where   parseJSON v = Trace <$> parseJSONList v --- $(deriveJSON defaultOptions 'Trace)- -- | The empty trace. singletonT :: Item a -> Trace a singletonT i = Trace [i]

@@ -0,0 +1,40 @@+{-# LANGUAGE TemplateHaskell #-}++-- |+-- Module      :  Mcmc.Verbosity+-- Description :  Be quiet! Or better not?+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Sat Jun 27 10:49:28 2020.+module Mcmc.Verbosity+  ( Verbosity (..),+    warn,+    info,+    debug,+  )+where++import Control.Monad+import Data.Aeson.TH++-- | Not much to say here.+data Verbosity = Quiet | Warn | Info | Debug deriving (Show, Eq, Ord)++$(deriveJSON defaultOptions ''Verbosity)++-- | Perform action if 'Verbosity' is 'Warn' or higher.+warn :: Applicative m => Verbosity -> m () -> m ()+warn v = when (v >= Warn)++-- | Perform action if 'Verbosity' is 'Info' or higher.+info :: Applicative m => Verbosity -> m () -> m ()+info v = when (v >= Info)++-- | Perform action if 'Verbosity' is 'Debug'.+debug :: Applicative m => Verbosity -> m () -> m ()+debug v = when (v == Debug)

@@ -1,34 +0,0 @@--- |--- Module      :  Mcmc.Move.SlideSpec--- Description :  Unit tests for Mcmc.Move.SlideSpec--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue May 19 12:07:43 2020.-module Mcmc.Move.SlideSpec-  ( spec,-  )-where--import Data.Maybe-import Mcmc.Move-import Mcmc.Move.Slide-import Test.Hspec-import Test.QuickCheck--prop_sym :: Eq b => (a -> a -> b) -> a -> a -> Bool-prop_sym f x y = f x y == f y x--slideSym :: Move Double-slideSym = slide "symmetric" 1 id 0 1.0 False--spec :: Spec-spec =-  describe "slide"-    $ it "has a symmetric proposal distribution if mean is 0"-    $ property-    $ prop_sym (fromJust $ mvDensity . mvSimple $ slideSym)

@@ -0,0 +1,43 @@+-- |+--   Module      :  Mcmc.ProposalSpec+--   Description :  Unit tests for Mcmc.Proposal+--   Copyright   :  (c) Dominik Schrempf, 2020+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Thu Jun 25 11:46:05 2020.+module Mcmc.ProposalSpec+  ( spec,+  )+where++import Mcmc.Proposal+import Mcmc.Proposal.Slide+import System.Random.MWC+import Test.Hspec++p1 :: Proposal Double+p1 = slideSymmetric "test1" 1 1.0 True++p2 :: Proposal Double+p2 = slideSymmetric "test2" 3 1.0 True++c :: Cycle Double+c = fromList [p1, p2]++spec :: Spec+spec =+  describe "getNCycles"+    $ it "returns the correct number of proposals in a cycle"+    $ do+      g <- create+      l1 <- length . head <$> getNCycles c 1 g+      l1 `shouldBe` 4+      l2 <- length . head <$> getNCycles (setOrder RandomReversibleO c) 1 g+      l2 `shouldBe` 8+      o3 <- head <$> getNCycles (setOrder SequentialReversibleO c) 1 g+      length o3 `shouldBe` 8+      o3 `shouldBe` [p1, p2, p2, p2, p2, p2, p2, p1]

@@ -14,13 +14,16 @@   ) where -import Mcmc hiding (save)+import Mcmc+import Mcmc.Save+import Mcmc.Status hiding (save) import Numeric.Log import Statistics.Distribution hiding   ( mean,     stdDev,   ) import Statistics.Distribution.Normal+import System.Directory import System.Random.MWC import Test.Hspec @@ -33,17 +36,17 @@ lh :: Double -> Log Double lh = Exp . logDensity (normalDistr trueMean trueStdDev) -moveCycle :: Cycle Double-moveCycle =+proposals :: Cycle Double+proposals =   fromList-    [ slideSymmetric "small" 5 id 0.1 True,-      slideSymmetric "medium" 2 id 1.0 True,-      slideSymmetric "large" 2 id 5.0 True,-      slide "skewed" 1 id 1.0 4.0 True+    [ slideSymmetric "small" 5 0.1 True,+      slideSymmetric "medium" 2 1.0 True,+      slideSymmetric "large" 2 5.0 True,+      slide "skewed" 1 1.0 4.0 True     ]  monStd :: MonitorStdOut Double-monStd = monitorStdOut [monitorRealFloat "mu" id] 10+monStd = monitorStdOut [monitorRealFloat "mu"] 10  mon :: Monitor Double mon = Monitor monStd [] []@@ -63,33 +66,35 @@     $ it "doesn't change the MCMC chain"     $ do       gen <- create-      let s = noSave $ status "SaveSpec" (const 1) lh moveCycle mon 0 nBurn nAutoTune nIter gen+      let s =+            force $ quiet $ noSave $+              status "SaveSpec" (const 1) lh proposals mon 0 nBurn nAutoTune nIter gen       saveStatus "SaveSpec.json" s-      s' <- loadStatus (const 1) lh moveCycle mon "SaveSpec.json"+      s' <- loadStatus (const 1) lh proposals mon "SaveSpec.json"       r <- mh s       r' <- mh s'+      removeFile "SaveSpec.json"       item r `shouldBe` item r'       iteration r `shouldBe` iteration r'       trace r `shouldBe` trace r'       g <- save $ generator r       g' <- save $ generator r'       g `shouldBe` g'--  -- -- TODO: This will only work with a splittable generator because getNCycles-  -- -- changes the generator.-  -- describe "mhContinue"-  --   $ it "mh 200 + mhContinue 200 == mh 400"-  --   $ do-  --     gen1 <- create-  --     let s1 = noSave $ status "SaveSpec" (const 1) likelihood moveCycle mon 0 nBurn nAutoTune 400 gen1-  --     r1 <- mh s1-  --     gen2 <- create-  --     let s2 = noSave $ status "SaveSpec" (const 1) likelihood moveCycle mon 0 nBurn nAutoTune 200 gen2-  --     r2' <- mh s2-  --     r2 <- mhContinue 200 r2'-  --     item r1 `shouldBe` item r2-  --     iteration r1 `shouldBe` iteration r2-  --     trace r1 `shouldBe` trace r2-  --     g <- save $ generator r1-  --     g' <- save $ generator r2-  --     g `shouldBe` g'+-- -- TODO: Splitmix. This will only work with a splittable generator+-- -- because getNCycles changes the generator.+-- describe "mhContinue"+--   $ it "mh 200 + mhContinue 200 == mh 400"+--   $ do+--     gen1 <- create+--     let s1 = noSave $ status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 400 gen1+--     r1 <- mh s1+--     gen2 <- create+--     let s2 = noSave $ status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 200 gen2+--     r2' <- mh s2+--     r2 <- mhContinue 200 r2'+--     item r1 `shouldBe` item r2+--     iteration r1 `shouldBe` iteration r2+--     trace r1 `shouldBe` trace r2+--     g <- save $ generator r1+--     g' <- save $ generator r2+--     g `shouldBe` g'