diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@

-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 diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -23,5 +23,6 @@ g <- create defaultMain [ bench "Normal" $ nfIO (normalBench g), + bench "NormalBactrian" $ nfIO (normalBactrianBench g), bench "Poisson" $ nfIO (poissonBench g) ] diff --git a/bench/Normal.hs b/bench/Normal.hs --- a/bench/Normal.hs +++ b/bench/Normal.hs @@ -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 diff --git a/bench/Poisson.hs b/bench/Poisson.hs --- a/bench/Poisson.hs +++ b/bench/Poisson.hs @@ -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 diff --git a/mcmc.cabal b/mcmc.cabal --- a/mcmc.cabal +++ b/mcmc.cabal @@ -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 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 diff --git a/src/Mcmc.hs b/src/Mcmc.hs --- a/src/Mcmc.hs +++ b/src/Mcmc.hs @@ -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) diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs --- a/src/Mcmc/Mcmc.hs +++ b/src/Mcmc/Mcmc.hs @@ -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 diff --git a/src/Mcmc/Metropolis.hs b/src/Mcmc/Metropolis.hs --- a/src/Mcmc/Metropolis.hs +++ b/src/Mcmc/Metropolis.hs @@ -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." diff --git a/src/Mcmc/Monitor.hs b/src/Mcmc/Monitor.hs --- a/src/Mcmc/Monitor.hs +++ b/src/Mcmc/Monitor.hs @@ -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) diff --git a/src/Mcmc/Monitor/Parameter.hs b/src/Mcmc/Monitor/Parameter.hs --- a/src/Mcmc/Monitor/Parameter.hs +++ b/src/Mcmc/Monitor/Parameter.hs @@ -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 #-} diff --git a/src/Mcmc/Monitor/ParameterBatch.hs b/src/Mcmc/Monitor/ParameterBatch.hs --- a/src/Mcmc/Monitor/ParameterBatch.hs +++ b/src/Mcmc/Monitor/ParameterBatch.hs @@ -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) diff --git a/src/Mcmc/Monitor/Time.hs b/src/Mcmc/Monitor/Time.hs --- a/src/Mcmc/Monitor/Time.hs +++ b/src/Mcmc/Monitor/Time.hs @@ -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 diff --git a/src/Mcmc/Move.hs b/src/Mcmc/Move.hs deleted file mode 100644 --- a/src/Mcmc/Move.hs +++ /dev/null @@ -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 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 diff --git a/src/Mcmc/Move/Generic.hs b/src/Mcmc/Move/Generic.hs deleted file mode 100644 --- a/src/Mcmc/Move/Generic.hs +++ /dev/null @@ -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 diff --git a/src/Mcmc/Move/Scale.hs b/src/Mcmc/Move/Scale.hs deleted file mode 100644 --- a/src/Mcmc/Move/Scale.hs +++ /dev/null @@ -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 diff --git a/src/Mcmc/Move/Slide.hs b/src/Mcmc/Move/Slide.hs deleted file mode 100644 --- a/src/Mcmc/Move/Slide.hs +++ /dev/null @@ -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 diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Proposal.hs @@ -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 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 diff --git a/src/Mcmc/Proposal/Bactrian.hs b/src/Mcmc/Proposal/Bactrian.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Proposal/Bactrian.hs @@ -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) diff --git a/src/Mcmc/Proposal/Generic.hs b/src/Mcmc/Proposal/Generic.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Proposal/Generic.hs @@ -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 diff --git a/src/Mcmc/Proposal/Scale.hs b/src/Mcmc/Proposal/Scale.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Proposal/Scale.hs @@ -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)) diff --git a/src/Mcmc/Proposal/Slide.hs b/src/Mcmc/Proposal/Slide.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Proposal/Slide.hs @@ -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) diff --git a/src/Mcmc/Save.hs b/src/Mcmc/Save.hs --- a/src/Mcmc/Save.hs +++ b/src/Mcmc/Save.hs @@ -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 diff --git a/src/Mcmc/Status.hs b/src/Mcmc/Status.hs --- a/src/Mcmc/Status.hs +++ b/src/Mcmc/Status.hs @@ -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} diff --git a/src/Mcmc/Tools/Shuffle.hs b/src/Mcmc/Tools/Shuffle.hs --- a/src/Mcmc/Tools/Shuffle.hs +++ b/src/Mcmc/Tools/Shuffle.hs @@ -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]] diff --git a/src/Mcmc/Trace.hs b/src/Mcmc/Trace.hs --- a/src/Mcmc/Trace.hs +++ b/src/Mcmc/Trace.hs @@ -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] diff --git a/src/Mcmc/Verbosity.hs b/src/Mcmc/Verbosity.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Verbosity.hs @@ -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) diff --git a/test/Mcmc/Move/SlideSpec.hs b/test/Mcmc/Move/SlideSpec.hs deleted file mode 100644 --- a/test/Mcmc/Move/SlideSpec.hs +++ /dev/null @@ -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) diff --git a/test/Mcmc/ProposalSpec.hs b/test/Mcmc/ProposalSpec.hs new file mode 100644 --- /dev/null +++ b/test/Mcmc/ProposalSpec.hs @@ -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] diff --git a/test/Mcmc/SaveSpec.hs b/test/Mcmc/SaveSpec.hs --- a/test/Mcmc/SaveSpec.hs +++ b/test/Mcmc/SaveSpec.hs @@ -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'