diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -5,6 +5,16 @@ ## Unreleased changes +## 0.5.0.0 + +- Marginal likelihood estimation using thermodynamic integration or stepping + stone sampling. +- Various changes of function names (e.g., metropologis -> mhg). +- Updated examples. +- Proper but minimal logging framework. +- Various other changes. + + ## 0.4.0.0 - Greatly improve documentation. diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -1,6 +1,6 @@ -- | -- Module : Main --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com diff --git a/bench/Normal.hs b/bench/Normal.hs --- a/bench/Normal.hs +++ b/bench/Normal.hs @@ -1,7 +1,7 @@ -- | -- Module : Normal -- Description : Benchmark Metropolis-Hastings-Green algorithm --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -31,7 +31,7 @@ lh = normal trueMean stdDev cc :: Cycle Double -cc = cycleFromList [slideSymmetric 1.0 (PName "Medium") (PWeight 1) Tune] +cc = cycleFromList [slideSymmetric 1.0 (PName "Medium") (pWeight 1) Tune] mons :: [MonitorParameter Double] mons = [monitorDouble "mu"] @@ -52,14 +52,15 @@ Overwrite Sequential NoSave + LogStdOutOnly Quiet - a <- mhg noPrior lh cc mon 0 g + a <- mhg noPrior lh cc mon TraceAuto 0 g void $ mcmc s a ccLarge :: Cycle Double ccLarge = cycleFromList - [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (PWeight 1) Tune | i <- [0 .. 100 :: Int]] + [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (pWeight 1) Tune | i <- [0 .. 100 :: Int]] -- Should have the same run time as 'normalSlide'. normalLargeCycleBench :: GenIO -> IO () @@ -72,12 +73,13 @@ Overwrite Sequential NoSave + LogStdOutOnly Quiet - a <- mhg noPrior lh ccLarge mon 0 g + a <- mhg noPrior lh ccLarge mon TraceAuto 0 g void $ mcmc s a ccBactrian :: Cycle Double -ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (PWeight 1) Tune] +ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (pWeight 1) Tune] normalBactrianBench :: GenIO -> IO () normalBactrianBench g = do @@ -89,8 +91,9 @@ Overwrite Sequential NoSave + LogStdOutOnly Quiet - a <- mhg noPrior lh ccBactrian mon 0 g + a <- mhg noPrior lh ccBactrian mon TraceAuto 0 g void $ mcmc s a normalMC3 :: GenIO -> Int -> IO () @@ -103,7 +106,8 @@ Overwrite Sequential NoSave + LogStdOutOnly Quiet mc3S = MC3Settings (NChains n) (SwapPeriod 2) (NSwaps 1) - a <- mc3 mc3S noPrior lh cc mon 0 g + a <- mc3 mc3S noPrior lh cc mon TraceAuto 0 g void $ mcmc mcmcS a diff --git a/bench/Poisson.hs b/bench/Poisson.hs --- a/bench/Poisson.hs +++ b/bench/Poisson.hs @@ -1,7 +1,7 @@ -- | -- Module : Poisson -- Description : Poisson regression model for airline fatalities --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -19,7 +19,6 @@ import Control.Monad import Lens.Micro import Mcmc -import Numeric.Log hiding (sum) import System.Random.MWC type I = (Double, Double) @@ -42,10 +41,10 @@ lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears] proposalAlpha :: Proposal I -proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (PWeight 1) NoTune +proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune proposalBeta :: Proposal I -proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (PWeight 1) NoTune +proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune proposals :: Cycle I proposals = cycleFromList [proposalAlpha, proposalBeta] @@ -75,6 +74,7 @@ Overwrite Sequential NoSave + LogStdOutOnly Quiet - a <- mhg noPrior lh proposals mon initial g + a <- mhg noPrior lh proposals mon TraceAuto initial g void $ mcmc s a diff --git a/mcmc.cabal b/mcmc.cabal --- a/mcmc.cabal +++ b/mcmc.cabal @@ -1,8 +1,8 @@ cabal-version: 2.2 name: mcmc -version: 0.4.0.0 +version: 0.5.0.0 license: GPL-3.0-or-later -copyright: Dominik Schrempf (2020) +copyright: Dominik Schrempf (2021) maintainer: dominik.schrempf@gmail.com author: Dominik Schrempf homepage: https://github.com/dschrempf/mcmc#readme @@ -25,19 +25,23 @@ exposed-modules: Mcmc Mcmc.Algorithm - Mcmc.Algorithm.Metropolis Mcmc.Algorithm.MC3 + Mcmc.Algorithm.MHG Mcmc.Chain.Chain Mcmc.Chain.Link Mcmc.Chain.Save Mcmc.Chain.Trace Mcmc.Environment + Mcmc.Likelihood + Mcmc.Logger + Mcmc.MarginalLikelihood Mcmc.Mcmc Mcmc.Monitor Mcmc.Monitor.Log Mcmc.Monitor.Parameter Mcmc.Monitor.ParameterBatch Mcmc.Monitor.Time + Mcmc.Posterior Mcmc.Prior Mcmc.Proposal Mcmc.Proposal.Bactrian @@ -46,6 +50,7 @@ Mcmc.Proposal.Slide Mcmc.Proposal.Simplex Mcmc.Settings + Mcmc.Statistics.Types hs-source-dirs: src other-modules: @@ -58,27 +63,27 @@ default-language: Haskell2010 ghc-options: -Wall -Wunused-packages build-depends: - aeson >=1.5.4.1 && <1.6, + aeson >=1.5.6.0, base >=4.7 && <5, - bytestring >=0.10.10.0 && <0.11, - circular >=0.3.1.1 && <0.4, - containers >=0.6.2.1 && <0.7, - data-default >=0.7.1.1 && <0.8, - deepseq >=1.4.4.0 && <1.5, - directory >=1.3.6.0 && <1.4, - dirichlet >=0.1.0.0 && <0.2, - double-conversion >=2.0.2.0 && <2.1, - log-domain ==0.13.*, - microlens >=0.4.11.2 && <0.5, - mwc-random >=0.14.0.0 && <0.15, - monad-parallel >=0.7.2.3 && <0.8, - pretty-show ==1.10.*, - primitive >=0.7.1.0 && <0.8, - statistics >=0.15.2.0 && <0.16, - time >=1.9.3 && <1.10, - transformers >=0.5.6.2 && <0.6, - vector >=0.12.1.2 && <0.13, - zlib >=0.6.2.2 && <0.7 + bytestring >=0.10.12.0, + circular >=0.4.0.0, + containers >=0.6.2.1, + data-default >=0.7.1.1, + deepseq >=1.4.4.0, + directory >=1.3.6.0, + dirichlet >=0.1.0.4, + double-conversion >=2.0.2.0, + log-domain >=0.13.1, + microlens >=0.4.12.0, + mwc-random >=0.15.0.1, + monad-parallel >=0.7.2.4, + pretty-show >=1.10, + primitive >=0.7.1.0, + statistics >=0.15.2.0, + time >=1.9.3, + transformers >=0.5.6.2, + vector >=0.12.3.0, + zlib >=0.6.2.3 test-suite mcmc-test type: exitcode-stdio-1.0 @@ -93,11 +98,11 @@ ghc-options: -Wall -Wunused-packages build-depends: base >=4.7 && <5, - hspec >=2.7.4 && <2.8, - log-domain ==0.13.*, + hspec >=2.7.10, + log-domain >=0.13.1, mcmc -any, - mwc-random >=0.14.0.0 && <0.15, - statistics >=0.15.2.0 && <0.16 + mwc-random >=0.15.0.1, + statistics >=0.15.2.0 benchmark mcmc-bench type: exitcode-stdio-1.0 @@ -112,8 +117,8 @@ ghc-options: -Wall -Wunused-packages build-depends: base >=4.7 && <5, - criterion >=1.5.7.0 && <1.6, - log-domain ==0.13.*, + criterion >=1.5.9.0, + log-domain >=0.13.1, mcmc -any, - microlens >=0.4.11.2 && <0.5, - mwc-random >=0.14.0.0 && <0.15 + microlens >=0.4.12.0, + mwc-random >=0.15.0.1 diff --git a/src/Mcmc.hs b/src/Mcmc.hs --- a/src/Mcmc.hs +++ b/src/Mcmc.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc -- Description : Markov chain Monte Carlo samplers, batteries included --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -76,9 +76,13 @@ -- Gibbs sampler, or a random sequence scan Gibbs sampler, respectively. See -- 'Order'. -- - -- Note that it is of utter importance that the given 'Cycle' enables - -- traversal of the complete state space. Otherwise, the Markov chain will - -- not converge to the correct stationary posterior distribution. + -- Notes: + -- - It is important that the given 'Cycle' enables traversal of the + -- complete state space. Otherwise, the Markov chain will not converge to + -- the correct stationary posterior distribution. + -- - Be careful when assigning proposals because acceptance ratios may have + -- to be amended using Jacobians. Please see an [example involving a pair + -- of numbers](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs). -- -- Proposals are named according to what they do, i.e., how they change the -- state of a Markov chain, and not according to the intrinsically used @@ -114,8 +118,12 @@ -- http://dx.doi.org/10.1093/sysbio/syw021 PName (..), PWeight (..), + pWeight, Proposal, + JacobianFunction, (@~), + liftProposal, + liftProposalWith, Tune (..), scale, scaleUnbiased, @@ -157,36 +165,39 @@ monitorFile, MonitorBatch, monitorBatch, + simpleMonitor, module Mcmc.Monitor.Parameter, module Mcmc.Monitor.ParameterBatch, - -- * Prior distributions - - -- | Convenience functions for computing priors. + -- * Prior, likelihood, and posterior values and functions module Mcmc.Prior, + module Mcmc.Likelihood, -- * MCMC samplers mcmc, mcmcContinue, - -- | See also 'settingsLoad', 'mhgLoad', and 'mc3Load'. -- * Algorithms - module Mcmc.Algorithm.Metropolis, + module Mcmc.Algorithm.MHG, module Mcmc.Algorithm.MC3, + -- * Marginal likelihood calculation + module Mcmc.MarginalLikelihood, - -- * Useful type synonyms - PriorFunction, - noPrior, - LikelihoodFunction, - noLikelihood, + -- * Useful types + module Mcmc.Statistics.Types, + + -- * Useful re-exports + module Numeric.Log, ) where import Mcmc.Algorithm.MC3 -import Mcmc.Algorithm.Metropolis +import Mcmc.Algorithm.MHG import Mcmc.Chain.Chain +import Mcmc.Likelihood +import Mcmc.MarginalLikelihood import Mcmc.Mcmc import Mcmc.Monitor import Mcmc.Monitor.Parameter @@ -198,3 +209,5 @@ import Mcmc.Proposal.Simplex import Mcmc.Proposal.Slide import Mcmc.Settings +import Mcmc.Statistics.Types +import Numeric.Log hiding (sum) diff --git a/src/Mcmc/Algorithm.hs b/src/Mcmc/Algorithm.hs --- a/src/Mcmc/Algorithm.hs +++ b/src/Mcmc/Algorithm.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Algorithm -- Description : Algortihms for Markov chain Monte Carlo samplers --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -25,6 +25,10 @@ -- | Current iteration. aIteration :: a -> Int + + -- | Check if the current state is invalid. At the moment this should just + -- check if the posterior probability is zero or NaN. + aIsInValidState :: a -> Bool -- | Sample the next state. aIterate :: ParallelizationMode -> a -> IO a diff --git a/src/Mcmc/Algorithm/MC3.hs b/src/Mcmc/Algorithm/MC3.hs --- a/src/Mcmc/Algorithm/MC3.hs +++ b/src/Mcmc/Algorithm/MC3.hs @@ -5,7 +5,7 @@ -- | -- Module : Mcmc.Algorithm.MC3 -- Description : Metropolis-coupled Markov chain Monte Carlo algorithm --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -19,8 +19,8 @@ -- Also known as parallel tempering. -- -- Like any other parallel MCMC algorithm, the 'MC3' algorithm is essentially an --- 'Mcmc.Algorithm.Metropolis.MHG' algorithm on the product space of all --- parallel chains. +-- 'Mcmc.Algorithm.MHG.MHG' algorithm on the product space of all parallel +-- chains. -- -- For example, see -- @@ -64,7 +64,7 @@ import Data.Word -- import Debug.Trace hiding (trace) import Mcmc.Algorithm -import Mcmc.Algorithm.Metropolis +import Mcmc.Algorithm.MHG import Mcmc.Chain.Chain import Mcmc.Chain.Link import Mcmc.Chain.Save @@ -72,11 +72,14 @@ import Mcmc.Internal.Random import Mcmc.Internal.Shuffle import Mcmc.Monitor +import Mcmc.Posterior import Mcmc.Proposal import Mcmc.Settings import Numeric.Log hiding (sum) import System.Random.MWC +-- import Debug.Trace hiding (trace) + -- | Total number of parallel chains. -- -- Must be two or larger. @@ -146,9 +149,17 @@ SavedMC3 a -> IO (MC3 a) fromSavedMC3 pr lh cc mn (SavedMC3 s scs bs i ac g') = do - mhgs <- V.mapM (fmap MHG . fromSavedChain pr lh cc mn) scs + mhgs <- + V.fromList + <$> sequence + [ MHG <$> fromSavedChain pf lf cc mn sc + | (sc, pf, lf) <- zip3 (V.toList scs) prs lhs + ] g <- loadGen g' return $ MC3 s mhgs bs i ac g + where + prs = map (heatFunction pr) $ U.toList bs + lhs = map (heatFunction lh) $ U.toList bs -- | The MC3 algorithm. data MC3 a = MC3 @@ -167,6 +178,7 @@ instance ToJSON a => Algorithm (MC3 a) where aName = const "Metropolis-coupled Markov chain Monte Carlo (MC3)" aIteration = mc3Iteration + aIsInValidState = mc3IsInValidState aIterate = mc3Iterate aAutoTune = mc3AutoTune aResetAcceptance = mc3ResetAcceptance @@ -177,6 +189,20 @@ aCloseMonitors = mc3CloseMonitors aSave = mc3Save +heatFunction :: + -- Cold Function. + (a -> Log Double) -> + -- Reciprocal temperature. + Double -> + -- The heated prior or likelihood function + (a -> Log Double) +heatFunction f b + | b <= 0 = error "heatFunction: Reciprocal temperature is zero or negative." + | b == 1.0 = f + | otherwise = (** b') . f + where + b' = Exp $ log b + -- The prior and likelihood values of the current link are updated. -- -- NOTE: The trace is not changed! In particular, the prior and likelihood @@ -191,7 +217,7 @@ Double -> MHG a -> MHG a -setReciprocalTemperature prf lhf beta a = +setReciprocalTemperature coldPrf coldLhf b a = MHG $ c { priorFunction = prf', @@ -200,15 +226,14 @@ } where c = fromMHG a - b' = Exp $ log beta -- We need twice the amount of computations compared to taking the power -- after calculating the posterior (pr x * lh x) ** b'. However, I don't -- think this is a serious problem. -- -- To minimize computations, it is key to avoid modification of the -- reciprocal temperature for the cold chain. - prf' = (** b') . prf - lhf' = (** b') . lhf + prf' = heatFunction coldPrf b + lhf' = heatFunction coldLhf b x = state $ link c initMHG :: @@ -224,14 +249,14 @@ IO (MHG a) initMHG prf lhf i beta a | i < 0 = error "initMHG: Chain index negative." - -- Do not temper with the cold chain. - | i == 0 = return a + -- Only set the id for the cold chain. + | i == 0 = return $ MHG $ c {chainId = Just 0} | otherwise = do -- We have to push the current link in the trace, since it is not set by -- 'setReciprocalTemperature'. The other links in the trace are still -- pointing to the link of the cold chain, but this has no effect. t' <- pushT l t - return $ MHG $ c {chainId = i, trace = t'} + return $ MHG $ c {chainId = Just i, trace = t'} where a' = setReciprocalTemperature prf lhf beta a c = fromMHG a' @@ -254,17 +279,18 @@ LikelihoodFunction a -> Cycle a -> Monitor a -> - a -> + TraceLength -> + InitialState a -> GenIO -> IO (MC3 a) -mc3 s pr lh cc mn i0 g +mc3 s pr lh cc mn tr i0 g | n < 2 = error "mc3: The number of chains must be two or larger." | sp < 1 = error "mc3: The swap period must be strictly positive." | sn < 1 || sn > n - 1 = error "mc3: The number of swaps must be in [1, NChains - 1]." | otherwise = do -- Split random number generators. gs <- V.fromList <$> splitGen n g - cs <- V.mapM (mhg pr lh cc mn i0) gs + cs <- V.mapM (mhg pr lh cc mn tr i0) gs hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs return $ MC3 s hcs bs 0 (emptyA [0 .. n - 2]) g where @@ -276,10 +302,10 @@ -- -- NOTE: Have to 'take n' elements, because vectors are not as lazy as -- lists. - bs = U.fromList $ take n $ iterate (* 0.92) 1.0 + bs = U.fromList $ take n $ iterate (* 0.97) 1.0 mc3Fn :: AnalysisName -> FilePath -mc3Fn (AnalysisName nm) = nm ++ ".mc3" +mc3Fn (AnalysisName nm) = nm ++ ".mcmc.mc3" -- | Save an MC3 algorithm. mc3Save :: @@ -315,7 +341,7 @@ -- Index j>=0, j/=i of right chain. Int -> MHGChains a -> - (MHGChains a, Log Double) + (MHGChains a, Posterior) swapWith i j xs | i < 0 = error "swapWith: Left index is negative." | j < 0 = error "swapWith: Right index is negative." @@ -356,12 +382,28 @@ Int -> IO (MC3 a) mc3ProposeSwap a i = do + let cs = mc3MHGChains a + -- -- Debug. + -- prL = prior $ link $ fromMHG $ cs V.! i + -- prR = prior $ link $ fromMHG $ cs V.! (i+1) + -- lhL = likelihood $ link $ fromMHG $ cs V.! i + -- lhR = likelihood $ link $ fromMHG $ cs V.! (i+1) -- 1. Sample new state and get the Metropolis ratio. - let (!y, !r) = swapWith i (i + 1) $ mc3MHGChains a + let (!y, !r) = swapWith i (i + 1) cs -- 2. Accept or reject. accept <- mhgAccept r g if accept then do + -- -- Debug. + -- traceIO $ "Swap accepted: " <> show i <> " <-> " <> show (i+1) + -- let prL' = prior $ link $ fromMHG $ y V.! i + -- prR' = prior $ link $ fromMHG $ y V.! (i+1) + -- lhL' = likelihood $ link $ fromMHG $ y V.! i + -- lhR' = likelihood $ link $ fromMHG $ y V.! (i+1) + -- traceIO $ "Log priors (left, right, before swap): " <> show (ln prL) <> " " <> show (ln prR) + -- traceIO $ "Log priors (left, right, after swap): " <> show (ln prL') <> " " <> show (ln prR') + -- traceIO $ "Log likelihoods (left, right, before swap): " <> show (ln lhL) <> " " <> show (ln lhR) + -- traceIO $ "Log likelihood (left, right, after swap): " <> show (ln lhL') <> " " <> show (ln lhR') let !ac' = pushA i True (mc3SwapAcceptance a) return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'} else do @@ -370,9 +412,15 @@ where g = mc3Generator a +mc3IsInValidState :: ToJSON a => MC3 a -> Bool +mc3IsInValidState a = V.any aIsInValidState mhgs + where + mhgs = mc3MHGChains a + -- TODO: Splimix. 'mc3Iterate' is actually not parallel, but concurrent because -- of the IO constraint. Use pure parallel code when we have a pure generator. --- However, we have honor the mutable traces. +-- +-- However, we to have honor the mutable traces. mc3Iterate :: ToJSON a => ParallelizationMode -> @@ -399,7 +447,7 @@ Parallel -> -- See 'Control.Monad.Parallel' of package 'monad-parallel'. Go via a -- list, and use 'forkIO'. - V.fromList <$> P.mapM (aIterate pm) (V.toList (mc3MHGChains a')) + V.fromList <$> P.mapM (aIterate pm) (V.toList $ mc3MHGChains a') let i = mc3Iteration a' return $ a' {mc3MHGChains = mhgs, mc3Iteration = succ i} @@ -432,7 +480,7 @@ mhgs' = V.map aAutoTune mhgs -- 2. Auto tune temperatures. optimalRate = getOptimalRate PDimensionUnknown - currentRates = acceptanceRates $ mc3SwapAcceptance a + mCurrentRates = acceptanceRates $ mc3SwapAcceptance a -- We assume that the acceptance rate of state swaps between two chains is -- roughly proportional to the ratio of the temperatures of the chains. -- Hence, we focus on temperature ratios, actually reciprocal temperature @@ -442,7 +490,9 @@ -- -- The factor (1/2) was determined by a few tests and is otherwise -- absolutely arbitrary. - xi i = exp $ (/ 2) $ (currentRates M.! i) - optimalRate + xi i = case mCurrentRates M.! i of + Nothing -> 1.0 + Just currentRate -> exp $ (/ 2) $ currentRate - optimalRate bs = mc3ReciprocalTemperatures a n = fromNChains $ mc3NChains $ mc3Settings a -- Do not change the temperature, and the prior and likelihood functions of @@ -481,9 +531,12 @@ [ "MC3: Cycle of cold chain.", coldMHGCycleSummary ] - ++ [ "MC3: Average acceptance rate across all chains: " <> BL.fromStrict (BC.toFixed 2 ar) - | not $ isNaN ar - ] + ++ case mAr of + Nothing -> [] + Just ar -> + [ "MC3: Average acceptance rate across all chains: " + <> BL.fromStrict (BC.toFixed 2 ar) + ] ++ [ "MC3: Reciprocal temperatures of the chains: " <> BL.intercalate ", " bsB <> ".", "MC3: Summary of state swaps.", "MC3: The swap period is " <> swapPeriodB <> ".", @@ -494,7 +547,7 @@ ++ [ summarizeProposal (PName $ show i ++ " <-> " ++ show (i + 1)) (PDescription "Swap states between chains") - (PWeight 1) + (pWeight 1) (Just $ bs U.! (i + 1)) PDimensionUnknown (acceptanceRate i swapAcceptance) @@ -505,9 +558,12 @@ mhgs = mc3MHGChains a coldMHGCycleSummary = aSummarizeCycle $ V.head mhgs cs = V.map fromMHG mhgs - as = V.map (acceptanceRates . acceptance) cs - vAr = V.map (\m -> sum m / fromIntegral (length m)) as - ar = V.sum vAr / fromIntegral (V.length vAr) + -- Acceptance rates may be 'Nothing' when no proposals have been undertaken. + -- The 'sequence' operations pull the 'Nothing's out of the inner + -- structures. + as = sequence $ V.map (sequence . acceptanceRates . acceptance) cs + mVecAr = V.map (\mp -> sum mp / fromIntegral (length mp)) <$> as + mAr = (\vec -> V.sum vec / fromIntegral (V.length vec)) <$> mVecAr bs = mc3ReciprocalTemperatures a bsB = map (BL.fromStrict . BC.toFixed 2) $ U.toList bs swapPeriod = fromSwapPeriod $ mc3SwapPeriod $ mc3Settings a diff --git a/src/Mcmc/Algorithm/MHG.hs b/src/Mcmc/Algorithm/MHG.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Algorithm/MHG.hs @@ -0,0 +1,283 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE OverloadedStrings #-} + +-- | +-- Module : Mcmc.Algorithm.MHG +-- Description : Metropolis-Hastings-Green algorithm +-- Copyright : (c) Dominik Schrempf 2021 +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : unstable +-- Portability : portable +-- +-- Creation date: Tue May 5 20:11:30 2020. +-- +-- The Metropolis-Hastings-Green ('MHG') algorithm. +-- +-- For example, see Geyer, C. J., Introduction to Markov chain Monte Carlo, In +-- Handbook of Markov Chain Monte Carlo (pp. 45) (2011). CRC press. +module Mcmc.Algorithm.MHG + ( MHG (..), + mhg, + mhgSave, + mhgLoad, + MHGRatio, + mhgAccept, + ) +where + +import Codec.Compression.GZip +import Control.Monad +import Control.Monad.IO.Class +import Data.Aeson +import qualified Data.ByteString.Lazy.Char8 as BL +import Data.Time +import Mcmc.Algorithm +import Mcmc.Chain.Chain +import Mcmc.Chain.Link +import Mcmc.Chain.Save +import Mcmc.Chain.Trace +import Mcmc.Monitor +import Mcmc.Posterior +import Mcmc.Proposal +import Mcmc.Settings +import Numeric.Log +import System.Random.MWC +import Text.Printf +import Prelude hiding (cycle) + +-- | The MHG algorithm. +newtype MHG a = MHG {fromMHG :: Chain a} + +instance ToJSON a => Algorithm (MHG a) where + aName = const "Metropolis-Hastings-Green (MHG)" + aIteration = iteration . fromMHG + aIsInValidState = mhgIsInValidState + aIterate = mhgIterate + aAutoTune = mhgAutoTune + aResetAcceptance = mhgResetAcceptance + aSummarizeCycle = mhgSummarizeCycle + aOpenMonitors = mhgOpenMonitors + aExecuteMonitors = mhgExecuteMonitors + aStdMonitorHeader = mhgStdMonitorHeader + aCloseMonitors = mhgCloseMonitors + aSave = mhgSave + +-- NOTE: IO is required because the trace is mutable. + +-- | Initialize an MHG algorithm. +mhg :: + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + TraceLength -> + InitialState a -> + GenIO -> + IO (MHG a) +mhg pr lh cc mn trLen i0 g = do + -- The trace is a mutable vector and the mutable state needs to be handled by + -- a monad. + tr <- replicateT traceLength l0 + return $ MHG $ Chain Nothing l0 0 tr ac g 0 pr lh cc mn + where + l0 = Link i0 (pr i0) (lh i0) + ac = emptyA $ ccProposals cc + batchMonitorSizes = map getMonitorBatchSize $ mBatches mn + minimumTraceLength = case trLen of + TraceAuto -> 1 + TraceMinimum n -> n + traceLength = maximum $ minimumTraceLength : batchMonitorSizes + +mhgFn :: AnalysisName -> FilePath +mhgFn (AnalysisName nm) = nm ++ ".mcmc.mhg" + +-- | Save an MHG algorithm. +mhgSave :: + ToJSON a => + AnalysisName -> + MHG a -> + IO () +mhgSave nm (MHG c) = do + savedChain <- toSavedChain c + BL.writeFile (mhgFn nm) $ compress $ encode savedChain + +-- | Load an MHG algorithm. +-- +-- See 'Mcmc.Mcmc.mcmcContinue'. +mhgLoad :: + FromJSON a => + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + AnalysisName -> + IO (MHG a) +mhgLoad pr lh cc mn nm = do + savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm) + chain <- either error (fromSavedChain pr lh cc mn) savedChain + return $ MHG chain + +-- | MHG ratios are stored in log domain. +type MHGRatio = Log Double + +-- The MHG ratio. This implementation has the following properties: +-- +-- - The ratio is 'Infinity' if fX is zero. In this case, the proposal is always +-- accepted. +-- +-- - The ratio 'NaN' if (fY or q or j) and fX are zero. In this case, the +-- proposal is always rejected. +-- +-- This means that a chain in a state with posterior probability zero (fX=0) can +-- only move if a state with non-zero posterior probability is proposed. +-- Otherwise it is stuck. Therefore, I print a warning when the posterior +-- probability is zero in the beginning of the MCMC run. This is probably not +-- the best behavior, but see below. +-- +-- There is a discrepancy between authors saying that one should (a) always +-- accept the new state when the current posterior is zero (Chapter 4 of [1], +-- [2]), or (b) almost surely reject the proposal when either fY or q are zero +-- (Chapter 1 of [1]). +-- +-- Since I trust the author of Chapter 1 (Charles Geyer) I choose to follow +-- option (b). However, Option (a) is more user-friendly. +-- +-- [1] Handbook of markov chain monte carlo (2011), CRC press. +-- +-- [2] Dellaportas, P., & Roberts, G. O., An introduction to mcmc, Lecture Notes +-- in Statistics, (), 1–41 (2003). +-- http://dx.doi.org/10.1007/978-0-387-21811-3_1. +mhgRatio :: Posterior -> Posterior -> KernelRatio -> Jacobian -> MHGRatio +-- q = qYX / qXY * jXY; see 'ProposalSimple'. +-- j = Jacobian. +mhgRatio fX fY q j = fY / fX * q * j +{-# INLINE mhgRatio #-} + +-- | Accept or reject a proposal with given MHG ratio? +mhgAccept :: MHGRatio -> GenIO -> IO Bool +mhgAccept r g + | ln r >= 0.0 = return True + | otherwise = do + b <- uniform g + return $ b < exp (ln r) + +mhgPropose :: MHG a -> Proposal a -> IO (MHG a) +mhgPropose (MHG c) p = do + -- 1. Sample new state. + (!y, !q, !j) <- liftIO $ s x g + -- 2. Calculate Metropolis-Hastings-Green ratio. + let !pY = pF y + !lY = lF y + !r = mhgRatio (pX * lX) (pY * lY) q j + -- 3. Accept or reject. + -- if ln r >= 0.0 + -- then do + -- let !ac' = pushA p True ac + -- return $ MHG $ c {link = Link y pY lY, acceptance = ac'} + -- else do + -- b <- uniform g + -- if b < exp (ln r) + -- then do + -- let !ac' = pushA p True ac + -- return $ MHG $ c {link = Link y pY lY, acceptance = ac'} + -- else do + -- let !ac' = pushA p False ac + -- return $ MHG $ c {acceptance = pushA p False ac'} + accept <- mhgAccept r g + if accept + then do + let !ac' = pushA p True ac + return $ MHG $ c {link = Link y pY lY, acceptance = ac'} + else do + let !ac' = pushA p False ac + return $ MHG $ c {acceptance = pushA p False ac'} + where + s = prSimple p + (Link x pX lX) = link c + pF = priorFunction c + lF = likelihoodFunction c + ac = acceptance c + g = generator c + +mhgPush :: MHG a -> IO (MHG a) +mhgPush (MHG c) = do + t' <- pushT i t + return $ MHG c {trace = t', iteration = succ n} + where + i = link c + t = trace c + n = iteration c + +mhgIsInValidState :: MHG a -> Bool +mhgIsInValidState a = ((p * l) == 0) || (p * l == (0 / 0)) + where + x = link $ fromMHG a + p = prior x + l = likelihood x + +-- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG +-- algorithm is just inherently sequential. Parallelization can be achieved by +-- having parallel prior and/or likelihood functions, or by using algorithms +-- running parallel chains such as 'MC3'. +mhgIterate :: ParallelizationMode -> MHG a -> IO (MHG a) +mhgIterate _ a = do + ps <- prepareProposals cc g + a' <- foldM mhgPropose a ps + mhgPush a' + where + c = fromMHG a + cc = cycle c + g = generator c + +mhgAutoTune :: MHG a -> MHG a +mhgAutoTune (MHG c) = MHG $ c {cycle = autoTuneCycle ac cc} + where + ac = acceptance c + cc = cycle c + +mhgResetAcceptance :: MHG a -> MHG a +mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac} + where + ac = acceptance c + +mhgSummarizeCycle :: MHG a -> BL.ByteString +mhgSummarizeCycle (MHG c) = summarizeCycle ac cc + where + cc = cycle c + ac = acceptance c + +mhgOpenMonitors :: AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a) +mhgOpenMonitors nm em (MHG c) = do + m' <- mOpen pre suf em m + return $ MHG c {monitor = m'} + where + m = monitor c + pre = fromAnalysisName nm + suf = maybe "" (printf "%02d") $ chainId c + +mhgExecuteMonitors :: + Verbosity -> + -- Starting time. + UTCTime -> + -- Total number of iterations. + Int -> + MHG a -> + IO (Maybe BL.ByteString) +mhgExecuteMonitors vb t0 iTotal (MHG c) = mExec vb i i0 t0 tr iTotal m + where + i = iteration c + i0 = start c + tr = trace c + m = monitor c + +mhgStdMonitorHeader :: MHG a -> BL.ByteString +mhgStdMonitorHeader (MHG c) = msHeader (mStdOut $ monitor c) + +mhgCloseMonitors :: MHG a -> IO (MHG a) +mhgCloseMonitors (MHG c) = do + m' <- mClose m + return $ MHG $ c {monitor = m'} + where + m = monitor c diff --git a/src/Mcmc/Algorithm/Metropolis.hs b/src/Mcmc/Algorithm/Metropolis.hs deleted file mode 100644 --- a/src/Mcmc/Algorithm/Metropolis.hs +++ /dev/null @@ -1,254 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE OverloadedStrings #-} - --- | --- Module : Mcmc.Algorithm.Metropolis --- Description : Metropolis-Hastings-Green algorithm --- Copyright : (c) Dominik Schrempf 2020 --- License : GPL-3.0-or-later --- --- Maintainer : dominik.schrempf@gmail.com --- Stability : unstable --- Portability : portable --- --- Creation date: Tue May 5 20:11:30 2020. --- --- The Metropolis-Hastings-Green ('MHG') algorithm. --- --- For example, see Geyer, C. J., Introduction to Markov chain Monte Carlo, In --- Handbook of Markov Chain Monte Carlo (pp. 45) (2011). CRC press. -module Mcmc.Algorithm.Metropolis - ( MHG (..), - mhg, - mhgSave, - mhgLoad, - mhgAccept, - ) -where - -import Codec.Compression.GZip -import Control.Monad -import Control.Monad.IO.Class -import Data.Aeson -import qualified Data.ByteString.Lazy.Char8 as BL -import Data.Time -import Mcmc.Algorithm -import Mcmc.Chain.Chain -import Mcmc.Chain.Link -import Mcmc.Chain.Save -import Mcmc.Chain.Trace -import Mcmc.Monitor -import Mcmc.Proposal -import Mcmc.Settings -import Numeric.Log -import System.Random.MWC -import Text.Printf -import Prelude hiding (cycle) - --- | The MHG algorithm. -newtype MHG a = MHG {fromMHG :: Chain a} - -instance ToJSON a => Algorithm (MHG a) where - aName = const "Metropolis-Hastings-Green (MHG)" - aIteration = iteration . fromMHG - aIterate = mhgIterate - aAutoTune = mhgAutoTune - aResetAcceptance = mhgResetAcceptance - aSummarizeCycle = mhgSummarizeCycle - aOpenMonitors = mhgOpenMonitors - aExecuteMonitors = mhgExecuteMonitors - aStdMonitorHeader = mhgStdMonitorHeader - aCloseMonitors = mhgCloseMonitors - aSave = mhgSave - --- NOTE: IO is required because the trace is mutable. - --- | Initialize an MHG algorithm. -mhg :: - PriorFunction a -> - LikelihoodFunction a -> - Cycle a -> - Monitor a -> - -- | The initial state in the state space @a@. - a -> - -- | A source of randomness. For reproducible runs, make sure to use - -- generators with the same seed. - GenIO -> - IO (MHG a) -mhg pr lh cc mn i0 g = do - -- The trace is a mutable vector and the mutable state needs to be handled by - -- a monad. - tr <- replicateT traceLength l0 - return $ MHG $ Chain 0 l0 0 tr ac g 0 pr lh cc mn - where - l0 = Link i0 (pr i0) (lh i0) - ac = emptyA $ ccProposals cc - batchMonitorSizes = map getMonitorBatchSize $ mBatches mn - traceLength = maximum $ 1 : batchMonitorSizes - -mhgFn :: AnalysisName -> FilePath -mhgFn (AnalysisName nm) = nm ++ ".mhg" - --- | Save an MHG algorithm. -mhgSave :: - ToJSON a => - AnalysisName -> - MHG a -> - IO () -mhgSave nm (MHG c) = do - savedChain <- toSavedChain c - BL.writeFile (mhgFn nm) $ compress $ encode savedChain - --- | Load an MHG algorithm. --- --- See 'Mcmc.Mcmc.mcmcContinue'. -mhgLoad :: - FromJSON a => - PriorFunction a -> - LikelihoodFunction a -> - Cycle a -> - Monitor a -> - AnalysisName -> - IO (MHG a) -mhgLoad pr lh cc mn nm = do - savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm) - chain <- either error (fromSavedChain pr lh cc mn) savedChain - return $ MHG chain - --- The MHG ratio. --- --- 'Infinity' if fX is zero. In this case, the proposal is always accepted. --- --- 'NaN' if (fY or q) and fX are zero. In this case, the proposal is always --- rejected. - --- There is a discrepancy between authors saying that one should (a) always --- accept the new state when the current posterior is zero (Chapter 4 of the --- Handbook of Markov Chain Monte Carlo), or (b) almost surely reject the --- proposal when either fY or q are zero (Chapter 1). Since I trust the author --- of Chapter 1 (Charles Geyer) I choose to follow option (b). -mhgRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double --- q = qYX / qXY * jXY; see 'ProposalSimple'. --- j = Jacobian. -mhgRatio fX fY q j = fY / fX * q * j -{-# INLINE mhgRatio #-} - --- | Accept or reject a proposal with given MHG ratio? -mhgAccept :: Log Double -> GenIO -> IO Bool -mhgAccept r g - | ln r >= 0.0 = return True - | otherwise = do - b <- uniform g - return $ b < exp (ln r) - -mhgPropose :: MHG a -> Proposal a -> IO (MHG a) -mhgPropose (MHG c) p = do - -- 1. Sample new state. - (!y, !q, !j) <- liftIO $ s x g - -- 2. Calculate Metropolis-Hastings-Green ratio. - let !pY = pF y - !lY = lF y - !r = mhgRatio (pX * lX) (pY * lY) q j - -- 3. Accept or reject. - -- if ln r >= 0.0 - -- then do - -- let !ac' = pushA p True ac - -- return $ MHG $ c {link = Link y pY lY, acceptance = ac'} - -- else do - -- b <- uniform g - -- if b < exp (ln r) - -- then do - -- let !ac' = pushA p True ac - -- return $ MHG $ c {link = Link y pY lY, acceptance = ac'} - -- else do - -- let !ac' = pushA p False ac - -- return $ MHG $ c {acceptance = pushA p False ac'} - accept <- mhgAccept r g - if accept - then do - let !ac' = pushA p True ac - return $ MHG $ c {link = Link y pY lY, acceptance = ac'} - else do - let !ac' = pushA p False ac - return $ MHG $ c {acceptance = pushA p False ac'} - where - s = pSimple p - (Link x pX lX) = link c - pF = priorFunction c - lF = likelihoodFunction c - ac = acceptance c - g = generator c - -mhgPush :: MHG a -> IO (MHG a) -mhgPush (MHG c) = do - t' <- pushT i t - return $ MHG c {trace = t', iteration = succ n} - where - i = link c - t = trace c - n = iteration c - --- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG --- algorithm is just inherently sequential. Parallelization can be achieved by --- having parallel prior and/or likelihood functions, or by using algorithms --- running parallel chains such as 'MC3'. -mhgIterate :: ParallelizationMode -> MHG a -> IO (MHG a) -mhgIterate _ a = do - ps <- orderProposals cc g - a' <- foldM mhgPropose a ps - mhgPush a' - where - c = fromMHG a - cc = cycle c - g = generator c - -mhgAutoTune :: MHG a -> MHG a -mhgAutoTune (MHG c) = MHG $ c {cycle = autoTuneCycle ac cc} - where - ac = acceptance c - cc = cycle c - -mhgResetAcceptance :: MHG a -> MHG a -mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac} - where - ac = acceptance c - -mhgSummarizeCycle :: MHG a -> BL.ByteString -mhgSummarizeCycle (MHG c) = summarizeCycle ac cc - where - cc = cycle c - ac = acceptance c - -mhgOpenMonitors :: AnalysisName -> ExecutionMode -> MHG a -> IO (MHG a) -mhgOpenMonitors nm em (MHG c) = do - m' <- mOpen pre suf em m - return $ MHG c {monitor = m'} - where - m = monitor c - pre = fromAnalysisName nm - suf = printf "%02d" $ chainId c - -mhgExecuteMonitors :: - Verbosity -> - -- Starting time. - UTCTime -> - -- Total number of iterations. - Int -> - MHG a -> - IO (Maybe BL.ByteString) -mhgExecuteMonitors vb t0 iTotal (MHG c) = mExec vb i i0 t0 tr iTotal m - where - i = iteration c - i0 = start c - tr = trace c - m = monitor c - -mhgStdMonitorHeader :: MHG a -> BL.ByteString -mhgStdMonitorHeader (MHG c) = msHeader (mStdOut $ monitor c) - -mhgCloseMonitors :: MHG a -> IO (MHG a) -mhgCloseMonitors (MHG c) = do - m' <- mClose m - return $ MHG $ c {monitor = m'} - where - m = monitor c diff --git a/src/Mcmc/Chain/Chain.hs b/src/Mcmc/Chain/Chain.hs --- a/src/Mcmc/Chain/Chain.hs +++ b/src/Mcmc/Chain/Chain.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Chain.Chain -- Description : Simple representation of a Markov chain --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -14,35 +14,29 @@ noPrior, LikelihoodFunction, noLikelihood, + InitialState, Chain (..), ) where --- Note: It is not necessary to add another type @b@ to store supplementary +-- NOTE: It is not necessary to add another type @b@ to store supplementary -- information about the chain. The information can just be stored in @a@ -- equally well. +-- TODO: First class parameters. Make a type class for parameter types (name, +-- lens, proposals, monitors). + import Mcmc.Chain.Link import Mcmc.Chain.Trace +import Mcmc.Likelihood import Mcmc.Monitor +import Mcmc.Prior import Mcmc.Proposal -import Numeric.Log import System.Random.MWC hiding (save) import Prelude hiding (cycle) --- | Prior function. -type PriorFunction a = a -> Log Double - --- | Flat prior function. Useful for testing and debugging. -noPrior :: PriorFunction a -noPrior = const 1.0 - --- | Likelihood function. -type LikelihoodFunction a = a -> Log Double - --- | Flat likelihood function. Useful for testing and debugging. -noLikelihood :: LikelihoodFunction a -noLikelihood = const 1.0 +-- | Type synonym to indicate the initial state. +type InitialState a = a -- | The chain contains all information to run an MCMC sampler. -- @@ -60,8 +54,9 @@ -- The 'Mcmc.Environment.Environment' of the chain is not stored externally. data Chain a = Chain { -- Variables; saved. + -- | Chain index; useful if more chains are run. - chainId :: Int, + chainId :: Maybe Int, -- | The current 'Link' of the chain combines the current state and the -- current likelihood. The link is updated after a proposal has been -- executed. diff --git a/src/Mcmc/Chain/Link.hs b/src/Mcmc/Chain/Link.hs --- a/src/Mcmc/Chain/Link.hs +++ b/src/Mcmc/Chain/Link.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Chain.Link -- Description : The state combined with auxiliary variables --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -18,6 +18,8 @@ import Data.Aeson import Data.Aeson.Types +import Mcmc.Likelihood +import Mcmc.Prior import Numeric.Log -- | Link of a Markov chain. For reasons of computational efficiency, each state @@ -25,10 +27,10 @@ data Link a = Link { -- | The current state in the state space @a@. state :: a, - -- | The current prior. - prior :: Log Double, - -- | The current likelihood. - likelihood :: Log Double + -- | The current prior value. + prior :: Prior, + -- | The current likelihood value. + likelihood :: Likelihood } deriving (Eq, Ord, Show, Read) diff --git a/src/Mcmc/Chain/Save.hs b/src/Mcmc/Chain/Save.hs --- a/src/Mcmc/Chain/Save.hs +++ b/src/Mcmc/Chain/Save.hs @@ -1,13 +1,6 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-} --- | --- Module : Mcmc.Chain.Save --- Description : Save and load a Markov chain --- Copyright : (c) Dominik Schrempf, 2020 --- License : GPL-3.0-or-later --- - -- Maintainer : dominik.schrempf@gmail.com -- Stability : unstable -- Portability : portable @@ -17,6 +10,12 @@ -- Save and load chains. It is easy to save and restore the current state and -- likelihood (or the trace), but it is not feasible to store all the proposals -- and so on, so they have to be provided again when continuing a run. + +-- | +-- Module : Mcmc.Chain.Save +-- Description : Save and load a Markov chain +-- Copyright : (c) Dominik Schrempf, 2021 +-- License : GPL-3.0-or-later module Mcmc.Chain.Save ( SavedChain (..), toSavedChain, @@ -30,6 +29,7 @@ import Data.List hiding (cycle) import qualified Data.Map as M import Data.Maybe +import qualified Data.Stack.Circular as C import qualified Data.Vector as VB import qualified Data.Vector.Unboxed as VU import Data.Word @@ -40,20 +40,18 @@ import Mcmc.Monitor import Mcmc.Proposal import Prelude hiding (cycle) -import qualified Data.Stack.Circular as C -- | Storable values of a Markov chain. -- -- See 'toSavedChain'. data SavedChain a = SavedChain - { - savedId :: Int, + { savedId :: Maybe Int, savedLink :: Link a, savedIteration :: Int, savedTrace :: C.Stack VB.Vector (Link a), savedAcceptance :: Acceptance Int, savedSeed :: VU.Vector Word32, - savedTuningParameters :: [Maybe Double] + savedTuningParameters :: [Maybe TuningParameter] } deriving (Eq, Read, Show) @@ -70,7 +68,7 @@ where ps = ccProposals cc ac' = transformKeysA ps [0 ..] ac - ts = [fmap tParam mt | mt <- map pTuner ps] + ts = [fmap tParam mt | mt <- map prTuner ps] -- | Load a saved chain. -- @@ -87,13 +85,19 @@ IO (Chain a) fromSavedChain pr lh cc mn (SavedChain ci it i tr ac' g' ts) | pr (state it) /= prior it = - error "fromSave: Provided prior function does not match the saved prior." + let msg = + unlines + [ "fromSave: Provided prior function does not match the saved prior.", + "fromSave: Current prior:" <> show (prior it) <> ".", + "fromSave: Given prior:" <> show (pr $ state it) <> "." + ] + in error msg | lh (state it) /= likelihood it = error "fromSave: Provided likelihood function does not match the saved likelihood." | otherwise = do - g <- loadGen g' - tr' <- thawT tr - return $ Chain ci it i tr' ac g i pr lh cc' mn + g <- loadGen g' + tr' <- thawT tr + return $ Chain ci it i tr' ac g i pr lh cc' mn where ac = transformKeysA [0 ..] (ccProposals cc) ac' getTuningF mt = case mt of diff --git a/src/Mcmc/Chain/Trace.hs b/src/Mcmc/Chain/Trace.hs --- a/src/Mcmc/Chain/Trace.hs +++ b/src/Mcmc/Chain/Trace.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Chain.Trace -- Description : History of a Markov chain --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com diff --git a/src/Mcmc/Environment.hs b/src/Mcmc/Environment.hs --- a/src/Mcmc/Environment.hs +++ b/src/Mcmc/Environment.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Environment -- Description : Runtime environment --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -12,38 +12,75 @@ module Mcmc.Environment ( Environment (..), initializeEnvironment, + closeEnvironment, ) where +import Control.Concurrent.MVar +import Control.Monad import Data.Time.Clock +import Mcmc.Logger import Mcmc.Settings import System.IO -- | The environment of an MCMC run. -data Environment = Environment - { settings :: Settings, - -- | We have to use 'Maybe' here, because we do not want to open any log - -- file when being 'Quiet'. - logHandle :: Maybe Handle, +data Environment s = Environment + { settings :: s, + -- | List will be empty if using 'Quiet'. If 'LogStdOutAndFile' is used + -- 'logHandles' contains two handles to the standard output and the log + -- file. + logHandles :: [Handle], + -- | MVar blocking output. + outLock :: MVar (), -- | Used to calculate the ETA. startingTime :: UTCTime } - deriving (Eq, Show) + deriving (Eq) +instance HasExecutionMode s => HasExecutionMode (Environment s) where + getExecutionMode = getExecutionMode . settings + +instance HasLock (Environment s) where + getLock = outLock + +instance HasLogHandles (Environment s) where + getLogHandles = logHandles + +instance HasStartingTime (Environment s) where + getStartingTime = startingTime + +instance HasLogMode s => HasLogMode (Environment s) where + getLogMode = getLogMode . settings + +instance HasVerbosity s => HasVerbosity (Environment s) where + getVerbosity = getVerbosity . settings + -- | Initialize the environment. -- -- Open log file, get current time. initializeEnvironment :: - Settings -> - IO Environment + (HasAnalysisName s, HasExecutionMode s, HasLogMode s, HasVerbosity s) => + s -> + IO (Environment s) initializeEnvironment s = do t <- getCurrentTime - mh <- case sVerbosity s of - Quiet -> return Nothing - _ -> do + mh <- case (getLogMode s, getVerbosity s) of + (_, Quiet) -> return [] + (LogStdOutAndFile, _) -> do h <- openWithExecutionMode em fn - return $ Just h - return $ Environment s mh t + return [stdout, h] + (LogFileOnly, _) -> do + h <- openWithExecutionMode em fn + return [h] + (LogStdOutOnly, _) -> return [stdout] + lock <- newMVar () + return $ Environment s mh lock t where - fn = fromAnalysisName (sAnalysisName s) ++ ".log" - em = sExecutionMode s + fn = fromAnalysisName (getAnalysisName s) ++ ".mcmc.log" + em = getExecutionMode s + +-- | Close file handles. +closeEnvironment :: Environment s -> IO () +closeEnvironment e = forM_ hs hClose + where + hs = filter (/= stdout) $ logHandles e diff --git a/src/Mcmc/Internal/ByteString.hs b/src/Mcmc/Internal/ByteString.hs --- a/src/Mcmc/Internal/ByteString.hs +++ b/src/Mcmc/Internal/ByteString.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Internal.ByteString -- Description : ByteString tools --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com diff --git a/src/Mcmc/Internal/Random.hs b/src/Mcmc/Internal/Random.hs --- a/src/Mcmc/Internal/Random.hs +++ b/src/Mcmc/Internal/Random.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Internal.Random -- Description : Tools for random calculations --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com diff --git a/src/Mcmc/Internal/Shuffle.hs b/src/Mcmc/Internal/Shuffle.hs --- a/src/Mcmc/Internal/Shuffle.hs +++ b/src/Mcmc/Internal/Shuffle.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Internal.Shuffle -- Description : Shuffle a list --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com diff --git a/src/Mcmc/Likelihood.hs b/src/Mcmc/Likelihood.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Likelihood.hs @@ -0,0 +1,29 @@ +-- | +-- Module : Mcmc.Likelihood +-- Description : Types and convenience functions for computing likelihoods +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Wed Mar 3 11:39:04 2021. +module Mcmc.Likelihood + ( Likelihood, + LikelihoodFunction, + noLikelihood, + ) +where + +import Numeric.Log + +-- | Likelihood values are stored in log domain. +type Likelihood = Log Double + +-- | Likelihood function. +type LikelihoodFunction a = a -> Log Double + +-- | Flat likelihood function. Useful for testing and debugging. +noLikelihood :: LikelihoodFunction a +noLikelihood = const 1.0 diff --git a/src/Mcmc/Logger.hs b/src/Mcmc/Logger.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Logger.hs @@ -0,0 +1,170 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} + +-- | +-- Module : Mcmc.Logger +-- Description : Minimal monad logger +-- Copyright : (c) Dominik Schrempf 2021 +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : unstable +-- Portability : portable +-- +-- Creation date: Tue Jan 12 09:03:04 2021. +module Mcmc.Logger + ( LogMode (..), + Verbosity (..), + HasLock (..), + HasLogHandles (..), + HasStartingTime (..), + HasLogMode (..), + HasVerbosity (..), + Logger, + logOutB, + logDebugB, + logDebugS, + logWarnB, + logWarnS, + logInfoB, + logInfoS, + logInfoHeader, + logInfoStartingTime, + logInfoEndTime, + ) +where + +import Control.Concurrent.MVar +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Trans.Reader +import Data.Aeson +import Data.Aeson.TH +import qualified Data.ByteString.Lazy.Char8 as BL +import Data.Time.Clock +import Mcmc.Monitor.Time +import System.IO +import Paths_mcmc (version) +import Data.Version (showVersion) + +-- TODO: Combine LogMode and Verbosity to: +-- +-- data Verbosity = Quiet | Warn LogMode | Info LogMode | Debug LogMode + +-- | Define where the log output should be directed to. +data LogMode = LogStdOutAndFile | LogStdOutOnly | LogFileOnly + deriving (Eq, Read, Show) + +$(deriveJSON defaultOptions ''LogMode) + +-- | Not much to say here. +data Verbosity = Quiet | Warn | Info | Debug + deriving (Eq, Ord, Read, Show) + +$(deriveJSON defaultOptions ''Verbosity) + +-- | Types with an output lock for concurrent output. +class HasLock e where + getLock :: e -> MVar () + +-- | Types with logging information. +class HasLogHandles e where + getLogHandles :: e -> [Handle] + +-- | Types with starting time. +class HasStartingTime s where + getStartingTime :: s -> UTCTime + +-- | Types with a log mode. +class HasLogMode s where + getLogMode :: s -> LogMode + +-- | Types with verbosity. +class HasVerbosity s where + getVerbosity :: s -> Verbosity + +-- | Reader transformer used for logging to a file and to standard output. +type Logger e a = ReaderT e IO a + +msgPrepare :: BL.ByteString -> BL.ByteString -> BL.ByteString +msgPrepare pref msg = BL.intercalate "\n" $ map (BL.append pref) $ BL.lines msg + +-- Make sure that concurrent output is not scrambled. +atomicAction :: HasLock e => IO () -> Logger e () +atomicAction a = do + l <- reader getLock + liftIO $ withMVar l (const a) + +-- | Write to standard output and maybe to log file. +logOutB :: + (HasLogHandles e, HasLock e) => + -- | Prefix. + BL.ByteString -> + -- | Message. + BL.ByteString -> + Logger e () +logOutB pref msg = do + hs <- reader getLogHandles + mapM_ (atomicAction . (`BL.hPutStrLn` msg')) hs + where + msg' = msgPrepare pref msg + +-- Perform debug action. +logDebugA :: (HasLock e, HasLogHandles e, HasVerbosity e) => Logger e () -> Logger e () +logDebugA a = reader getVerbosity >>= \v -> when (v >= Debug) a + +-- | Log debug message. +logDebugB :: (HasLock e, HasLogHandles e, HasVerbosity e) => BL.ByteString -> Logger e () +logDebugB = logDebugA . logOutB "D: " + +-- | Log debug message. +logDebugS :: (HasLock e, HasLogHandles e, HasVerbosity e) => String -> Logger e () +logDebugS = logDebugB . BL.pack + +-- Perform warning action. +logWarnA :: (HasLogHandles e, HasVerbosity e) => Logger e () -> Logger e () +logWarnA a = reader getVerbosity >>= \v -> when (v >= Warn) a + +-- | Log warning message. +logWarnB :: (HasLock e, HasLogHandles e, HasVerbosity e) => BL.ByteString -> Logger e () +logWarnB = logWarnA . logOutB "W: " + +-- | Log warning message. +logWarnS :: (HasLock e, HasLogHandles e, HasVerbosity e) => String -> Logger e () +logWarnS = logWarnB . BL.pack + +-- Perform info action. +logInfoA :: (HasLogHandles e, HasVerbosity e) => Logger e () -> Logger e () +logInfoA a = reader getVerbosity >>= \v -> when (v >= Info) a + +-- | Log info message. +logInfoB :: (HasLock e, HasLogHandles e, HasVerbosity e) => BL.ByteString -> Logger e () +logInfoB = logInfoA . logOutB " " + +-- | Log info message. +logInfoS :: (HasLock e, HasLogHandles e, HasVerbosity e) => String -> Logger e () +logInfoS = logInfoB . BL.pack + +-- | Log info header. +logInfoHeader :: (HasLock e, HasLogHandles e, HasVerbosity e) => Logger e () +logInfoHeader = do + logInfoS (replicate 70 '-') + logInfoS ("MCMC sampler; version " ++ showVersion version <> ".") + logInfoS "Developed by: Dominik Schrempf." + logInfoS "License: GPL-3.0-or-later." + logInfoS (replicate 70 '-') + +-- | Log starting time. +logInfoStartingTime :: (HasLock e, HasLogHandles e, HasStartingTime e, HasVerbosity e) => Logger e () +logInfoStartingTime = do + ti <- reader getStartingTime + logInfoS $ "Starting time: " <> renderTime ti + +-- | Log end time. +logInfoEndTime :: (HasLock e, HasLogHandles e, HasStartingTime e, HasVerbosity e) => Logger e () +logInfoEndTime = do + ti <- reader getStartingTime + te <- liftIO getCurrentTime + let dt = te `diffUTCTime` ti + logInfoB $ "Wall clock run time: " <> renderDuration dt <> "." + logInfoS $ "End time: " <> renderTime te diff --git a/src/Mcmc/MarginalLikelihood.hs b/src/Mcmc/MarginalLikelihood.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/MarginalLikelihood.hs @@ -0,0 +1,404 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE OverloadedStrings #-} + +-- | +-- Module : Mcmc.MarginalLikelihood +-- Description : Calculate the marginal likelihood +-- Copyright : (c) Dominik Schrempf 2021 +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : unstable +-- Portability : portable +-- +-- Creation date: Mon Jan 11 16:34:18 2021. +module Mcmc.MarginalLikelihood + ( MarginalLikelihood, + NPoints (..), + MLAlgorithm (..), + MLSettings (..), + marginalLikelihood, + ) +where + +import Control.Monad +import Control.Monad.IO.Class +import qualified Control.Monad.Parallel as P +import Control.Monad.Trans.Reader +import Data.Aeson +import Data.List hiding (cycle) +import qualified Data.Map.Strict as M +import qualified Data.Vector as VB +import qualified Data.Vector.Unboxed as VU +import Mcmc.Algorithm.MHG +import Mcmc.Chain.Chain +import Mcmc.Chain.Link +import Mcmc.Chain.Trace +import Mcmc.Environment +import Mcmc.Internal.Random +import Mcmc.Likelihood +import Mcmc.Logger +import Mcmc.Mcmc +import Mcmc.Monitor +import Mcmc.Proposal +import Mcmc.Settings +import Numeric.Log hiding (sum) +import System.Directory +import System.Random.MWC +import Text.Printf +import Text.Show.Pretty +import Prelude hiding (cycle) + +-- | Marginal likelihood values are stored in log domain. +type MarginalLikelihood = Log Double + +-- Reciprocal temperature value traversed along the path integral. +type Point = Double + +-- | The number of points used to approximate the path integral. +newtype NPoints = NPoints {fromNPoints :: Int} + deriving (Eq, Read, Show) + +-- | Algorithms to calculate the marginal likelihood. +data MLAlgorithm + = -- | Use a classical path integral. Also known as thermodynamic integration. + -- In particular, /Annealing-Melting Integration/ is used. + -- + -- See Lartillot, N., & Philippe, H., Computing Bayes Factors Using + -- Thermodynamic Integration, Systematic Biology, 55(2), 195–207 (2006). + -- http://dx.doi.org/10.1080/10635150500433722 + ThermodynamicIntegration + | -- | Use stepping stone sampling. + -- + -- See Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M., Improving + -- marginal likelihood estimation for Bayesian phylogenetic model selection, + -- Systematic Biology, 60(2), 150–160 (2010). + -- http://dx.doi.org/10.1093/sysbio/syq085 + -- + -- Or Fan, Y., Wu, R., Chen, M., Kuo, L., & Lewis, P. O., Choosing among + -- partition models in bayesian phylogenetics, Molecular Biology and + -- Evolution, 28(1), 523–532 (2010). http://dx.doi.org/10.1093/molbev/msq224 + SteppingStoneSampling + deriving (Eq, Read, Show) + +-- | Settings of the marginal likelihood estimation. +data MLSettings = MLSettings + { mlAnalysisName :: AnalysisName, + mlAlgorithm :: MLAlgorithm, + mlNPoints :: NPoints, + -- | Initial burn in at the starting point of the path. + mlInitialBurnIn :: BurnInSpecification, + -- | Repetitive burn in at each point on the path. + mlPointBurnIn :: BurnInSpecification, + -- | The number of iterations performed at each point. + mlIterations :: Iterations, + mlExecutionMode :: ExecutionMode, + mlLogMode :: LogMode, + mlVerbosity :: Verbosity + } + deriving (Eq, Read, Show) + +instance HasAnalysisName MLSettings where + getAnalysisName = mlAnalysisName + +instance HasExecutionMode MLSettings where + getExecutionMode = mlExecutionMode + +instance HasLogMode MLSettings where + getLogMode = mlLogMode + +instance HasVerbosity MLSettings where + getVerbosity = mlVerbosity + +type ML a = ReaderT (Environment MLSettings) IO a + +-- See 'getPoints'. Alpha=0.3 is the standard choice. +alpha :: Double +alpha = 0.3 + +-- Distribute the points according to a skewed beta distribution with given +-- 'alpha' value. If alpha is below 1.0, more points at lower values, which is +-- desired. It is inconvenient that the reciprocal temperatures are denoted as +-- beta, and we also use the beta distribution :). Don't mix them up! +-- +-- See discussion in Xie, W., Lewis, P. O., Fan, Y., Kuo, L., & Chen, M., +-- Improving marginal likelihood estimation for bayesian phylogenetic model +-- selection, Systematic Biology, 60(2), 150–160 (2010). +-- http://dx.doi.org/10.1093/sysbio/syq085 +-- +-- Or Figure 1 in Höhna, S., Landis, M. J., & Huelsenbeck, J. P., Parallel power +-- posterior analyses for fast computation of marginal likelihoods in +-- phylogenetics (2017). http://dx.doi.org/10.1101/104422 +getPoints :: NPoints -> [Point] +getPoints x = [f i ** (1.0 / alpha) | i <- [0 .. k1]] + where + k = fromNPoints x + k1 = pred k + f j = fromIntegral j / fromIntegral k1 + +sampleAtPoint :: + ToJSON a => + Point -> + Settings -> + LikelihoodFunction a -> + MHG a -> + ML (MHG a) +sampleAtPoint x ss lhf a = do + a'' <- liftIO $ mcmc ss' a' + let ch'' = fromMHG a'' + ac = acceptance ch'' + mAr = sequence $ acceptanceRates ac + logDebugB "sampleAtPoint: Summarize cycle." + logDebugB $ summarizeCycle ac $ cycle ch'' + case mAr of + Nothing -> logWarnB "Some acceptance rates are unavailable. The tuning period may be too small." + Just ar -> do + unless (M.null $ M.filter (<= 0.1) ar) $ logWarnB "Some acceptance rates are below 0.1." + unless (M.null $ M.filter (>= 0.9) ar) $ logWarnB "Some acceptance rates are above 0.9." + return a'' + where + -- For debugging set a proper analysis name. + nm = sAnalysisName ss + getName :: Point -> AnalysisName + getName y = nm <> AnalysisName ("/" <> printf "point%.8f" y) + ss' = ss {sAnalysisName = getName x} + -- Amend the likelihood function. Don't calculate the likelihood when the + -- point is 0.0. + lhf' = if x == 0.0 then const 1.0 else (** Exp (log x)) . lhf + -- Amend the MHG algorithm. + ch = fromMHG a + l = link ch + ch' = + ch + { -- Important: Update the likelihood using the new likelihood function. + link = l {likelihood = lhf' $ state l}, + iteration = 0, + start = 0, + likelihoodFunction = lhf' + } + a' = MHG ch' + +traversePoints :: + ToJSON a => + -- Current point. + Int -> + NPoints -> + [Point] -> + Settings -> + LikelihoodFunction a -> + MHG a -> + -- For each point a vector of obtained likelihoods stored in the log domain. + ML [VU.Vector Likelihood] +traversePoints _ _ [] _ _ _ = return [] +traversePoints i k (b : bs) ss lhf a = do + logInfoS $ "Point " <> show i <> " of " <> show k' <> ": " <> show b <> "." + a' <- sampleAtPoint b ss lhf a + -- Get the links samples at this point. + ls <- liftIO $ takeT n $ trace $ fromMHG a' + -- Extract the likelihoods. + -- + -- NOTE: This could be sped up by mapping (** -b) on the power likelihoods. + -- + -- NOTE: This bang is an important one, because if the lhs are not strictly + -- calculated here, the complete MCMC runs are dragged along before doing so + -- resulting in a severe memory leak. + let !lhs = VU.convert $ VB.map (lhf . state) ls + -- Sample the other points. + lhss <- traversePoints (i + 1) k bs ss lhf a' + return $ lhs : lhss + where + n = fromIterations $ sIterations ss + (NPoints k') = k + +mlRun :: + ToJSON a => + NPoints -> + [Point] -> + ExecutionMode -> + Verbosity -> + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + a -> + GenIO -> + -- For each point a vector of likelihoods stored in log domain. + ML [VU.Vector Likelihood] +mlRun k xs em vb prf lhf cc mn i0 g = do + logDebugB "mlRun: Begin." + s <- reader settings + let nm = mlAnalysisName s + is = mlIterations s + biI = mlInitialBurnIn s + biP = mlPointBurnIn s + -- Only log sub MCMC samplers when debugging. + vb' = case vb of + Debug -> Debug + _ -> Quiet + ssI = Settings nm biI (Iterations 0) em Sequential NoSave LogFileOnly vb' + ssP = Settings nm biP is em Sequential NoSave LogFileOnly vb' + trLen = TraceMinimum $ fromIterations is + logDebugB "mlRun: Initialize MHG algorithm." + a0 <- liftIO $ mhg prf lhf cc mn trLen i0 g + logDebugS $ "mlRun: Perform initial burn in at first point " <> show x0 <> "." + a1 <- sampleAtPoint x0 ssI lhf a0 + logDebugB "mlRun: Traverse points." + traversePoints 1 k xs ssP lhf a1 + where + x0 = head xs + +-- Use lists since the number of points is expected to be low. +integrateSimpsonTriangle :: + -- X values. + [Point] -> + -- Y values. + [Double] -> + -- Integral. + Double +integrateSimpsonTriangle xs ys = 0.5 * go xs ys + where + go (p0 : p1 : ps) (z0 : z1 : zs) = (z0 + z1) * (p1 - p0) + go (p1 : ps) (z1 : zs) + go _ _ = 0 + +tiWrapper :: + ToJSON a => + MLSettings -> + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + a -> + GenIO -> + ML MarginalLikelihood +tiWrapper s prf lhf cc mn i0 g = do + logInfoB "Path integral (thermodynamic integration)." + [g0, g1] <- splitGen 2 g + + -- Parallel execution of both path integrals. + [lhssForward, lhssBackward] <- + P.sequence + [ mlRun k bsForward em vb prf lhf cc mn i0 g0, + mlRun k bsBackward em vb prf lhf cc mn i0 g1 + ] + logInfoEndTime + + logDebugB "tiWrapper: Calculate mean log likelihoods." + -- It is important to average across the log likelihoods here (and not the + -- likelihoods). I am not exactly sure why this is. + let getMeanLogLhs = map (\x -> VU.sum (VU.map ln x) / fromIntegral (VU.length x)) + mlForward = integrateSimpsonTriangle bsForward (getMeanLogLhs lhssForward) + mlBackward = negate $ integrateSimpsonTriangle bsBackward (getMeanLogLhs lhssBackward) + logDebugS $ "tiWrapper: Marginal log likelihood of forward integral: " ++ show mlForward + logDebugS $ "tiWrapper: Marginal log likelihood of backward integral: " ++ show mlBackward + let mean = 0.5 * (mlForward + mlBackward) + logDebugS $ "tiWrapper: The mean is: " ++ show mean + return $ Exp mean + where + k = mlNPoints s + bsForward = getPoints k + bsBackward = reverse bsForward + em = mlExecutionMode s + vb = mlVerbosity s + +-- Helper function to exponentiate log domain values with a double value. +pow' :: Log Double -> Double -> Log Double +pow' x p = Exp $ ln x * p + +-- See Xie2010 p. 153, bottom left. +sssCalculateMarginalLikelihood :: [Point] -> [VU.Vector Likelihood] -> MarginalLikelihood +sssCalculateMarginalLikelihood xs lhss = product $ zipWith3 f xs (tail xs) lhss + where + f :: Point -> Point -> VU.Vector Likelihood -> MarginalLikelihood + -- f beta_{k-1} beta_k lhs_{k-1} + f bkm1 bk lhs = n1 * VU.sum lhsPowered + where + n1 = recip $ fromIntegral $ VU.length lhs + dbeta = bk - bkm1 + lhsPowered = VU.map (`pow'` dbeta) lhs + +-- -- Numerical stability by factoring out lhMax. But no observed +-- -- improvement towards the standard version. +-- +-- f bkm1 bk lhs = n1 * pow' lhMax dbeta * VU.sum lhsNormedPowered +-- where n1 = recip $ fromIntegral $ VU.length lhs +-- lhMax = VU.maximum lhs +-- dbeta = bk - bkm1 +-- lhsNormed = VU.map (/lhMax) lhs +-- lhsNormedPowered = VU.map (`pow'` dbeta) lhsNormed + +-- -- Computation of the log of the marginal likelihood. According to the paper, +-- -- this estimator is biased and I did not observe any improvements compared +-- -- to the direct estimator implemented above. +-- +-- -- See Xie2010 p. 153, top right. +-- sssCalculateMarginalLikelihood' :: [Point] -> [VU.Vector Likelihood] -> MarginalLikelihood +-- sssCalculateMarginalLikelihood' xs lhss = Exp $ sum $ zipWith3 f xs (tail xs) lhss +-- where f :: Point -> Point -> VU.Vector Likelihood -> Double +-- -- f beta_{k-1} beta_k lhs_{k-1} +-- f bkm1 bk lhs = dbeta * llhMax + log (n1 * VU.sum lhsNormedPowered) +-- where dbeta = bk - bkm1 +-- llhMax = ln $ VU.maximum lhs +-- n1 = recip $ fromIntegral $ VU.length lhs +-- llhs = VU.map ln lhs +-- llhsNormed = VU.map (\x -> x - llhMax) llhs +-- lhsNormedPowered = VU.map (\x -> exp $ dbeta * x) llhsNormed +sssWrapper :: + ToJSON a => + MLSettings -> + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + a -> + GenIO -> + ML MarginalLikelihood +sssWrapper s prf lhf cc mn i0 g = do + logInfoB "Stepping stone sampling." + logLhss <- mlRun k bsForward' em vb prf lhf cc mn i0 g + logInfoB "The last point does not need to be sampled with stepping stone sampling." + logDebugB "sssWrapper: Calculate marginal likelihood." + return $ sssCalculateMarginalLikelihood bsForward logLhss + where + k = mlNPoints s + bsForward = getPoints k + bsForward' = init bsForward + em = mlExecutionMode s + vb = mlVerbosity s + +-- | Estimate the marginal likelihood. +marginalLikelihood :: + ToJSON a => + MLSettings -> + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + InitialState a -> + -- | A source of randomness. For reproducible runs, make sure to use + -- generators with the same seed. + GenIO -> + IO MarginalLikelihood +marginalLikelihood s prf lhf cc mn i0 g = do + -- Initialize. + e <- initializeEnvironment s + + when (mlVerbosity s == Debug) $ do + let n = fromAnalysisName $ mlAnalysisName s + createDirectoryIfMissing True n + + -- Run. + runReaderT + ( do + logInfoStartingTime + logInfoB "Estimate marginal likelihood." + logDebugB "marginalLikelihood: The marginal likelihood settings are:" + logDebugS $ ppShow s + val <- case mlAlgorithm s of + ThermodynamicIntegration -> tiWrapper s prf lhf cc mn i0 g + SteppingStoneSampling -> sssWrapper s prf lhf cc mn i0 g + logInfoS $ "Marginal log likelihood: " ++ show (ln val) + -- TODO: Simulation variance. + logInfoS "The simulation variance is not yet available." + return val + ) + e diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs --- a/src/Mcmc/Mcmc.hs +++ b/src/Mcmc/Mcmc.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Mcmc -- Description : Framework for running Markov chain Monte Carlo samplers --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -23,14 +23,10 @@ import Control.Monad import Control.Monad.IO.Class --- import Control.Monad.Trans.RWS.CPS import Control.Monad.Trans.Reader -import qualified Data.ByteString.Lazy.Char8 as BL -import Data.Maybe -import Data.Time.Clock import Mcmc.Algorithm import Mcmc.Environment -import Mcmc.Monitor.Time +import Mcmc.Logger import Mcmc.Settings import System.IO import Text.Show.Pretty @@ -38,79 +34,24 @@ -- The MCMC algorithm has read access to an environment and uses an algorithm -- transforming the state @a@. -type MCMC a = ReaderT Environment IO a - -msgPrepare :: BL.ByteString -> BL.ByteString -> BL.ByteString -msgPrepare pref msg = BL.intercalate "\n" $ map (BL.append pref) $ BL.lines msg - --- Write to standard output and log file. -mcmcOutB :: BL.ByteString -> BL.ByteString -> MCMC () -mcmcOutB pref msg = do - h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> reader logHandle - liftIO $ BL.putStrLn msg' >> BL.hPutStrLn h msg' - where - msg' = msgPrepare pref msg - --- -- Perform warning action. --- mcmcWarnA :: MCMC a () -> MCMC a () --- mcmcWarnA a = reader (verbosity . settings) >>= \v -> when (v >= Warn) a - --- -- Print warning message. --- mcmcWarnB :: BL.ByteString -> MCMC a () --- mcmcWarnB = mcmcWarnA . mcmcOutB . msgPrepare 'W' - --- -- Print warning message. --- mcmcWarnS :: String -> MCMC a () --- mcmcWarnS = mcmcWarnB . BL.pack - --- Perform info action. -mcmcInfoA :: MCMC () -> MCMC () -mcmcInfoA a = reader (sVerbosity . settings) >>= \v -> when (v >= Info) a - --- Print info message. -mcmcInfoB :: BL.ByteString -> MCMC () -mcmcInfoB = mcmcInfoA . mcmcOutB "I: " - --- Print info message. -mcmcInfoS :: String -> MCMC () -mcmcInfoS = mcmcInfoB . BL.pack - --- Perform debug action. -mcmcDebugA :: MCMC () -> MCMC () -mcmcDebugA a = reader (sVerbosity . settings) >>= \v -> when (v == Debug) a - --- Print debug message. -mcmcDebugB :: BL.ByteString -> MCMC () -mcmcDebugB = mcmcDebugA . mcmcOutB "D: " - --- Print debug message. -mcmcDebugS :: String -> MCMC () -mcmcDebugS = mcmcDebugB . BL.pack - -mcmcReportTime :: MCMC () -mcmcReportTime = do - mcmcDebugB "Report time." - ti <- reader startingTime - mcmcInfoS $ "Starting time of MCMC sampler: " <> renderTime ti +type MCMC a = ReaderT (Environment Settings) IO a mcmcExecute :: Algorithm a => a -> MCMC a mcmcExecute a = do - mcmcDebugB "Executing MCMC run." + logDebugB "Executing MCMC run." s <- reader settings a' <- case sExecutionMode s of Fail -> mcmcNewRun a Overwrite -> mcmcNewRun a Continue -> mcmcContinueRun a - mcmcDebugB "Executed MCMC run." + logDebugB "Executed MCMC run." return a' --- Reset acceptance counts. mcmcResetAcceptance :: Algorithm a => a -> MCMC a mcmcResetAcceptance a = do - mcmcDebugB "Reset acceptance rates." + logDebugB "Reset acceptance rates." return $ aResetAcceptance a --- Execute the monitors of the chain. mcmcExecuteMonitors :: Algorithm a => a -> MCMC () mcmcExecuteMonitors a = do e <- ask @@ -119,7 +60,7 @@ t0 = startingTime e iTotal = burnInIterations (sBurnIn s) + fromIterations (sIterations s) mStdLog <- liftIO (aExecuteMonitors vb t0 iTotal a) - forM_ mStdLog (mcmcOutB " ") + forM_ mStdLog (logOutB " ") mcmcIterate :: Algorithm a => Int -> a -> MCMC a mcmcIterate n a @@ -134,30 +75,30 @@ mcmcNewRun :: Algorithm a => a -> MCMC a mcmcNewRun a = do s <- reader settings - mcmcInfoB "Start new MCMC sampler." - mcmcInfoB "Initial state." - mcmcInfoB $ aStdMonitorHeader a + logInfoB "Start new MCMC sampler." + logInfoB "Initial state." + logInfoB $ aStdMonitorHeader a mcmcExecuteMonitors a - mcmcInfoB $ aSummarizeCycle a + when (aIsInValidState a) (logWarnB "The initial state is invalid!") + logInfoB $ aSummarizeCycle a a' <- mcmcBurnIn a - a'' <- mcmcResetAcceptance a' let i = fromIterations $ sIterations s - mcmcInfoS $ "Run chain for " ++ show i ++ " iterations." - mcmcInfoB $ aStdMonitorHeader a'' - mcmcIterate i a'' + logInfoS $ "Run chain for " ++ show i ++ " iterations." + logInfoB $ aStdMonitorHeader a' + mcmcIterate i a' mcmcContinueRun :: Algorithm a => a -> MCMC a mcmcContinueRun a = do s <- reader settings let iTotal = fromIterations (sIterations s) + burnInIterations (sBurnIn s) - mcmcInfoB "Continuation of MCMC sampler." + logInfoB "Continuation of MCMC sampler." let iCurrent = aIteration a - mcmcInfoS $ "Current iteration: " ++ show iCurrent ++ "." - mcmcInfoS $ "Total iterations: " ++ show iTotal ++ "." + logInfoS $ "Current iteration: " ++ show iCurrent ++ "." + logInfoS $ "Total iterations: " ++ show iTotal ++ "." let di = iTotal - iCurrent - mcmcInfoB $ aSummarizeCycle a - mcmcInfoS $ "Run chain for " ++ show di ++ " iterations." - mcmcInfoB $ aStdMonitorHeader a + logInfoB $ aSummarizeCycle a + logInfoS $ "Run chain for " ++ show di ++ " iterations." + logInfoB $ aStdMonitorHeader a mcmcIterate di a mcmcBurnIn :: Algorithm a => a -> MCMC a @@ -165,53 +106,67 @@ s <- reader settings case sBurnIn s of NoBurnIn -> do - mcmcInfoS "No burn in." + logInfoS "No burn in." return a BurnInWithoutAutoTuning n -> do - mcmcInfoS $ "Burn in for " <> show n <> " iterations." - mcmcInfoS "Auto tuning is disabled." - mcmcInfoB $ aStdMonitorHeader a + logInfoS $ "Burn in for " <> show n <> " iterations." + logInfoS "Auto tuning is disabled." + logInfoB $ aStdMonitorHeader a a' <- mcmcIterate n a - mcmcInfoB $ aSummarizeCycle a' - mcmcInfoB "Burn in finished." - return a' + logInfoB $ aSummarizeCycle a' + a'' <- mcmcResetAcceptance a' + logInfoB "Burn in finished." + return a'' BurnInWithAutoTuning n t -> do - mcmcInfoS $ "Burn in for " ++ show n ++ " iterations." - mcmcInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "." - mcmcInfoB $ aStdMonitorHeader a - a' <- mcmcBurnInWithAutoTuning n t a - mcmcInfoB "Burn in finished." + logInfoS $ "Burn in for " ++ show n ++ " iterations." + logInfoS $ "Auto tuning is enabled with a period of " ++ show t ++ "." + logInfoB $ aStdMonitorHeader a + let (m, r) = n `divMod` t + -- Don't add if 0. Because then we auto tune without acceptance counts + -- and get NaNs. + xs = replicate m t <> [r | r > 0] + a' <- mcmcBurnInWithAutoTuning xs a + logInfoB "Burn in finished." return a' + BurnInWithCustomAutoTuning xs -> do + logInfoS $ "Burn in for " ++ show (sum xs) ++ " iterations." + logInfoS $ "Custom auto tuning is enabled with periods " ++ show xs ++ "." + logInfoB $ aStdMonitorHeader a + a' <- mcmcBurnInWithAutoTuning xs a + logInfoB "Burn in finished." + return a' -- Auto tune the proposals. mcmcAutotune :: Algorithm a => a -> MCMC a mcmcAutotune a = do - mcmcDebugB "Auto tune." + logDebugB "Auto tune." return $ aAutoTune a -mcmcBurnInWithAutoTuning :: Algorithm a => Int -> Int -> a -> MCMC a -mcmcBurnInWithAutoTuning b t a - | b > t = do - a' <- mcmcResetAcceptance a - a'' <- mcmcIterate t a' - mcmcDebugB $ aSummarizeCycle a'' - a''' <- mcmcAutotune a'' - mcmcDebugB $ aStdMonitorHeader a'' - mcmcBurnInWithAutoTuning (b - t) t a''' - | otherwise = do - a' <- mcmcResetAcceptance a - a'' <- mcmcIterate b a' - mcmcInfoB $ aSummarizeCycle a'' - mcmcInfoS $ "Acceptance rates calculated over the last " <> show b <> " iterations." - return a'' +mcmcBurnInWithAutoTuning :: Algorithm a => [Int] -> a -> MCMC a +mcmcBurnInWithAutoTuning [] _ = error "mcmcBurnInWithAutoTuning: Empty lisst." +mcmcBurnInWithAutoTuning [x] a = do + -- Last round. + a' <- mcmcIterate x a + a'' <- mcmcAutotune a' + logInfoB $ aSummarizeCycle a'' + logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations." + mcmcResetAcceptance a'' +mcmcBurnInWithAutoTuning (x:xs) a = do + a' <- mcmcIterate x a + a'' <- mcmcAutotune a' + logDebugB $ aSummarizeCycle a'' + logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations." + logDebugB $ aStdMonitorHeader a'' + a''' <- mcmcResetAcceptance a'' + mcmcBurnInWithAutoTuning xs a''' mcmcInitialize :: Algorithm a => a -> MCMC a mcmcInitialize a = do - mcmcInfoS $ aName a ++ " algorithm." + logInfoS $ aName a ++ " algorithm." s <- settings <$> ask - mcmcDebugB "Opening monitors." + logDebugB "Opening monitors." a' <- liftIO $ aOpenMonitors (sAnalysisName s) (sExecutionMode s) a - mcmcDebugB "Monitors opened." + logDebugB "Monitors opened." return a' -- Save the MCMC run. @@ -219,42 +174,41 @@ mcmcSave a = do s <- reader settings case sSaveMode s of - NoSave -> mcmcInfoB "Do not save the MCMC analysis." + NoSave -> logInfoB "Do not save the MCMC analysis." Save -> do - mcmcInfoB "Save settings." + logInfoB "Save settings." liftIO $ settingsSave s let nm = sAnalysisName s - mcmcInfoB "Save compressed MCMC analysis." - mcmcInfoB "For long traces, or complex objects, this may take a while." + logInfoB "Save compressed MCMC analysis." + logInfoB "For long traces, or complex objects, this may take a while." liftIO $ aSave nm a - mcmcInfoB "Markov chain saved." + logInfoB "Markov chain saved." -- Report and finish up. mcmcClose :: Algorithm a => a -> MCMC a mcmcClose a = do - mcmcDebugB "Closing MCMC run." - mcmcInfoB $ aSummarizeCycle a - mcmcInfoS $ aName a ++ " algorithm finished." + logDebugB "Closing MCMC run." + logInfoB $ aSummarizeCycle a + logInfoS $ aName a ++ " algorithm finished." mcmcSave a - ti <- reader startingTime - te <- liftIO getCurrentTime - let dt = te `diffUTCTime` ti - mcmcInfoB $ "Wall clock run time: " <> renderDuration dt <> "." - mcmcInfoS $ "End time: " <> renderTime te + logInfoEndTime a' <- liftIO $ aCloseMonitors a - h <- reader logHandle - liftIO $ forM_ h hClose + e <- ask + liftIO $ closeEnvironment e return a' -- Initialize the run, execute the run, and close the run. mcmcRun :: Algorithm a => a -> MCMC a mcmcRun a = do - mcmcDebugB "The settings are:" - reader settings >>= mcmcDebugS . ppShow + logInfoHeader + -- Debug settings. + logDebugB "The MCMC settings are:" + reader settings >>= logDebugS . ppShow + -- Initialize. a' <- mcmcInitialize a - mcmcReportTime + logInfoStartingTime -- Execute. a'' <- mcmcExecute a' @@ -274,15 +228,15 @@ -- Currently, it is only possible to continue MCMC algorithms that have -- completed successfully. This restriction is necessary, because for parallel -- chains, it is hardly possible to ensure all chains are synchronized when the --- process is killed. +-- process is killed or fails. -- -- See: -- --- - 'Mcmc.Algorithm.Metropolis.mhgLoad' +-- - 'Mcmc.Algorithm.MHG.mhgLoad' -- -- - 'Mcmc.Algorithm.MC3.mc3Load' -mcmcContinue :: Algorithm a => Int -> Settings -> a -> IO a +mcmcContinue :: Algorithm a => Iterations -> Settings -> a -> IO a mcmcContinue dn s = mcmc s' where - n' = Iterations $ fromIterations (sIterations s) + dn + n' = Iterations $ fromIterations (sIterations s) + fromIterations dn s' = s {sIterations = n', sExecutionMode = Continue} diff --git a/src/Mcmc/Monitor.hs b/src/Mcmc/Monitor.hs --- a/src/Mcmc/Monitor.hs +++ b/src/Mcmc/Monitor.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Monitor -- Description : Monitor a Markov chain --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -14,11 +14,14 @@ module Mcmc.Monitor ( -- * Create monitors Monitor (..), + Period, + simpleMonitor, MonitorStdOut, monitorStdOut, msHeader, MonitorFile, monitorFile, + BatchSize, MonitorBatch, monitorBatch, getMonitorBatchSize, @@ -64,21 +67,31 @@ mBatches :: [MonitorBatch a] } +-- | Monitor period. +type Period = Int + +-- | Do not monitor parameters. +-- +-- Monitor prior and likelihood with given period. +simpleMonitor :: Period -> Monitor a +simpleMonitor p + | p < 1 = error "simpleMonitor: Monitor period must be 1 or larger." + | otherwise = + Monitor (MonitorStdOut [] p) [] [] + -- | Monitor to standard output; constructed with 'monitorStdOut'. data MonitorStdOut a = MonitorStdOut { msParams :: [MonitorParameter a], - msPeriod :: Int + msPeriod :: Period } -- | Monitor to standard output. monitorStdOut :: - -- | Instructions about which parameters to log. [MonitorParameter a] -> - -- | Logging period. - Int -> + Period -> MonitorStdOut a monitorStdOut ps p - | p < 1 = error "monitorStdOut: Monitor period has to be 1 or larger." + | p < 1 = error "monitorStdOut: Monitor period must be 1 or larger." | otherwise = MonitorStdOut ps p msIWidth :: Int @@ -151,20 +164,18 @@ { mfName :: String, mfHandle :: Maybe Handle, mfParams :: [MonitorParameter a], - mfPeriod :: Int + mfPeriod :: Period } -- | Monitor parameters to a file. monitorFile :: -- | Name; used as part of the file name. String -> - -- | Instructions about which parameters to log. [MonitorParameter a] -> - -- | Logging period. - Int -> + Period -> MonitorFile a monitorFile n ps p - | p < 1 = error "monitorFile: Monitor period has to be 1 or larger." + | p < 1 = error "monitorFile: Monitor period must be 1 or larger." | otherwise = MonitorFile n Nothing ps p mfRenderRow :: [BL.ByteString] -> BL.ByteString @@ -216,6 +227,9 @@ Just h -> hClose h Nothing -> error $ "mfClose: File was not opened for monitor " <> mfName m <> "." +-- | Batch size. +type BatchSize = Int + -- | Batch monitor to a file. -- -- Calculate summary statistics over the last given number of iterations (batch @@ -224,26 +238,24 @@ { mbName :: String, mbHandle :: Maybe Handle, mbParams :: [MonitorParameterBatch a], - mbSize :: Int + mbSize :: BatchSize } -- | Batch monitor parameters to a file, see 'MonitorBatch'. monitorBatch :: -- | Name; used as part of the file name. String -> - -- | Instructions about how to calculate the summary statistics. [MonitorParameterBatch a] -> - -- | Batch size. - Int -> + BatchSize -> MonitorBatch a -monitorBatch n ps p - | p < 2 = error "monitorBatch: Batch size has to be 2 or larger." - | otherwise = MonitorBatch n Nothing ps p +monitorBatch n ps b + | b < 2 = error "monitorBatch: Batch size must be 2 or larger." + | otherwise = MonitorBatch n Nothing ps b -- | Batch monitor size. -- -- Useful to determine the trace length. -getMonitorBatchSize :: MonitorBatch a -> Int +getMonitorBatchSize :: MonitorBatch a -> BatchSize getMonitorBatchSize = mbSize mbOpen :: String -> String -> ExecutionMode -> MonitorBatch a -> IO (MonitorBatch a) @@ -322,7 +334,6 @@ -- | Execute monitors; print status information to files and return text to be -- printed to standard output and log file. mExec :: - -- | Verbosity Verbosity -> -- | Iteration. Int -> @@ -330,11 +341,9 @@ Int -> -- | Starting time. UTCTime -> - -- | Trace of Markov chain. Trace a -> -- | Total number of iterations; to calculate ETA. Int -> - -- | The monitor. Monitor a -> IO (Maybe BL.ByteString) mExec v i ss st xs j (Monitor s fs bs) = do diff --git a/src/Mcmc/Monitor/Log.hs b/src/Mcmc/Monitor/Log.hs --- a/src/Mcmc/Monitor/Log.hs +++ b/src/Mcmc/Monitor/Log.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Monitor.Log -- Description : Monitor logarithmic values --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com 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 @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Monitor.Parameter -- Description : Monitor parameters --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com 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 @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Monitor.ParameterBatch -- Description : Batch monitor parameters --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com 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 @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Monitor.Time -- Description : Print time related values --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com diff --git a/src/Mcmc/Posterior.hs b/src/Mcmc/Posterior.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Posterior.hs @@ -0,0 +1,20 @@ +-- | +-- Module : Mcmc.Posterior +-- Description : Types for posterior values and functions +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Fri May 28 12:26:35 2021. +module Mcmc.Posterior + ( Posterior, + ) +where + +import Numeric.Log + +-- | Posterior values are stored in log domain. +type Posterior = Log Double diff --git a/src/Mcmc/Prior.hs b/src/Mcmc/Prior.hs --- a/src/Mcmc/Prior.hs +++ b/src/Mcmc/Prior.hs @@ -2,8 +2,8 @@ -- | -- Module : Prior --- Description : Convenience functions for computing priors --- Copyright : (c) Dominik Schrempf, 2020 +-- Description : Types and convenience functions for computing priors +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -12,7 +12,11 @@ -- -- Creation date: Thu Jul 23 13:26:14 2020. module Mcmc.Prior - ( -- * Improper priors + ( Prior, + PriorFunction, + noPrior, + + -- * Improper priors largerThan, positive, lowerThan, @@ -21,6 +25,10 @@ -- * Continuous priors exponential, gamma, + gammaMeanVariance, + gammaMeanOne, + gammaShapeScaleToMeanVariance, + gammaMeanVarianceToShapeScale, normal, uniform, @@ -34,6 +42,7 @@ import Control.Monad import Data.Maybe (fromMaybe) +import Mcmc.Statistics.Types import Numeric.Log import qualified Statistics.Distribution as S import qualified Statistics.Distribution.Exponential as S @@ -41,80 +50,98 @@ import qualified Statistics.Distribution.Normal as S import qualified Statistics.Distribution.Poisson as S +-- | Prior values are stored in log domain. +type Prior = Log Double + +-- | Prior function. +type PriorFunction a = a -> Prior + +-- | Flat prior function. Useful for testing and debugging. +noPrior :: PriorFunction a +noPrior = const 1.0 + -- | Improper uniform prior; strictly larger than a given value. -largerThan :: Double -> Double -> Log Double +largerThan :: LowerBoundary -> PriorFunction Double largerThan a x | x <= a = 0 | otherwise = 1 -- | Improper uniform prior; strictly larger than zero. -positive :: Double -> Log Double +positive :: PriorFunction Double positive = largerThan 0 -- | Improper uniform prior; strictly lower than a given value. -lowerThan :: Double -> Double -> Log Double +lowerThan :: UpperBoundary -> PriorFunction Double lowerThan b x | x >= b = 0 | otherwise = 1 -- | Improper uniform prior; strictly lower than zero. -negative :: Double -> Log Double +negative :: PriorFunction Double negative = lowerThan 0 -- | Exponential distributed prior. -exponential :: - -- | Rate. - Double -> - Double -> - Log Double -exponential l x = Exp $ S.logDensity d x +exponential :: Rate -> PriorFunction Double +exponential l = Exp . S.logDensity d where d = S.exponential l -- | Gamma distributed prior. -gamma :: - -- | Shape. - Double -> - -- | Scale. - Double -> - Double -> - Log Double -gamma k t x = Exp $ S.logDensity d x +gamma :: Shape -> Scale -> PriorFunction Double +gamma k t = Exp . S.logDensity d where d = S.gammaDistr k t +-- | See 'gamma' but parametrized using mean and variance. +gammaMeanVariance :: Mean -> Variance -> PriorFunction Double +gammaMeanVariance m v = Exp . S.logDensity d + where (k, th) = gammaMeanVarianceToShapeScale m v + d = S.gammaDistr k th + +-- | Gamma disstributed prior with given shape and mean 1.0. +gammaMeanOne :: Shape -> PriorFunction Double +gammaMeanOne k = Exp . S.logDensity d + where + d = S.gammaDistr k (recip k) + +-- The mean and variance of the gamma distribution are +-- +-- m = k*t +-- +-- v = k*t*t +-- +-- Hence, the shape and scale are +-- +-- k = m^2/v +-- +-- t = v/m + +-- | Calculate mean and variance of the gamma distribution given the shape and +-- the scale. +gammaShapeScaleToMeanVariance :: Shape -> Scale -> (Mean, Variance) +gammaShapeScaleToMeanVariance k t = let m = k * t in (m, m * t) + +-- | Calculate shape and scale of the gamma distribution given the mean and +-- the variance. +gammaMeanVarianceToShapeScale :: Mean -> Variance -> (Shape, Scale) +gammaMeanVarianceToShapeScale m v = (m * m / v, v / m) + -- | Normal distributed prior. -normal :: - -- | Mean. - Double -> - -- | Standard deviation. - Double -> - Double -> - Log Double -normal m s x = Exp $ S.logDensity d x +normal :: Mean -> StandardDeviation -> PriorFunction Double +normal m s = Exp . S.logDensity d where d = S.normalDistr m s -- | Uniform prior on [a, b]. -uniform :: - -- | Lower bound a. - Double -> - -- | Upper bound b. - Double -> - Double -> - Log Double +uniform :: LowerBoundary -> UpperBoundary -> PriorFunction Double uniform a b x | x <= a = 0 | x >= b = 0 | otherwise = Exp 0 -- | Poisson distributed prior. -poisson :: - -- | Rate. - Double -> - Int -> - Log Double -poisson l x = Exp $ S.logProbability d x +poisson :: Rate -> PriorFunction Int +poisson l = Exp . S.logProbability d where d = S.poisson l diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs --- a/src/Mcmc/Proposal.hs +++ b/src/Mcmc/Proposal.hs @@ -6,7 +6,7 @@ -- | -- Module : Mcmc.Proposal -- Description : Proposals are instruction to move around the state space --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -15,34 +15,39 @@ -- -- Creation date: Wed May 20 13:42:53 2020. module Mcmc.Proposal - ( -- * Proposal + ( -- * Proposals and types PName (..), PDescription (..), - PWeight (..), + PWeight (fromPWeight), + pWeight, PDimension (..), Proposal (..), + KernelRatio, + Jacobian, + JacobianFunction, (@~), + liftProposal, + liftProposalWith, ProposalSimple, - Tuner (tParam, tFunc), + Tuner (..), Tune (..), createProposal, + TuningParameter, + tuningParameterMin, + tuningParameterMax, tune, getOptimalRate, - proposalHeader, - proposalHLine, - summarizeProposal, - -- * Cycle + -- * Cycles Order (..), Cycle (ccProposals), cycleFromList, setOrder, - orderProposals, + prepareProposals, tuneCycle, autoTuneCycle, - summarizeCycle, - -- * Acceptance + -- * Acceptance rates Acceptance (fromAcceptance), emptyA, pushA, @@ -50,6 +55,12 @@ transformKeysA, acceptanceRate, acceptanceRates, + + -- * Output + proposalHeader, + proposalHLine, + summarizeProposal, + summarizeCycle, ) where @@ -79,11 +90,16 @@ newtype PDescription = PDescription {fromPDescription :: String} deriving (Show, Eq, Ord) --- | The weight determines how often a 'Proposal' is executed per iteration of --- the Markov chain. +-- | The positive weight determines how often a 'Proposal' is executed per +-- iteration of the Markov chain. newtype PWeight = PWeight {fromPWeight :: Int} deriving (Show, Eq, Ord) +-- | Check if the weight is positive. +pWeight :: Int -> PWeight +pWeight n | n <= 0 = error "pWeight: Proposal weight is zero or negative." + | otherwise = PWeight n + -- | Proposal dimension. -- -- The number of affected, independent parameters. @@ -127,48 +143,90 @@ -- -- A 'Proposal' may be tuneable in that it contains information about how to enlarge -- or shrink the step size to tune the acceptance rate. +-- +-- Predefined proposals are provided. To create custom proposals, one may use +-- the convenience function 'createProposal'. data Proposal a = Proposal { -- | Name of the affected variable. - pName :: PName, + prName :: PName, -- | Description of the proposal type and parameters. - pDescription :: PDescription, + prDescription :: PDescription, -- | Dimension of the proposal. The dimension is used to calculate the -- optimal acceptance rate, and does not have to be exact. - pDimension :: PDimension, + prDimension :: PDimension, -- | The weight determines how often a 'Proposal' is executed per iteration of -- the Markov chain. - pWeight :: PWeight, + prWeight :: PWeight, -- | Simple proposal without name, weight, and tuning information. - pSimple :: ProposalSimple a, + prSimple :: ProposalSimple a, -- | Tuning is disabled if set to 'Nothing'. - pTuner :: Maybe (Tuner a) + prTuner :: Maybe (Tuner a) } instance Eq (Proposal a) where - m == n = pName m == pName n && pDescription m == pDescription n + m == n = prName m == prName n && prDescription m == prDescription n instance Ord (Proposal a) where - compare = compare `on` (\p -> (pDescription p, pName p, pWeight p)) + compare = compare `on` (\p -> (prDescription p, prName p, prWeight p)) --- | Convert a proposal from one data type to another using a lens. +-- | Ratio of the proposal kernels. -- +-- Part of the MHG acceptance ratio. +-- +-- See also 'Jacobian'. +-- +-- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However, +-- it is more declarative to have them separate. It is a constant reminder: Is +-- the Jacobian modifier different from 1.0? +type KernelRatio = Log Double + +-- | Absolute value of the determinant of the Jacobian matrix. +-- +-- Part of the MHG acceptance ratio. +-- +-- See also 'Jacobian'. +type Jacobian = Log Double + +-- | Function calculating the 'Jacobian' of a proposal. +type JacobianFunction a = a -> Jacobian + +-- | Lift a proposal from one data type to another. +-- +-- Assume the Jacobian is 1.0 (see also 'liftProposal' and 'liftProposalWith'). +-- -- For example: -- -- @ -- scaleFirstEntryOfTuple = _1 @~ scale -- @ +infixl 7 @~ (@~) :: Lens' b a -> Proposal a -> Proposal b -(@~) l (Proposal n r d w s t) = Proposal n r d w (convertProposalSimple l s) (convertTuner l <$> t) +(@~) = liftProposal +-- | Lift a proposal from one data type to another. +-- +-- Assume the Jacobian is 1.0 (see also '(@~)' and 'liftProposalWith'). +liftProposal :: Lens' b a -> Proposal a -> Proposal b +liftProposal = liftProposalWith (const 1.0) + +-- | Lift a proposal from one data type to another. +-- +-- A function to calculate the Jacobian has to be provided (see also +-- 'liftProposal'). +-- +-- For further reference, please see the [example +-- @Pair@](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs). +liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b +liftProposalWith jf l (Proposal n r d w s t) = Proposal n r d w (convertProposalSimple jf l s) (convertTuner jf l <$> t) + -- | Simple proposal without tuning information. -- -- Instruction about randomly moving from the current state to a new state, -- given some source of randomness. -- -- In order to calculate the Metropolis-Hastings-Green ratio, we need to know --- the ratio of the backward to forward kernels (i.e., the probability masses or --- probability densities) and the absolute value of the determinant of the --- Jacobian matrix. +-- the ratio of the backward to forward kernels (the 'KernelRatio' or the +-- probability masses or probability densities) and the 'Jacobian'. -- -- For unbiased proposals, these values are 1.0 such that -- @@ -179,38 +237,45 @@ -- For biased proposals, the kernel ratio is qYX / qXY, where qXY is the -- probability density to move from X to Y, and the absolute value of the -- determinant of the Jacobian matrix differs from 1.0. -type ProposalSimple a = a -> GenIO -> IO (a, Log Double, Log Double) +type ProposalSimple a = a -> GenIO -> IO (a, KernelRatio, Jacobian) -convertProposalSimple :: Lens' b a -> ProposalSimple a -> ProposalSimple b -convertProposalSimple l s = s' +convertProposalSimple :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b +convertProposalSimple jf l s = s' where - s' v g = do - (x', r, j) <- s (v ^. l) g - return (set l x' v, r, j) + s' y g = do + (x', r, j) <- s (y ^. l) g + let y' = set l x' y + jxy = jf y + jyx = jf y' + j' = j * jyx / jxy + return (y', r, j') -- | Tune the acceptance rate of a 'Proposal'; see 'tune', or 'autoTuneCycle'. data Tuner a = Tuner - { tParam :: Double, - tFunc :: Double -> ProposalSimple a + { tParam :: TuningParameter, + tFunc :: TuningParameter -> ProposalSimple a } -convertTuner :: Lens' b a -> Tuner a -> Tuner b -convertTuner l (Tuner p f) = Tuner p f' +convertTuner :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b +convertTuner jf l (Tuner p f) = Tuner p f' where - f' x = convertProposalSimple l $ f x + f' x = convertProposalSimple jf l $ f x -- | Tune the proposal? data Tune = Tune | NoTune deriving (Show, Eq) +-- | Tuning parameter. +type TuningParameter = Double + -- | Create a tuneable proposal. createProposal :: -- | Description of the proposal type and parameters. PDescription -> -- | Function creating a simple proposal for a given tuning parameter. The - -- larger the tuning parameter, the larger the proposal (and the lower the - -- expected acceptance rate), and vice versa. - (Double -> ProposalSimple a) -> + -- larger the tuning parameter, the larger the proposal and the lower the + -- expected acceptance rate; and vice versa. + (TuningParameter -> ProposalSimple a) -> -- | Dimension. PDimension -> -- | Name. @@ -223,19 +288,34 @@ createProposal r f d n w Tune = Proposal n r d w (f 1.0) (Just $ Tuner 1.0 f) createProposal r f d n w NoTune = Proposal n r d w (f 1.0) Nothing --- Minimal tuning parameter; subject to change. -tuningParamMin :: Double -tuningParamMin = 1e-12 +-- | Minimal tuning parameter; @1e-12@, subject to change. +-- +-- >>> tuningParameterMin +-- 1e-5 +tuningParameterMin :: TuningParameter +tuningParameterMin = 1e-5 --- | Tune a 'Proposal'. Return 'Nothing' if 'Proposal' is not tuneable. The size --- of the proposal is proportional to the tuning parameter. Negative tuning --- parameters are not allowed. -tune :: (Double -> Double) -> Proposal a -> Maybe (Proposal a) +-- | Maximal tuning parameter; @1e12@, subject to change. +-- >>> tuningParameterMax +-- 1e3 +tuningParameterMax :: TuningParameter +tuningParameterMax = 1e3 + +-- | Tune a 'Proposal'. +-- +-- The size of the proposal is proportional to the tuning parameter which has a +-- positive lower bound of 'tuningParameterMin'. +-- +-- The tuning function maps the current tuning parameter to a new one. +-- +-- Return 'Nothing' if 'Proposal' is not tuneable. +tune :: (TuningParameter -> TuningParameter) -> Proposal a -> Maybe (Proposal a) tune f m = do - (Tuner t g) <- pTuner m - -- Ensure that the tuning parameter is strictly positive. - let t' = max tuningParamMin (f t) - return $ m {pSimple = g t', pTuner = Just $ Tuner t' g} + (Tuner t g) <- prTuner m + -- Ensure that the tuning parameter is strictly positive and well bounded. + let t' = max tuningParameterMin (f t) + t'' = min tuningParameterMax t' + return $ m {prSimple = g t'', prTuner = Just $ Tuner t'' g} -- | See 'PDimension'. getOptimalRate :: PDimension -> Double @@ -329,8 +409,8 @@ setOrder o c = c {ccOrder = o} -- | Replicate 'Proposal's according to their weights and possibly shuffle them. -orderProposals :: Cycle a -> GenIO -> IO [Proposal a] -orderProposals (Cycle xs o) g = case o of +prepareProposals :: Cycle a -> GenIO -> IO [Proposal a] +prepareProposals (Cycle xs o) g = case o of RandomO -> shuffle ps g SequentialO -> return ps RandomReversibleO -> do @@ -338,7 +418,7 @@ return $ psR ++ reverse psR SequentialReversibleO -> return $ ps ++ reverse ps where - !ps = concat [replicate (fromPWeight $ pWeight p) p | p <- xs] + !ps = concat [replicate (fromPWeight $ prWeight p) p | p <- xs] -- The number of proposals depends on the order. getNProposalsPerCycle :: Cycle a -> Int @@ -348,10 +428,10 @@ RandomReversibleO -> 2 * once SequentialReversibleO -> 2 * once where - once = sum $ map (fromPWeight . pWeight) xs + once = sum $ map (fromPWeight . prWeight) xs -- | Tune 'Proposal's in the 'Cycle'. See 'tune'. -tuneCycle :: M.Map (Proposal a) (Double -> Double) -> Cycle a -> Cycle a +tuneCycle :: M.Map (Proposal a) (TuningParameter -> TuningParameter) -> Cycle a -> Cycle a tuneCycle m c = if sort (M.keys m) == sort ps then c {ccProposals = map tuneF ps} @@ -368,9 +448,11 @@ autoTuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a autoTuneCycle a = tuneCycle (M.mapWithKey tuningF $ acceptanceRates a) where - tuningF proposal currentRate currentTuningParam = - let optimalRate = getOptimalRate (pDimension proposal) - in exp (currentRate - optimalRate) * currentTuningParam + tuningF proposal mCurrentRate currentTuningParam = case mCurrentRate of + Nothing -> currentTuningParam + Just currentRate -> + let optimalRate = getOptimalRate (prDimension proposal) + in exp (2 * (currentRate - optimalRate)) * currentTuningParam renderRow :: BL.ByteString -> @@ -395,6 +477,67 @@ tp = alignRight 20 tuneParam mt = alignRight 30 manualAdjustment +-- | For each key @k@, store the number of accepted and rejected proposals. +newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)} + deriving (Eq, Read, Show) + +instance ToJSONKey k => ToJSON (Acceptance k) where + toJSON (Acceptance m) = toJSON m + toEncoding (Acceptance m) = toEncoding m + +instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where + parseJSON v = Acceptance <$> parseJSON v + +-- | In the beginning there was the Word. +-- +-- Initialize an empty storage of accepted/rejected values. +emptyA :: Ord k => [k] -> Acceptance k +emptyA ks = Acceptance $ M.fromList [(k, (0, 0)) | k <- ks] + +-- | For key @k@, prepend an accepted (True) or rejected (False) proposal. +pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k +pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance +pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance +{-# INLINEABLE pushA #-} + +-- | Reset acceptance storage. +resetA :: Ord k => Acceptance k -> Acceptance k +resetA = emptyA . M.keys . fromAcceptance + +transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> M.Map k1 v -> M.Map k2 v +transformKeys ks1 ks2 m = foldl' insrt M.empty $ zip ks1 ks2 + where + insrt m' (k1, k2) = M.insert k2 (m M.! k1) m' + +-- | Transform keys using the given lists. Keys not provided will not be present +-- in the new 'Acceptance' variable. +transformKeysA :: (Ord k1, Ord k2) => [k1] -> [k2] -> Acceptance k1 -> Acceptance k2 +transformKeysA ks1 ks2 = Acceptance . transformKeys ks1 ks2 . fromAcceptance + +-- | Acceptance counts and rate for a specific proposal. +-- +-- Return 'Nothing' if no proposals have been accepted or rejected (division by +-- zero). +acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, Double) +acceptanceRate k a = case fromAcceptance a M.!? k of + Just (0, 0) -> Nothing + Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs)) + Nothing -> error "acceptanceRate: Key not found in map." + +-- | Acceptance rates for all proposals. +-- +-- Set rate to 'Nothing' if no proposals have been accepted or rejected +-- (division by zero). +acceptanceRates :: Acceptance k -> M.Map k (Maybe Double) +acceptanceRates = + M.map + ( \(as, rs) -> + if as + rs == 0 + then Nothing + else Just $ fromIntegral as / fromIntegral (as + rs) + ) + . fromAcceptance + -- | Header of proposal summaries. proposalHeader :: BL.ByteString proposalHeader = @@ -418,12 +561,11 @@ PName -> PDescription -> PWeight -> - -- Tuning parameter. - Maybe Double -> + Maybe TuningParameter -> PDimension -> Maybe (Int, Int, Double) -> BL.ByteString -summarizeProposal name description weight tuningParam dimension r = +summarizeProposal name description weight tuningParameter dimension ar = renderRow (BL.pack $ fromPName name) (BL.pack $ fromPDescription description) @@ -436,16 +578,27 @@ manualAdjustmentStr where weightStr = BB.toLazyByteString $ BB.intDec $ fromPWeight weight - nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) r - nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) r - acceptRate = BL.fromStrict $ maybe "" (BC.toFixed 2 . (^. _3)) r + nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) ar + nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) ar + acceptRate = BL.fromStrict $ maybe "" (BC.toFixed 2 . (^. _3)) ar optimalRate = BL.fromStrict $ BC.toFixed 2 $ getOptimalRate dimension - tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) tuningParam - check v - | v < rateMin = "rate too low" - | v > rateMax = "rate too high" - | otherwise = "" - manualAdjustmentStr = BL.fromStrict $ maybe "" (check . (^. _3)) r + tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) tuningParameter + checkRate rate + | rate < rateMin = Just "rate too low" + | rate > rateMax = Just "rate too high" + | otherwise = Nothing + checkTuningParam tp + | tp <= (1.1 * tuningParameterMin) = Just "tuning parameter too low" + | tp >= (0.9 * tuningParameterMax) = Just "tuning parameter too high" + | otherwise = Nothing + tps = checkTuningParam =<< tuningParameter + ars = (checkRate . (^. _3)) =<< ar + manualAdjustmentStr = + let + in case (ars, tps) of + (Nothing, Nothing) -> "" + (Just s, _) -> s + (_, Just s) -> s -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates. summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString @@ -458,11 +611,11 @@ proposalHLine ] ++ [ summarizeProposal - (pName p) - (pDescription p) - (pWeight p) - (tParam <$> pTuner p) - (pDimension p) + (prName p) + (prDescription p) + (prWeight p) + (tParam <$> prTuner p) + (prDimension p) (ar p) | p <- ps ] @@ -475,51 +628,3 @@ 1 -> nProposalsStr <> " proposal is performed per iteration." _ -> nProposalsStr <> " proposals are performed per iterations." ar m = acceptanceRate m a - --- | For each key @k@, store the number of accepted and rejected proposals. -newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)} - deriving (Eq, Read, Show) - -instance ToJSONKey k => ToJSON (Acceptance k) where - toJSON (Acceptance m) = toJSON m - toEncoding (Acceptance m) = toEncoding m - -instance (Ord k, FromJSONKey k) => FromJSON (Acceptance k) where - parseJSON v = Acceptance <$> parseJSON v - --- | In the beginning there was the Word. --- --- Initialize an empty storage of accepted/rejected values. -emptyA :: Ord k => [k] -> Acceptance k -emptyA ks = Acceptance $ M.fromList [(k, (0, 0)) | k <- ks] - --- | For key @k@, prepend an accepted (True) or rejected (False) proposal. -pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k -pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance -pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance -{-# INLINEABLE pushA #-} - --- | Reset acceptance storage. -resetA :: Ord k => Acceptance k -> Acceptance k -resetA = emptyA . M.keys . fromAcceptance - -transformKeys :: (Ord k1, Ord k2) => [k1] -> [k2] -> M.Map k1 v -> M.Map k2 v -transformKeys ks1 ks2 m = foldl' insrt M.empty $ zip ks1 ks2 - where - insrt m' (k1, k2) = M.insert k2 (m M.! k1) m' - --- | Transform keys using the given lists. Keys not provided will not be present --- in the new 'Acceptance' variable. -transformKeysA :: (Ord k1, Ord k2) => [k1] -> [k2] -> Acceptance k1 -> Acceptance k2 -transformKeysA ks1 ks2 = Acceptance . transformKeys ks1 ks2 . fromAcceptance - --- | Acceptance counts and rate for a specific proposal. -acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, Double) -acceptanceRate k a = case fromAcceptance a M.!? k of - Just (0, 0) -> Nothing - Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs)) - Nothing -> error "acceptanceRate: Key not found in map." - --- | Acceptance rates for all proposals. -acceptanceRates :: Acceptance k -> M.Map k Double -acceptanceRates = M.map (\(as, rs) -> fromIntegral as / fromIntegral (as + rs)) . fromAcceptance diff --git a/src/Mcmc/Proposal/Bactrian.hs b/src/Mcmc/Proposal/Bactrian.hs --- a/src/Mcmc/Proposal/Bactrian.hs +++ b/src/Mcmc/Proposal/Bactrian.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Proposal.Bactrian -- Description : Bactrian proposals --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -12,21 +12,26 @@ -- -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. module Mcmc.Proposal.Bactrian - ( slideBactrian, + ( SpikeParameter, + slideBactrian, scaleBactrian, ) where import Mcmc.Proposal +import Mcmc.Statistics.Types import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Normal import System.Random.MWC import System.Random.MWC.Distributions +-- | Type synonym indicating the spike parameter. +type SpikeParameter = Double + genBactrian :: - Double -> - Double -> + SpikeParameter -> + StandardDeviation -> GenIO -> IO Double genBactrian m s g = do @@ -37,7 +42,7 @@ b <- bernoulli 0.5 g return $ if b then x else - x -logDensityBactrian :: Double -> Double -> Double -> Log Double +logDensityBactrian :: SpikeParameter -> StandardDeviation -> Double -> Log Double logDensityBactrian m s x = Exp $ log $ kernel1 + kernel2 where mn = m * s @@ -48,20 +53,18 @@ kernel2 = density dist2 x bactrianAdditive :: - Double -> - Double -> - Double -> - GenIO -> - IO (Double, Log Double, Log Double) + SpikeParameter -> + StandardDeviation -> + ProposalSimple Double bactrianAdditive m s x g = do dx <- genBactrian m s g return (x + dx, 1.0, 1.0) -- bactrianSimple lens spike stdDev tune forwardOp backwardOp bactrianAdditiveSimple :: - Double -> - Double -> - Double -> + SpikeParameter -> + StandardDeviation -> + TuningParameter -> ProposalSimple Double bactrianAdditiveSimple m s t | m < 0 = error "bactrianAdditiveSimple: Spike parameter negative." @@ -81,15 +84,10 @@ -- -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. slideBactrian :: - -- | Spike parameter \(m\). - Double -> - -- | Standard deviation \(s\). - Double -> - -- | Name. + SpikeParameter -> + StandardDeviation -> PName -> - -- | Weight. PWeight -> - -- | Enable tuning. Tune -> Proposal Double slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) (PDimension 1) @@ -106,11 +104,9 @@ fInv dx = recip (1 - dx) - 1 bactrianMult :: - Double -> - Double -> - Double -> - GenIO -> - IO (Double, Log Double, Log Double) + SpikeParameter -> + StandardDeviation -> + ProposalSimple Double bactrianMult m s x g = do du <- genBactrian m s g let qXY = logDensityBactrian m s du @@ -119,7 +115,7 @@ jac = Exp $ log $ recip u return (x * u, qYX / qXY, jac) -bactrianMultSimple :: Double -> Double -> Double -> ProposalSimple Double +bactrianMultSimple :: SpikeParameter -> StandardDeviation -> TuningParameter -> ProposalSimple Double bactrianMultSimple m s t | m < 0 = error "bactrianMultSimple: Spike parameter negative." | m >= 1 = error "bactrianMultSimple: Spike parameter 1.0 or larger." @@ -129,15 +125,10 @@ -- | Multiplicative proposal with kernel similar to the silhouette of a Bactrian -- camel. See 'slideBactrian'. scaleBactrian :: - -- | Spike parameter. - Double -> - -- | Standard deviation. - Double -> - -- | Name. + SpikeParameter -> + StandardDeviation -> PName -> - -- | Weight. PWeight -> - -- | Enable tuning. Tune -> Proposal Double scaleBactrian m s = createProposal description (bactrianMultSimple m s) (PDimension 1) diff --git a/src/Mcmc/Proposal/Generic.hs b/src/Mcmc/Proposal/Generic.hs --- a/src/Mcmc/Proposal/Generic.hs +++ b/src/Mcmc/Proposal/Generic.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.Proposal.Generic -- Description : Generic interface for creating proposals --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -51,7 +51,7 @@ -- -- Conversion to log domain is necessary, because some determinants of -- Jacobians are very small (or large). - Maybe (a -> Double -> Log Double) -> + Maybe (a -> Double -> Jacobian) -> ProposalSimple a genericContinuous d f mInv mJac x g = do u <- genContVar d g diff --git a/src/Mcmc/Proposal/Scale.hs b/src/Mcmc/Proposal/Scale.hs --- a/src/Mcmc/Proposal/Scale.hs +++ b/src/Mcmc/Proposal/Scale.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Proposal.Scale -- Description : Multiplicative proposals --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -20,12 +20,13 @@ import Mcmc.Proposal import Mcmc.Proposal.Generic +import Mcmc.Statistics.Types import Numeric.Log import Statistics.Distribution.Gamma -- The actual proposal with tuning parameter. The tuning parameter does not -- change the mean. -scaleSimple :: Double -> Double -> Double -> ProposalSimple Double +scaleSimple :: Shape -> Scale -> TuningParameter -> ProposalSimple Double scaleSimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) @@ -37,10 +38,8 @@ -- | Multiplicative proposal with gamma distributed kernel. scale :: - -- | Shape. - Double -> - -- | Scale. - Double -> + Shape -> + Scale -> PName -> PWeight -> Tune -> @@ -54,8 +53,7 @@ -- The scale of the gamma distribution is set to (shape)^{-1}, so that the mean -- of the gamma distribution is 1.0. scaleUnbiased :: - -- | Shape. - Double -> + Shape -> PName -> PWeight -> Tune -> @@ -64,7 +62,7 @@ where description = PDescription $ "Scale unbiased; shape: " ++ show k -scaleContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double) +scaleContrarilySimple :: Shape -> Scale -> TuningParameter -> ProposalSimple (Double, Double) scaleContrarilySimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) @@ -80,10 +78,8 @@ -- The two values are scaled contrarily so that their product stays constant. -- Contrary proposals are useful when parameters are confounded. scaleContrarily :: - -- | Shape. - Double -> - -- | Scale. - Double -> + Shape -> + Scale -> PName -> PWeight -> Tune -> diff --git a/src/Mcmc/Proposal/Simplex.hs b/src/Mcmc/Proposal/Simplex.hs --- a/src/Mcmc/Proposal/Simplex.hs +++ b/src/Mcmc/Proposal/Simplex.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Proposal.Simplex -- Description : Proposals on simplices --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -27,6 +27,7 @@ import Data.Aeson.TH import qualified Data.Vector.Unboxed as V import Mcmc.Proposal +import Mcmc.Statistics.Types import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Beta @@ -77,11 +78,12 @@ -- | Create the uniform element of the K-dimensional simplex. -- -- Set all values to \(1/D\). -simplexUniform :: Int -> Simplex +simplexUniform :: Dimension -> Simplex +-- Should never fail. simplexUniform k = either error id $ simplexFromVector $ V.replicate k (1.0 / fromIntegral k) -- Tuning function is inverted (high alpha means small steps). -getTuningFunction :: Double -> (Double -> Double) +getTuningFunction :: TuningParameter -> (TuningParameter -> TuningParameter) getTuningFunction t = (/ t'') where -- Start with small steps. @@ -98,7 +100,7 @@ -- The values determining the proposal size have been set using an example -- analysis. They are good values for this analysis, but may fail for other -- analyses. -dirichletSimple :: Double -> ProposalSimple Simplex +dirichletSimple :: TuningParameter -> ProposalSimple Simplex dirichletSimple t (SimplexUnsafe xs) g = do -- If @t@ is high and above 1.0, the parameter vector will be low, and the -- variance will be high. If @t@ is low and below 1.0, the parameter vector @@ -135,7 +137,7 @@ -- The proposal dimension, which is the dimension of the simplex, is used to -- determine the optimal acceptance rate. -- --- For high dimensional simplices, this proposal may have low acceptance ratios. +-- For high dimensional simplices, this proposal may have low acceptance rates. -- In this case, please see the coordinate wise 'beta' proposal. dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple @@ -147,7 +149,7 @@ -- analyses. -- -- See also the 'dirichlet' proposal. -betaSimple :: Int -> Double -> ProposalSimple Simplex +betaSimple :: Dimension -> TuningParameter -> ProposalSimple Simplex betaSimple i t (SimplexUnsafe xs) g = do -- Shape parameters of beta distribution. Do not assume that the sum of the -- elements of 'xs' is 1.0, because then repeated proposals let the sum of the @@ -203,7 +205,7 @@ -- -- This proposal has been assigned a dimension of 2. See the discussion at -- 'PDimension'. -beta :: Int -> PName -> PWeight -> Tune -> Proposal Simplex +beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex beta i = createProposal description (betaSimple i) (PDimension 2) where description = PDescription $ "Beta; coordinate: " ++ show i diff --git a/src/Mcmc/Proposal/Slide.hs b/src/Mcmc/Proposal/Slide.hs --- a/src/Mcmc/Proposal/Slide.hs +++ b/src/Mcmc/Proposal/Slide.hs @@ -3,7 +3,7 @@ -- | -- Module : Mcmc.Proposal.Slide -- Description : Additive proposals --- Copyright : (c) Dominik Schrempf 2020 +-- Copyright : (c) Dominik Schrempf 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -21,25 +21,21 @@ import Mcmc.Proposal import Mcmc.Proposal.Generic +import Mcmc.Statistics.Types import Statistics.Distribution.Normal import Statistics.Distribution.Uniform -- The actual proposal with tuning parameter. -slideSimple :: Double -> Double -> Double -> ProposalSimple Double +slideSimple :: Mean -> StandardDeviation -> TuningParameter -> ProposalSimple Double slideSimple m s t = genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing -- | Additive proposal with normally distributed kernel. slide :: - -- | Mean. - Double -> - -- | Standard deviation. - Double -> - -- | Name. + Mean -> + StandardDeviation -> PName -> - -- | Weight. PWeight -> - -- | Enable tuning. Tune -> Proposal Double slide m s = createProposal description (slideSimple m s) (PDimension 1) @@ -47,7 +43,7 @@ description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s -- The actual proposal with tuning parameter. -slideSymmetricSimple :: Double -> Double -> ProposalSimple Double +slideSymmetricSimple :: StandardDeviation -> TuningParameter -> ProposalSimple Double slideSymmetricSimple s t = genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing @@ -55,13 +51,9 @@ -- proposal is very fast, because the Metropolis-Hastings-Green ratio does not -- include calculation of the forwards and backwards kernels. slideSymmetric :: - -- | Standard deviation. - Double -> - -- | Name. + StandardDeviation -> PName -> - -- | Weight. PWeight -> - -- | Enable tuning. Tune -> Proposal Double slideSymmetric s = createProposal description (slideSymmetricSimple s) (PDimension 1) @@ -69,7 +61,7 @@ description = PDescription $ "Slide symmetric; sd: " ++ show s -- The actual proposal with tuning parameter. -slideUniformSimple :: Double -> Double -> ProposalSimple Double +slideUniformSimple :: Size -> TuningParameter -> ProposalSimple Double slideUniformSimple d t = genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing @@ -77,13 +69,9 @@ -- proposal is very fast, because the Metropolis-Hastings-Green ratio does not -- include calculation of the forwards and backwards kernels. slideUniformSymmetric :: - -- | Delta. - Double -> - -- | Name. + Size -> PName -> - -- | Weight. PWeight -> - -- | Enable tuning. Tune -> Proposal Double slideUniformSymmetric d = createProposal description (slideUniformSimple d) (PDimension 1) @@ -93,7 +81,11 @@ contra :: (Double, Double) -> Double -> (Double, Double) contra (x, y) u = (x + u, y - u) -slideContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double) +slideContrarilySimple :: + Mean -> + StandardDeviation -> + TuningParameter -> + ProposalSimple (Double, Double) slideContrarilySimple m s t = genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing @@ -102,15 +94,10 @@ -- The two values are slid contrarily so that their sum stays constant. Contrary -- proposals are useful when parameters are confounded. slideContrarily :: - -- | Mean. - Double -> - -- | Standard deviation. - Double -> - -- | Name. + Mean -> + StandardDeviation -> PName -> - -- | Weight. PWeight -> - -- | Enable tuning. Tune -> Proposal (Double, Double) slideContrarily m s = createProposal description (slideContrarilySimple m s) (PDimension 2) diff --git a/src/Mcmc/Settings.hs b/src/Mcmc/Settings.hs --- a/src/Mcmc/Settings.hs +++ b/src/Mcmc/Settings.hs @@ -1,10 +1,11 @@ +{-# LANGUAGE DerivingVia #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE TemplateHaskell #-} -- | -- Module : Mcmc.Settings -- Description : Settings of Markov chain Monte Carlo samplers --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -15,13 +16,17 @@ module Mcmc.Settings ( -- * Data types AnalysisName (..), + HasAnalysisName (..), BurnInSpecification (..), burnInIterations, Iterations (..), + TraceLength (..), ExecutionMode (..), + HasExecutionMode (..), openWithExecutionMode, ParallelizationMode (..), SaveMode (..), + LogMode (..), Verbosity (..), -- * Settings @@ -35,15 +40,21 @@ import Data.Aeson import Data.Aeson.TH import qualified Data.ByteString.Lazy.Char8 as BL +import Mcmc.Logger import System.Directory import System.IO -- | Analysis name of the MCMC sampler. newtype AnalysisName = AnalysisName {fromAnalysisName :: String} deriving (Eq, Read, Show) + deriving (Monoid, Semigroup) via String $(deriveJSON defaultOptions ''AnalysisName) +-- | Types with analysis names. +class HasAnalysisName s where + getAnalysisName :: s -> AnalysisName + -- | Burn in specification. data BurnInSpecification = -- | No burn in. @@ -53,15 +64,32 @@ | -- | Burn in for a given number of iterations. Enable auto tuning with a -- given period. BurnInWithAutoTuning Int Int + | -- | Burn in with the given list of auto tuning periods. + -- + -- For example, @BurnInWithCustomAutoTuning [100,200]@ performs 300 + -- iterations with two auto tuning steps. One after 100 iterations, the + -- second one after 200 more iterations. + -- + -- Usually it is useful to auto tune more frequently in the beginning of the + -- MCMC run. + BurnInWithCustomAutoTuning [Int] deriving (Eq, Read, Show) $(deriveJSON defaultOptions ''BurnInSpecification) +-- Check if the burn in specification is valid. +burnInValid :: BurnInSpecification -> Bool +burnInValid NoBurnIn = True +burnInValid (BurnInWithoutAutoTuning n) = n > 0 +burnInValid (BurnInWithAutoTuning n t) = n > 0 && t > 0 +burnInValid (BurnInWithCustomAutoTuning xs) = not (null xs) && all (> 0) xs + -- | Get the number of burn in iterations. burnInIterations :: BurnInSpecification -> Int burnInIterations NoBurnIn = 0 burnInIterations (BurnInWithoutAutoTuning n) = n burnInIterations (BurnInWithAutoTuning n _) = n +burnInIterations (BurnInWithCustomAutoTuning xs) = sum xs -- | Number of normal iterations after burn in. -- @@ -71,6 +99,17 @@ $(deriveJSON defaultOptions ''Iterations) +-- | The length of the stored "Mcmc.Chain.Trace". +-- +-- Be careful, this setting determines the memory requirement of the MCMC chain. +data TraceLength + = -- | Automatically determine the length of the trace. The value is + -- determined by the 'Mcmc.Monitor.MonitorBatch' with largest batch size. + TraceAuto + | -- | Store a given minimum number of iterations of the chain. Store more + -- iterations if required (see 'TraceAuto'). + TraceMinimum Int + -- | Execution mode. data ExecutionMode = -- | Perform new run. @@ -89,6 +128,10 @@ $(deriveJSON defaultOptions ''ExecutionMode) +-- | Types with execution modes. +class HasExecutionMode s where + getExecutionMode :: s -> ExecutionMode + -- | Open a file honoring the execution mode. -- -- Call 'error' if execution mode is @@ -121,17 +164,16 @@ -- | Parallelization mode. -- -- Parallel execution of the chains is only beneficial when the algorithm allows --- for parallelization, and if computation of the next iteration takes a long +-- for parallelization, and if computation of the next iteration takes some -- time. If the calculation of the next state is fast, sequential execution is --- usually beneficial, even for algorithms involving parallel chains. If the --- calculation of the next state is slow, parallel execution may be beneficial. +-- usually beneficial, even for algorithms involving parallel chains. -- --- - The "Mcmc.Algorithm.Metropolis" algorithm is inherently sequential. +-- - The "Mcmc.Algorithm.MHG" algorithm is inherently sequential. -- -- - The "Mcmc.Algorithm.MC3" algorithm works well with parallelization. -- -- Of course, also the prior or likelihood functions can be computed in --- parallel. However, this library is not aware of how these functions are +-- parallel. However, this library is unaware about how these functions are -- computed. data ParallelizationMode = Sequential @@ -140,20 +182,18 @@ $(deriveJSON defaultOptions ''ParallelizationMode) --- | Should the MCMC run be saved at the end of the run? -data SaveMode = NoSave | Save +-- | Define information stored on disk. +data SaveMode + = -- | Do not save the MCMC analysis. The analysis can not be continued. + NoSave + | -- | Save the MCMC analysis so that it can be continued. This can be slow, + -- if the trace is long, or if the states are large objects. See + -- 'TraceLength'. + Save deriving (Eq, Read, Show) $(deriveJSON defaultOptions ''SaveMode) --- $(deriveJSON defaultOptions ''CleaningMode) - --- | Not much to say here. -data Verbosity = Quiet | Warn | Info | Debug - deriving (Eq, Ord, Read, Show) - -$(deriveJSON defaultOptions ''Verbosity) - -- | Settings of an MCMC sampler. data Settings = Settings { sAnalysisName :: AnalysisName, @@ -162,14 +202,27 @@ sExecutionMode :: ExecutionMode, sParallelizationMode :: ParallelizationMode, sSaveMode :: SaveMode, + sLogMode :: LogMode, sVerbosity :: Verbosity } deriving (Eq, Show) +instance HasAnalysisName Settings where + getAnalysisName = sAnalysisName + +instance HasExecutionMode Settings where + getExecutionMode = sExecutionMode + +instance HasLogMode Settings where + getLogMode = sLogMode + +instance HasVerbosity Settings where + getVerbosity = sVerbosity + $(deriveJSON defaultOptions ''Settings) settingsFn :: String -> FilePath -settingsFn n = n ++ ".settings" +settingsFn n = n ++ ".mcmc.settings" -- | Save settings to a file determined by the analysis name. settingsSave :: Settings -> IO () @@ -217,10 +270,10 @@ -- | Current iteration. Int -> IO () -settingsCheck s@(Settings nm bi i em _ _ _) iCurrent +settingsCheck s@(Settings nm bi i em _ _ _ _) iCurrent | null (fromAnalysisName nm) = serr "Analysis name is the empty string." | burnInIterations bi < 0 = serr "Number of burn in iterations is negative." - | not $ burnInAutoTuningPeriodValid bi = serr "Auto tuning period is zero or negative." + | not $ burnInValid bi = serr $ "Burn in setting invalid: " <> show bi | fromIterations i < 0 = serr "Number of iterations is negative." | burnInIterations bi + fromIterations i - iCurrent < 0 = serr "Current iteration is larger than the total number of iterations." @@ -231,6 +284,3 @@ | otherwise = return () where serr = settingsError s iCurrent - burnInAutoTuningPeriodValid :: BurnInSpecification -> Bool - burnInAutoTuningPeriodValid (BurnInWithAutoTuning _ t) = t > 0 - burnInAutoTuningPeriodValid _ = True diff --git a/src/Mcmc/Statistics/Types.hs b/src/Mcmc/Statistics/Types.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Statistics/Types.hs @@ -0,0 +1,56 @@ +-- | +-- Module : Mcmc.Statistics.Types +-- Description : Types indicating properties of distributions +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Tue Feb 23 14:51:06 2021. +module Mcmc.Statistics.Types + ( Mean, + StandardDeviation, + Variance, + Shape, + Scale, + Rate, + Dimension, + Size, + LowerBoundary, + UpperBoundary, + ) +where + +-- | Mean of a distribution. +type Mean = Double + +-- | Standard deviation of a distribution. +type StandardDeviation = Double + +-- | Variance of a distribution. +type Variance = Double + +-- | Shape of a distribution. +type Shape = Double + +-- | Scale of a distribution. +type Scale = Double + +-- | Rate of a distribution. +type Rate = Double + +-- | Dimension of a distribution. +type Dimension = Int + +-- | Size of a distribution. +-- +-- For example, the size of the interval of the uniform distribution. +type Size = Double + +-- | Lower boundary of a distribution. +type LowerBoundary = Double + +-- | Upper boundary of a distribution. +type UpperBoundary = Double diff --git a/test/Mcmc/ProposalSpec.hs b/test/Mcmc/ProposalSpec.hs --- a/test/Mcmc/ProposalSpec.hs +++ b/test/Mcmc/ProposalSpec.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.ProposalSpec -- Description : Unit tests for Mcmc.Proposal --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -20,24 +20,24 @@ import Test.Hspec p1 :: Proposal Double -p1 = slideSymmetric 1.0 (PName "Test 1") (PWeight 1) Tune +p1 = slideSymmetric 1.0 (PName "Test 1") (pWeight 1) Tune p2 :: Proposal Double -p2 = slideSymmetric 1.0 (PName "Test 2") (PWeight 3) Tune +p2 = slideSymmetric 1.0 (PName "Test 2") (pWeight 3) Tune c :: Cycle Double c = cycleFromList [p1, p2] spec :: Spec spec = - describe "orderProposals" $ + describe "prepareProposals" $ it "returns the correct number of proposals in a cycle" $ do g <- create - l1 <- length <$> orderProposals c g + l1 <- length <$> prepareProposals c g l1 `shouldBe` 4 - l2 <- length <$> orderProposals (setOrder RandomReversibleO c) g + l2 <- length <$> prepareProposals (setOrder RandomReversibleO c) g l2 `shouldBe` 8 - o3 <- orderProposals (setOrder SequentialReversibleO c) g + o3 <- prepareProposals (setOrder SequentialReversibleO c) g length o3 `shouldBe` 8 o3 == [p1, p2, p2, p2, p2, p2, p2, p1] `shouldBe` True diff --git a/test/Mcmc/SaveSpec.hs b/test/Mcmc/SaveSpec.hs --- a/test/Mcmc/SaveSpec.hs +++ b/test/Mcmc/SaveSpec.hs @@ -1,7 +1,7 @@ -- | -- Module : Mcmc.SaveSpec -- Description : Unit tests for Mcmc.Save --- Copyright : (c) Dominik Schrempf, 2020 +-- Copyright : (c) Dominik Schrempf, 2021 -- License : GPL-3.0-or-later -- -- Maintainer : dominik.schrempf@gmail.com @@ -18,7 +18,6 @@ import Mcmc.Chain.Chain import Mcmc.Chain.Save import Mcmc.Chain.Trace -import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Normal import qualified System.Random.MWC as R @@ -36,10 +35,10 @@ proposals :: Cycle Double proposals = cycleFromList - [ slideSymmetric 0.1 (PName "Small") (PWeight 5) Tune, - slideSymmetric 1.0 (PName "Medium") (PWeight 2) Tune, - slideSymmetric 5.0 (PName "Large") (PWeight 2) Tune, - slide 1.0 4.0 (PName "Skewed") (PWeight 1) Tune + [ slideSymmetric 0.1 (PName "Small") (pWeight 5) Tune, + slideSymmetric 1.0 (PName "Medium") (pWeight 2) Tune, + slideSymmetric 5.0 (PName "Large") (pWeight 2) Tune, + slide 1.0 4.0 (PName "Skewed") (pWeight 1) Tune ] monStd :: MonitorStdOut Double @@ -62,8 +61,9 @@ Overwrite Sequential NoSave + LogStdOutOnly Quiet - c <- fromMHG <$> mhg noPrior lh proposals mon 0 gen + c <- fromMHG <$> mhg noPrior lh proposals mon TraceAuto 0 gen savedChain <- toSavedChain c c' <- fromSavedChain noPrior lh proposals mon savedChain putStrLn "@load . save@ should be @id@." @@ -87,7 +87,7 @@ g2' <- R.save $ generator r' g2 `shouldBe` g2' --- -- TODO. +-- -- TODO: 'mhContinue'. -- describe "mhContinue" -- $ it "mh 200 + mhContinue 200 == mh 400" -- $ do