diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -5,6 +5,18 @@ ## Unreleased changes +## 0.6.0.0 + +- Improve documentation. +- Generalized priors allowing automatic differentiation. +- Hamiltonian proposal. + + +### mcmc-tree + +- Moved to another repository: . + + ## 0.5.0.0 - Marginal likelihood estimation using thermodynamic integration or stepping diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -10,6 +10,7 @@ - Metropolis-Hastings-Green 1; - Metropolis-coupled Markov chain Monte Carlo (also known as parallel tempering) 2 , 3. +- Hamilton Monte Carlo proposal 4. ## Documentation @@ -20,8 +21,8 @@ ## Examples -[Example MCMC analyses](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples) can be built with [Stack](https://docs.haskellstack.org/en/stable/README/) and are attached to this -repository. +[Example MCMC analyses](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples) can be built with [cabal-install](https://cabal.readthedocs.io/en/latest/cabal-commands.html#) or [Stack](https://docs.haskellstack.org/en/stable/README/) and are attached +to this repository. git clone https://github.com/dschrempf/mcmc.git cd mcmc @@ -31,7 +32,9 @@ stack exec archery +For a more involved example, have a look at the [phylogenetic dating project](https://github.com/dschrempf/mcmc-dating). + # Footnotes 1 Geyer, C. J., Introduction to Markov chain Monte Carlo, In Handbook of @@ -39,8 +42,11 @@ 2 Geyer, C. J., Markov chain monte carlo maximum likelihood, Computing Science and Statistics, Proceedings of the 23rd Symposium on the Interface, -(), (1991). +(1991). 3 Altekar, G., Dwarkadas, S., Huelsenbeck, J. P., & Ronquist, F., Parallel metropolis coupled markov chain monte carlo for bayesian phylogenetic inference, Bioinformatics, 20(3), 407–415 (2004). + +4 Neal, R. M., Mcmc Using Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. +Jones, & X. Meng (Eds.), Handbook of Markov Chain Monte Carlo (2011). CRC press. diff --git a/bench/Normal.hs b/bench/Normal.hs --- a/bench/Normal.hs +++ b/bench/Normal.hs @@ -49,12 +49,13 @@ (AnalysisName "Normal") (BurnInWithAutoTuning 2000 500) (Iterations 20000) + TraceAuto Overwrite Sequential NoSave LogStdOutOnly Quiet - a <- mhg noPrior lh cc mon TraceAuto 0 g + a <- mhg s noPrior lh cc mon 0 g void $ mcmc s a ccLarge :: Cycle Double @@ -70,12 +71,13 @@ (AnalysisName "Normal") (BurnInWithAutoTuning 20 5) (Iterations 200) + TraceAuto Overwrite Sequential NoSave LogStdOutOnly Quiet - a <- mhg noPrior lh ccLarge mon TraceAuto 0 g + a <- mhg s noPrior lh ccLarge mon 0 g void $ mcmc s a ccBactrian :: Cycle Double @@ -88,12 +90,13 @@ (AnalysisName "NormalBactrian") (BurnInWithAutoTuning 2000 200) (Iterations 20000) + TraceAuto Overwrite Sequential NoSave LogStdOutOnly Quiet - a <- mhg noPrior lh ccBactrian mon TraceAuto 0 g + a <- mhg s noPrior lh ccBactrian mon 0 g void $ mcmc s a normalMC3 :: GenIO -> Int -> IO () @@ -103,11 +106,12 @@ (AnalysisName "MC3") (BurnInWithAutoTuning 200 20) (Iterations 2000) + TraceAuto Overwrite Sequential NoSave LogStdOutOnly Quiet mc3S = MC3Settings (NChains n) (SwapPeriod 2) (NSwaps 1) - a <- mc3 mc3S noPrior lh cc mon TraceAuto 0 g + a <- mc3 mc3S mcmcS noPrior lh cc mon 0 g void $ mcmc mcmcS a diff --git a/bench/Poisson.hs b/bench/Poisson.hs --- a/bench/Poisson.hs +++ b/bench/Poisson.hs @@ -32,7 +32,7 @@ ys = [1976.0 .. 1985.0] m = sum ys / fromIntegral (length ys) -f :: Int -> Double -> I -> Log Double +f :: Int -> Double -> I -> Likelihood f ft yr (a, b) = poisson l (fromIntegral ft) where l = exp $ a + b * yr @@ -71,10 +71,11 @@ (AnalysisName "Poisson") (BurnInWithAutoTuning 2000 200) (Iterations 10000) + TraceAuto Overwrite Sequential NoSave LogStdOutOnly Quiet - a <- mhg noPrior lh proposals mon TraceAuto initial g + a <- mhg s noPrior lh proposals mon initial g void $ mcmc s a diff --git a/mcmc.cabal b/mcmc.cabal --- a/mcmc.cabal +++ b/mcmc.cabal @@ -1,124 +1,139 @@ -cabal-version: 2.2 -name: mcmc -version: 0.5.0.0 -license: GPL-3.0-or-later -copyright: Dominik Schrempf (2021) -maintainer: dominik.schrempf@gmail.com -author: Dominik Schrempf -homepage: https://github.com/dschrempf/mcmc#readme -bug-reports: https://github.com/dschrempf/mcmc/issues -synopsis: Sample from a posterior using Markov chain Monte Carlo -description: - Please see the README on GitHub at +cabal-version: 2.2 +name: mcmc +version: 0.6.0.0 +synopsis: Sample from a posterior using Markov chain Monte Carlo +description: Please see the README on GitHub at +category: Math, Statistics +homepage: https://github.com/dschrempf/mcmc#readme +bug-reports: https://github.com/dschrempf/mcmc/issues +author: Dominik Schrempf +maintainer: dominik.schrempf@gmail.com +copyright: Dominik Schrempf (2021) +license: GPL-3.0-or-later +build-type: Simple -category: Math, Statistics -build-type: Simple extra-source-files: README.md ChangeLog.md source-repository head - type: git - location: https://github.com/dschrempf/mcmc + type: git + location: https://github.com/dschrempf/mcmc -library - exposed-modules: - Mcmc - Mcmc.Algorithm - 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 - Mcmc.Proposal.Generic - Mcmc.Proposal.Scale - Mcmc.Proposal.Slide - Mcmc.Proposal.Simplex - Mcmc.Settings - Mcmc.Statistics.Types +common lib + ghc-options: -Wall -Wunused-packages + default-language: Haskell2010 - hs-source-dirs: src - other-modules: - Mcmc.Internal.ByteString - Mcmc.Internal.Random - Mcmc.Internal.Shuffle - Paths_mcmc +common test + ghc-options: -Wall -Wunused-packages + default-language: Haskell2010 - autogen-modules: Paths_mcmc - default-language: Haskell2010 - ghc-options: -Wall -Wunused-packages - build-depends: - aeson >=1.5.6.0, - base >=4.7 && <5, - 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 +common bench + ghc-options: -Wall -Wunused-packages + default-language: Haskell2010 -test-suite mcmc-test - type: exitcode-stdio-1.0 - main-is: Spec.hs - hs-source-dirs: test - other-modules: - Mcmc.ProposalSpec - Mcmc.SaveSpec - Paths_mcmc +common exec + ghc-options: -Wall -Wunused-packages -threaded -rtsopts -with-rtsopts=-N + default-language: Haskell2010 - default-language: Haskell2010 - ghc-options: -Wall -Wunused-packages - build-depends: - base >=4.7 && <5, - hspec >=2.7.10, - log-domain >=0.13.1, - mcmc -any, - mwc-random >=0.15.0.1, - statistics >=0.15.2.0 +library + import: lib + exposed-modules: + Mcmc + Mcmc.Acceptance + Mcmc.Algorithm + Mcmc.Algorithm.MC3 + Mcmc.Algorithm.MHG + Mcmc.Chain.Chain + Mcmc.Chain.Link Mcmc.Chain.Save + Mcmc.Chain.Trace + Mcmc.Cycle + Mcmc.Environment + Mcmc.Internal.Gamma + 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 + Mcmc.Proposal.Generic + Mcmc.Proposal.Hamiltonian + Mcmc.Proposal.Scale + Mcmc.Proposal.Slide + Mcmc.Proposal.Simplex + Mcmc.Settings + Mcmc.Statistics.Types + other-modules: + Mcmc.Internal.ByteString + Mcmc.Internal.Random + Mcmc.Internal.Shuffle + Paths_mcmc + autogen-modules: + Paths_mcmc + hs-source-dirs: src + build-depends: + aeson + , base >=4.7 && <5 + , bytestring + , circular + , containers + , data-default + , deepseq + , directory + , dirichlet + , double-conversion + , log-domain + , math-functions + , matrices + , microlens + , mwc-random + , monad-parallel + , pretty-show + , primitive + , statistics + , time + , transformers + , vector + , zlib -benchmark mcmc-bench - type: exitcode-stdio-1.0 - main-is: Bench.hs - hs-source-dirs: bench - other-modules: - Normal - Poisson - Paths_mcmc +test-suite mcmc-test + import: test + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + Mcmc.ProposalSpec + Mcmc.SaveSpec + Paths_mcmc + hs-source-dirs: test + build-depends: + base >=4.7 && <5 + , hspec + , log-domain + , mcmc + , mwc-random + , statistics - default-language: Haskell2010 - ghc-options: -Wall -Wunused-packages - build-depends: - base >=4.7 && <5, - criterion >=1.5.9.0, - log-domain >=0.13.1, - mcmc -any, - microlens >=0.4.12.0, - mwc-random >=0.15.0.1 +benchmark mcmc-bench + import: bench + type: exitcode-stdio-1.0 + main-is: Bench.hs + other-modules: + Normal + Poisson + Paths_mcmc + autogen-modules: + Paths_mcmc + hs-source-dirs: bench + build-depends: + base >=4.7 && <5 + , criterion + , mcmc + , microlens + , mwc-random diff --git a/src/Mcmc.hs b/src/Mcmc.hs --- a/src/Mcmc.hs +++ b/src/Mcmc.hs @@ -134,7 +134,10 @@ slideUniformSymmetric, slideContrarily, slideBactrian, + module Mcmc.Proposal.Hamiltonian, module Mcmc.Proposal.Simplex, + + -- ** Cycles Cycle, cycleFromList, Order (..), @@ -172,6 +175,7 @@ -- * Prior, likelihood, and posterior values and functions module Mcmc.Prior, module Mcmc.Likelihood, + module Mcmc.Posterior, -- * MCMC samplers mcmc, @@ -195,16 +199,18 @@ import Mcmc.Algorithm.MC3 import Mcmc.Algorithm.MHG -import Mcmc.Chain.Chain +import Mcmc.Cycle import Mcmc.Likelihood import Mcmc.MarginalLikelihood import Mcmc.Mcmc import Mcmc.Monitor import Mcmc.Monitor.Parameter import Mcmc.Monitor.ParameterBatch +import Mcmc.Posterior import Mcmc.Prior import Mcmc.Proposal import Mcmc.Proposal.Bactrian +import Mcmc.Proposal.Hamiltonian import Mcmc.Proposal.Scale import Mcmc.Proposal.Simplex import Mcmc.Proposal.Slide diff --git a/src/Mcmc/Acceptance.hs b/src/Mcmc/Acceptance.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Acceptance.hs @@ -0,0 +1,93 @@ +-- | +-- Module : Mcmc.Acceptance +-- Description : Handle acceptance rates +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Thu Jul 8 18:12:07 2021. +module Mcmc.Acceptance + ( -- * Acceptance rates + AcceptanceRate, + Acceptance (fromAcceptance), + emptyA, + pushA, + resetA, + transformKeysA, + acceptanceRate, + acceptanceRates, + ) +where + +import Control.DeepSeq +import Data.Aeson +import Data.Bifunctor +import Data.Foldable +import qualified Data.Map.Strict as M + +-- | Acceptance rate. +type AcceptanceRate = Double + +-- | 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, AcceptanceRate) +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 AcceptanceRate) +acceptanceRates = + M.map + ( \(as, rs) -> + if as + rs == 0 + then Nothing + else Just $ fromIntegral as / fromIntegral (as + rs) + ) + . fromAcceptance diff --git a/src/Mcmc/Algorithm.hs b/src/Mcmc/Algorithm.hs --- a/src/Mcmc/Algorithm.hs +++ b/src/Mcmc/Algorithm.hs @@ -33,8 +33,11 @@ -- | Sample the next state. aIterate :: ParallelizationMode -> a -> IO a - -- | Auto tune all proposals. - aAutoTune :: a -> a + -- | Auto tune all proposals over the last N iterations. + -- + -- NOTE: Computation in the 'IO' Monad is necessary because the trace is + -- mutable. + aAutoTune :: Int -> a -> IO a -- | Reset acceptance counts. aResetAcceptance :: a -> 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 @@ -62,17 +62,20 @@ import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U import Data.Word --- import Debug.Trace hiding (trace) +import Mcmc.Acceptance import Mcmc.Algorithm import Mcmc.Algorithm.MHG import Mcmc.Chain.Chain import Mcmc.Chain.Link import Mcmc.Chain.Save import Mcmc.Chain.Trace +import Mcmc.Cycle import Mcmc.Internal.Random import Mcmc.Internal.Shuffle +import Mcmc.Likelihood import Mcmc.Monitor import Mcmc.Posterior +import Mcmc.Prior import Mcmc.Proposal import Mcmc.Settings import Numeric.Log hiding (sum) @@ -275,28 +278,28 @@ -- - The swap period is zero or negative. mc3 :: MC3Settings -> + Settings -> PriorFunction a -> LikelihoodFunction a -> Cycle a -> Monitor a -> - TraceLength -> InitialState a -> GenIO -> IO (MC3 a) -mc3 s pr lh cc mn tr i0 g +mc3 sMc3 s pr lh cc mn 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 tr i0) gs + cs <- V.mapM (mhg s pr lh cc mn i0) gs hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs - return $ MC3 s hcs bs 0 (emptyA [0 .. n - 2]) g + return $ MC3 sMc3 hcs bs 0 (emptyA [0 .. n - 2]) g where - n = fromNChains $ mc3NChains s - sp = fromSwapPeriod $ mc3SwapPeriod s - sn = fromNSwaps $ mc3NSwaps s + n = fromNChains $ mc3NChains sMc3 + sp = fromSwapPeriod $ mc3SwapPeriod sMc3 + sn = fromNSwaps $ mc3NSwaps sMc3 -- NOTE: The initial choice of reciprocal temperatures is based on a few -- tests but otherwise pretty arbitrary. -- @@ -417,10 +420,10 @@ where mhgs = mc3MHGChains a --- TODO: Splimix. 'mc3Iterate' is actually not parallel, but concurrent because +-- TODO: Splitmix. 'mc3Iterate' is actually not parallel, but concurrent because -- of the IO constraint. Use pure parallel code when we have a pure generator. -- --- However, we to have honor the mutable traces. +-- However, we have to take care of the mutable traces. mc3Iterate :: ToJSON a => ParallelizationMode -> @@ -472,41 +475,40 @@ rNew = (brOld / blOld) ** xi brNew = blNew * rNew -mc3AutoTune :: ToJSON a => MC3 a -> MC3 a -mc3AutoTune a = a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'} - where - mhgs = mc3MHGChains a - -- 1. Auto tune all chains. - mhgs' = V.map aAutoTune mhgs - -- 2. Auto tune temperatures. - optimalRate = getOptimalRate PDimensionUnknown - 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 - -- ratios, which is the same. Also, by working with ratios in (0,1) of - -- neighboring chains, we ensure the monotonicity of the reciprocal - -- temperatures. - -- - -- The factor (1/2) was determined by a few tests and is otherwise - -- absolutely arbitrary. - 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 - -- the cold chain. - bs' = foldl' (\xs j -> tuneBeta bs j (xi j) xs) bs [0 .. n - 2] - coldChain = fromMHG $ V.head mhgs' - coldPrF = priorFunction coldChain - coldLhF = likelihoodFunction coldChain - mhgs'' = - V.head mhgs' - `V.cons` V.zipWith - (setReciprocalTemperature coldPrF coldLhF) - (V.convert $ U.tail bs') - (V.tail mhgs') +mc3AutoTune :: ToJSON a => Int -> MC3 a -> IO (MC3 a) +mc3AutoTune l a = do + -- 1. Auto tune all chains. + mhgs' <- V.mapM (aAutoTune l) $ mc3MHGChains a + -- 2. Auto tune temperatures. + let optimalRate = getOptimalRate PDimensionUnknown + 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 + -- ratios, which is the same. Also, by working with ratios in (0,1) of + -- neighboring chains, we ensure the monotonicity of the reciprocal + -- temperatures. + -- + -- The factor (1/2) was determined by a few tests and is otherwise + -- absolutely arbitrary. + 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 + -- the cold chain. + bs' = foldl' (\xs j -> tuneBeta bs j (xi j) xs) bs [0 .. n - 2] + coldChain = fromMHG $ V.head mhgs' + coldPrF = priorFunction coldChain + coldLhF = likelihoodFunction coldChain + mhgs'' = + V.head mhgs' + `V.cons` V.zipWith + (setReciprocalTemperature coldPrF coldLhF) + (V.convert $ U.tail bs') + (V.tail mhgs') + return $ a {mc3MHGChains = mhgs'', mc3ReciprocalTemperatures = bs'} mc3ResetAcceptance :: ToJSON a => MC3 a -> MC3 a mc3ResetAcceptance a = a' @@ -536,6 +538,7 @@ 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.", diff --git a/src/Mcmc/Algorithm/MHG.hs b/src/Mcmc/Algorithm/MHG.hs --- a/src/Mcmc/Algorithm/MHG.hs +++ b/src/Mcmc/Algorithm/MHG.hs @@ -33,13 +33,18 @@ import Data.Aeson import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time +import qualified Data.Vector as VB +import Mcmc.Acceptance import Mcmc.Algorithm import Mcmc.Chain.Chain import Mcmc.Chain.Link import Mcmc.Chain.Save import Mcmc.Chain.Trace +import Mcmc.Cycle +import Mcmc.Likelihood import Mcmc.Monitor import Mcmc.Posterior +import Mcmc.Prior hiding (uniform) import Mcmc.Proposal import Mcmc.Settings import Numeric.Log @@ -64,19 +69,20 @@ aCloseMonitors = mhgCloseMonitors aSave = mhgSave --- NOTE: IO is required because the trace is mutable. - -- | Initialize an MHG algorithm. +-- +-- NOTE: Computation in the 'IO' Monad is necessary because the trace is +-- mutable. mhg :: + Settings -> PriorFunction a -> LikelihoodFunction a -> Cycle a -> Monitor a -> - TraceLength -> InitialState a -> GenIO -> IO (MHG a) -mhg pr lh cc mn trLen i0 g = do +mhg s 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 @@ -85,10 +91,14 @@ l0 = Link i0 (pr i0) (lh i0) ac = emptyA $ ccProposals cc batchMonitorSizes = map getMonitorBatchSize $ mBatches mn - minimumTraceLength = case trLen of + minimumTraceLength = case sTraceLength s of TraceAuto -> 1 TraceMinimum n -> n - traceLength = maximum $ minimumTraceLength : batchMonitorSizes + bi = case sBurnIn s of + BurnInWithAutoTuning _ n -> n + BurnInWithCustomAutoTuning ns -> maximum ns + _ -> 0 + traceLength = maximum $ minimumTraceLength : bi : batchMonitorSizes mhgFn :: AnalysisName -> FilePath mhgFn (AnalysisName nm) = nm ++ ".mcmc.mhg" @@ -210,12 +220,17 @@ t = trace c n = iteration c +-- Check if the current state is invalid. +-- +-- At the moment this just checks whether the prior, likelihood, or posterior +-- are NaN or infinite. mhgIsInValidState :: MHG a -> Bool -mhgIsInValidState a = ((p * l) == 0) || (p * l == (0 / 0)) +mhgIsInValidState a = check p || check l || check (p * l) where x = link $ fromMHG a p = prior x l = likelihood x + check v = let v' = ln v in isNaN v' || isInfinite v' || v' == 0 -- 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 @@ -231,8 +246,10 @@ cc = cycle c g = generator c -mhgAutoTune :: MHG a -> MHG a -mhgAutoTune (MHG c) = MHG $ c {cycle = autoTuneCycle ac cc} +mhgAutoTune :: Int -> MHG a -> IO (MHG a) +mhgAutoTune n (MHG c) = do + tr <- VB.map state <$> takeT n (trace c) + return $ MHG $ c {cycle = autoTuneCycle ac tr cc} where ac = acceptance c cc = cycle 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 @@ -10,24 +10,51 @@ -- -- Creation date: Tue May 5 18:01:15 2020. module Mcmc.Chain.Chain - ( PriorFunction, - noPrior, - LikelihoodFunction, - noLikelihood, - InitialState, + ( InitialState, Chain (..), ) where --- 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. +-- NOTE: Auxiliary data. +-- +-- It is not necessary to add another type @b@ to store auxiliary data about the +-- chain. The information can just be stored in @a@ equally well. +-- +-- I am not sure if this is the case. If @b@ is affected by a change in @a@, it +-- has to be recomputed. This is difficult to implement at the proposal step. +-- Maybe the new value is not even accepted. On the other hand, one could trust +-- the lazyness of Haskell, and recompute @b@. The computation is done only when +-- the value is accessed. +-- +-- That is, I have to think about how to implement auxiliary data A. Prior and +-- likelihood functions can then act on the state space I and A, e.g., +-- +-- > type PriorFunction a b = a -> b -> Log Double +-- +-- where a is the type of the state, and b is the auxiliary data type. --- TODO: First class parameters. Make a type class for parameter types (name, --- lens, proposals, monitors). +-- NOTE: First class parameters. +-- +-- I thought a lot about implementing a type class for parameters. For example, +-- parameters should have a name, a lens, possibly proposals, monitors, and so +-- on. +-- +-- However, this is really difficult. If I use a type class, I need different +-- data types for each parameter which is cumbersome (and slow?). Of course, one +-- could use a data type such as +-- +-- > data ParamterSpec a = ParameterSpec { name :: ByteString, pMonitor :: (a -> ByteString) } +-- +-- But even in this case we run into problems: There are proposals and monitors +-- acting on a combination of parameters. Even setting the name doesn't make +-- sense in this case. +-- +-- I decided to let this idea rest. +import Mcmc.Acceptance import Mcmc.Chain.Link import Mcmc.Chain.Trace +import Mcmc.Cycle import Mcmc.Likelihood import Mcmc.Monitor import Mcmc.Prior 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 @@ -27,17 +27,20 @@ import Data.Aeson import Data.Aeson.TH 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 +import Mcmc.Acceptance import Mcmc.Chain.Chain import Mcmc.Chain.Link import Mcmc.Chain.Trace +import Mcmc.Cycle import Mcmc.Internal.Random +import Mcmc.Likelihood import Mcmc.Monitor +import Mcmc.Prior import Mcmc.Proposal import Prelude hiding (cycle) @@ -51,7 +54,7 @@ savedTrace :: C.Stack VB.Vector (Link a), savedAcceptance :: Acceptance Int, savedSeed :: VU.Vector Word32, - savedTuningParameters :: [Maybe TuningParameter] + savedTuningParameters :: [Maybe (TuningParameter, AuxiliaryTuningParameters)] } deriving (Eq, Read, Show) @@ -68,7 +71,10 @@ where ps = ccProposals cc ac' = transformKeysA ps [0 ..] ac - ts = [fmap tParam mt | mt <- map prTuner ps] + ts = + [ (\t -> (tGetTuningParameter t, tGetAuxiliaryTuningParameters t)) <$> mt + | mt <- map prTuner ps + ] -- | Load a saved chain. -- @@ -100,13 +106,9 @@ 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 - Nothing -> const 1.0 - Just t -> const t - cc' = - tuneCycle - ( M.map getTuningF $ - M.fromList $ - zip (ccProposals cc) ts - ) - cc + tunePs mt p = case mt of + Nothing -> p + Just (x, xs) -> either (error . (<> err)) id $ tuneWithTuningParameters x xs p + err = error "\nfromSavedChain: Proposal with stored tuning parameters is not tunable." + ps = ccProposals cc + cc' = cc {ccProposals = zipWith tunePs ts ps} diff --git a/src/Mcmc/Cycle.hs b/src/Mcmc/Cycle.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Cycle.hs @@ -0,0 +1,187 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE OverloadedStrings #-} + +-- | +-- Module : Mcmc.Cycle +-- Description : A cycle is a list of proposals +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Thu Jul 8 17:56:03 2021. +module Mcmc.Cycle + ( -- * Cycles + Order (..), + Cycle (ccProposals), + cycleFromList, + setOrder, + prepareProposals, + autoTuneCycle, + + -- * Output + proposalHLine, + summarizeCycle, + ) +where + +import qualified Data.ByteString.Builder as BB +import qualified Data.ByteString.Lazy.Char8 as BL +import Data.Default +import Data.Either +import Data.List +import qualified Data.Map.Strict as M +import qualified Data.Vector as VB +import Mcmc.Acceptance +import Mcmc.Internal.Shuffle +import Mcmc.Proposal +import System.Random.MWC + +-- | Define the order in which 'Proposal's are executed in a 'Cycle'. The total +-- number of 'Proposal's per 'Cycle' may differ between 'Order's (e.g., compare +-- 'RandomO' and 'RandomReversibleO'). +data Order + = -- | Shuffle the 'Proposal's in the 'Cycle'. The 'Proposal's are replicated + -- according to their weights and executed in random order. If a 'Proposal' has + -- weight @w@, it is executed exactly @w@ times per iteration. + RandomO + | -- | The 'Proposal's are executed sequentially, in the order they appear in the + -- 'Cycle'. 'Proposal's with weight @w>1@ are repeated immediately @w@ times + -- (and not appended to the end of the list). + SequentialO + | -- | Similar to 'RandomO'. However, a reversed copy of the list of + -- shuffled 'Proposal's is appended such that the resulting Markov chain is + -- reversible. + -- Note: the total number of 'Proposal's executed per cycle is twice the number + -- of 'RandomO'. + RandomReversibleO + | -- | Similar to 'SequentialO'. However, a reversed copy of the list of + -- sequentially ordered 'Proposal's is appended such that the resulting Markov + -- chain is reversible. + SequentialReversibleO + deriving (Eq, Show) + +instance Default Order where def = RandomO + +-- Describe the order. +describeOrder :: Order -> BL.ByteString +describeOrder RandomO = "The proposals are executed in random order." +describeOrder SequentialO = "The proposals are executed sequentially." +describeOrder RandomReversibleO = + BL.intercalate + "\n" + [ describeOrder RandomO, + "A reversed copy of the shuffled proposals is appended to ensure reversibility." + ] +describeOrder SequentialReversibleO = + BL.intercalate + "\n" + [ describeOrder SequentialO, + "A reversed copy of the sequential proposals is appended to ensure reversibility." + ] + +-- | In brief, a 'Cycle' is a list of proposals. +-- +-- The state of the Markov chain will be logged only after all 'Proposal's in +-- the 'Cycle' have been completed, and the iteration counter will be increased +-- by one. The order in which the 'Proposal's are executed is specified by +-- 'Order'. The default is 'RandomO'. +-- +-- No proposals with the same name and description are allowed in a 'Cycle', so +-- that they can be uniquely identified. +data Cycle a = Cycle + { ccProposals :: [Proposal a], + ccOrder :: Order + } + +-- | Create a 'Cycle' from a list of 'Proposal's. +cycleFromList :: [Proposal a] -> Cycle a +cycleFromList [] = + error "cycleFromList: Received an empty list but cannot create an empty Cycle." +cycleFromList xs = + if length uniqueXs == length xs + then Cycle xs def + else error $ "\n" ++ msg ++ "cycleFromList: Proposals are not unique." + where + uniqueXs = nub xs + removedXs = xs \\ uniqueXs + removedNames = map (show . prName) removedXs + removedDescriptions = map (show . prDescription) removedXs + removedMsgs = zipWith (\n d -> n ++ " " ++ d) removedNames removedDescriptions + msg = unlines removedMsgs + +-- | Set the order of 'Proposal's in a 'Cycle'. +setOrder :: Order -> Cycle a -> Cycle a +setOrder o c = c {ccOrder = o} + +-- | Replicate 'Proposal's according to their weights and possibly shuffle them. +prepareProposals :: Cycle a -> GenIO -> IO [Proposal a] +prepareProposals (Cycle xs o) g = case o of + RandomO -> shuffle ps g + SequentialO -> return ps + RandomReversibleO -> do + psR <- shuffle ps g + return $ psR ++ reverse psR + SequentialReversibleO -> return $ ps ++ reverse ps + where + !ps = concat [replicate (fromPWeight $ prWeight p) p | p <- xs] + +-- The number of proposals depends on the order. +getNProposalsPerCycle :: Cycle a -> Int +getNProposalsPerCycle (Cycle xs o) = case o of + RandomO -> once + SequentialO -> once + RandomReversibleO -> 2 * once + SequentialReversibleO -> 2 * once + where + once = sum $ map (fromPWeight . prWeight) xs + +-- | Calculate acceptance rates and auto tunes the 'Proposal's in the 'Cycle'. +-- +-- Do not change 'Proposal's that are not tuneable. +autoTuneCycle :: Acceptance (Proposal a) -> VB.Vector a -> Cycle a -> Cycle a +autoTuneCycle a xs c = + if sort (M.keys ar) == sort ps + then c {ccProposals = map tuneF ps} + else error "autoTuneCycle: Proposals in map and cycle do not match." + where + ar = acceptanceRates a + ps = ccProposals c + tuneF p = case ar M.!? p of + Just (Just x) -> fromRight p (tuneWithChainParameters x xs p) + _ -> p + +-- | Horizontal line of proposal summaries. +proposalHLine :: BL.ByteString +proposalHLine = BL.replicate (BL.length proposalHeader) '-' + +-- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance rates. +summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString +summarizeCycle a c = + BL.intercalate "\n" $ + [ "Summary of proposal(s) in cycle.", + nProposalsFullStr, + describeOrder (ccOrder c), + proposalHeader, + proposalHLine + ] + ++ [ summarizeProposal + (prName p) + (prDescription p) + (prWeight p) + (tGetTuningParameter <$> prTuner p) + (prDimension p) + (ar p) + | p <- ps + ] + ++ [proposalHLine] + where + ps = ccProposals c + nProposals = getNProposalsPerCycle c + nProposalsStr = BB.toLazyByteString $ BB.intDec nProposals + nProposalsFullStr = case nProposals of + 1 -> nProposalsStr <> " proposal is performed per iteration." + _ -> nProposalsStr <> " proposals are performed per iterations." + ar m = acceptanceRate m a diff --git a/src/Mcmc/Environment.hs b/src/Mcmc/Environment.hs --- a/src/Mcmc/Environment.hs +++ b/src/Mcmc/Environment.hs @@ -18,7 +18,7 @@ import Control.Concurrent.MVar import Control.Monad -import Data.Time.Clock +import Data.Time import Mcmc.Logger import Mcmc.Settings import System.IO diff --git a/src/Mcmc/Internal/Gamma.hs b/src/Mcmc/Internal/Gamma.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Internal/Gamma.hs @@ -0,0 +1,195 @@ +-- | +-- Module : Mcmc.Internal.Gamma +-- Description : Generalized gamma function for automatic differentiation +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Tue Jul 13 12:53:09 2021. +-- +-- The code is taken from "Numeric.SpecFunctions". +module Mcmc.Internal.Gamma + ( logGammaG, + ) +where + +import qualified Data.Vector as VB +import Numeric.Polynomial + +mSqrtEps :: RealFloat a => a +mSqrtEps = 1.4901161193847656e-8 + +mEulerMascheroni :: RealFloat a => a +mEulerMascheroni = 0.5772156649015328606065121 + +-- | See 'Numeric.SpecFunctions.logGamma'. +logGammaG :: RealFloat a => a -> a +logGammaG z + | z <= 0 = 1 / 0 + | z < mSqrtEps = log (1 / z - mEulerMascheroni) + | z < 0.5 = lgamma1_15 z (z - 1) - log z + | z < 1 = lgamma15_2 z (z - 1) - log z + | z <= 1.5 = lgamma1_15 (z - 1) (z - 2) + | z < 2 = lgamma15_2 (z - 1) (z - 2) + | z < 15 = lgammaSmall z + | otherwise = lanczosApprox z +{-# SPECIALIZE logGammaG :: Double -> Double #-} + +lgamma1_15 :: RealFloat a => a -> a -> a +lgamma1_15 zm1 zm2 = + r * y + r + * ( evaluatePolynomial zm1 tableLogGamma_1_15P + / evaluatePolynomial zm1 tableLogGamma_1_15Q + ) + where + r = zm1 * zm2 + y = 0.52815341949462890625 + +tableLogGamma_1_15P :: RealFloat a => VB.Vector a +tableLogGamma_1_15P = + VB.fromList + [ 0.490622454069039543534e-1, + -0.969117530159521214579e-1, + -0.414983358359495381969e0, + -0.406567124211938417342e0, + -0.158413586390692192217e0, + -0.240149820648571559892e-1, + -0.100346687696279557415e-2 + ] +{-# NOINLINE tableLogGamma_1_15P #-} + +tableLogGamma_1_15Q :: RealFloat a => VB.Vector a +tableLogGamma_1_15Q = + VB.fromList + [ 1, + 0.302349829846463038743e1, + 0.348739585360723852576e1, + 0.191415588274426679201e1, + 0.507137738614363510846e0, + 0.577039722690451849648e-1, + 0.195768102601107189171e-2 + ] +{-# NOINLINE tableLogGamma_1_15Q #-} + +lgamma15_2 :: RealFloat a => a -> a -> a +lgamma15_2 zm1 zm2 = + r * y + r + * ( evaluatePolynomial (- zm2) tableLogGamma_15_2P + / evaluatePolynomial (- zm2) tableLogGamma_15_2Q + ) + where + r = zm1 * zm2 + y = 0.452017307281494140625 + +tableLogGamma_15_2P :: RealFloat a => VB.Vector a +tableLogGamma_15_2P = + VB.fromList + [ -0.292329721830270012337e-1, + 0.144216267757192309184e0, + -0.142440390738631274135e0, + 0.542809694055053558157e-1, + -0.850535976868336437746e-2, + 0.431171342679297331241e-3 + ] +{-# NOINLINE tableLogGamma_15_2P #-} + +tableLogGamma_15_2Q :: RealFloat a => VB.Vector a +tableLogGamma_15_2Q = + VB.fromList + [ 1, + -0.150169356054485044494e1, + 0.846973248876495016101e0, + -0.220095151814995745555e0, + 0.25582797155975869989e-1, + -0.100666795539143372762e-2, + -0.827193521891290553639e-6 + ] +{-# NOINLINE tableLogGamma_15_2Q #-} + +lgammaSmall :: RealFloat a => a -> a +lgammaSmall = go 0 + where + go acc z + | z < 3 = acc + lgamma2_3 z + | otherwise = go (acc + log zm1) zm1 + where + zm1 = z - 1 + +lgamma2_3 :: RealFloat a => a -> a +lgamma2_3 z = + r * y + r + * ( evaluatePolynomial zm2 tableLogGamma_2_3P + / evaluatePolynomial zm2 tableLogGamma_2_3Q + ) + where + r = zm2 * (z + 1) + zm2 = z - 2 + y = 0.158963680267333984375e0 + +tableLogGamma_2_3P :: RealFloat a => VB.Vector a +tableLogGamma_2_3P = + VB.fromList + [ -0.180355685678449379109e-1, + 0.25126649619989678683e-1, + 0.494103151567532234274e-1, + 0.172491608709613993966e-1, + -0.259453563205438108893e-3, + -0.541009869215204396339e-3, + -0.324588649825948492091e-4 + ] +{-# NOINLINE tableLogGamma_2_3P #-} + +tableLogGamma_2_3Q :: RealFloat a => VB.Vector a +tableLogGamma_2_3Q = + VB.fromList + [ 1, + 0.196202987197795200688e1, + 0.148019669424231326694e1, + 0.541391432071720958364e0, + 0.988504251128010129477e-1, + 0.82130967464889339326e-2, + 0.224936291922115757597e-3, + -0.223352763208617092964e-6 + ] +{-# NOINLINE tableLogGamma_2_3Q #-} + +lanczosApprox :: RealFloat a => a -> a +lanczosApprox z = + (log (z + g - 0.5) - 1) * (z - 0.5) + + log (evalRatio tableLanczos z) + where + g = 6.024680040776729583740234375 + +tableLanczos :: RealFloat a => VB.Vector (a, a) +tableLanczos = + VB.fromList + [ (56906521.91347156388090791033559122686859, 0), + (103794043.1163445451906271053616070238554, 39916800), + (86363131.28813859145546927288977868422342, 120543840), + (43338889.32467613834773723740590533316085, 150917976), + (14605578.08768506808414169982791359218571, 105258076), + (3481712.15498064590882071018964774556468, 45995730), + (601859.6171681098786670226533699352302507, 13339535), + (75999.29304014542649875303443598909137092, 2637558), + (6955.999602515376140356310115515198987526, 357423), + (449.9445569063168119446858607650988409623, 32670), + (19.51992788247617482847860966235652136208, 1925), + (0.5098416655656676188125178644804694509993, 66), + (0.006061842346248906525783753964555936883222, 1) + ] +{-# NOINLINE tableLanczos #-} + +data L a = L !a !a + +evalRatio :: RealFloat a => VB.Vector (a, a) -> a -> a +evalRatio coef x + | x > 1 = fini $ VB.foldl' stepL (L 0 0) coef + | otherwise = fini $ VB.foldr' stepR (L 0 0) coef + where + fini (L num den) = num / den + stepR (a, b) (L num den) = L (num * x + a) (den * x + b) + stepL (L num den) (a, b) = L (num * rx + a) (den * rx + b) + rx = recip x diff --git a/src/Mcmc/Likelihood.hs b/src/Mcmc/Likelihood.hs --- a/src/Mcmc/Likelihood.hs +++ b/src/Mcmc/Likelihood.hs @@ -11,7 +11,9 @@ -- Creation date: Wed Mar 3 11:39:04 2021. module Mcmc.Likelihood ( Likelihood, + LikelihoodG, LikelihoodFunction, + LikelihoodFunctionG, noLikelihood, ) where @@ -21,9 +23,16 @@ -- | Likelihood values are stored in log domain. type Likelihood = Log Double +-- | Generalized likelihood. +type LikelihoodG a = Log a + -- | Likelihood function. -type LikelihoodFunction a = a -> Log Double +type LikelihoodFunction a = LikelihoodFunctionG a Double +-- | Generalized likelihood function. +type LikelihoodFunctionG a b = a -> LikelihoodG b + -- | Flat likelihood function. Useful for testing and debugging. -noLikelihood :: LikelihoodFunction a +noLikelihood :: RealFloat b => LikelihoodFunctionG a b noLikelihood = const 1.0 +{-# SPECIALIZE noLikelihood :: LikelihoodFunction Double #-} diff --git a/src/Mcmc/Logger.hs b/src/Mcmc/Logger.hs --- a/src/Mcmc/Logger.hs +++ b/src/Mcmc/Logger.hs @@ -42,16 +42,14 @@ import Data.Aeson.TH import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time.Clock +import Data.Version (showVersion) 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 +import System.IO -- | Define where the log output should be directed to. +-- +-- Logging is disabled if 'Verbosity' is set to 'Quiet'. data LogMode = LogStdOutAndFile | LogStdOutOnly | LogFileOnly deriving (Eq, Read, Show) diff --git a/src/Mcmc/MarginalLikelihood.hs b/src/Mcmc/MarginalLikelihood.hs --- a/src/Mcmc/MarginalLikelihood.hs +++ b/src/Mcmc/MarginalLikelihood.hs @@ -30,17 +30,19 @@ import qualified Data.Map.Strict as M import qualified Data.Vector as VB import qualified Data.Vector.Unboxed as VU +import Mcmc.Acceptance import Mcmc.Algorithm.MHG import Mcmc.Chain.Chain import Mcmc.Chain.Link import Mcmc.Chain.Trace +import Mcmc.Cycle import Mcmc.Environment import Mcmc.Internal.Random import Mcmc.Likelihood import Mcmc.Logger import Mcmc.Mcmc import Mcmc.Monitor -import Mcmc.Proposal +import Mcmc.Prior import Mcmc.Settings import Numeric.Log hiding (sum) import System.Directory @@ -87,9 +89,9 @@ mlAlgorithm :: MLAlgorithm, mlNPoints :: NPoints, -- | Initial burn in at the starting point of the path. - mlInitialBurnIn :: BurnInSpecification, + mlInitialBurnIn :: BurnInSettings, -- | Repetitive burn in at each point on the path. - mlPointBurnIn :: BurnInSpecification, + mlPointBurnIn :: BurnInSettings, -- | The number of iterations performed at each point. mlIterations :: Iterations, mlExecutionMode :: ExecutionMode, @@ -235,11 +237,11 @@ 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 + ssI = Settings nm biI (Iterations 0) trLen em Sequential NoSave LogFileOnly vb' + ssP = Settings nm biP is trLen em Sequential NoSave LogFileOnly vb' logDebugB "mlRun: Initialize MHG algorithm." - a0 <- liftIO $ mhg prf lhf cc mn trLen i0 g + a0 <- liftIO $ mhg ssI prf lhf cc mn i0 g logDebugS $ "mlRun: Perform initial burn in at first point " <> show x0 <> "." a1 <- sampleAtPoint x0 ssI lhf a0 logDebugB "mlRun: Traverse points." diff --git a/src/Mcmc/Mcmc.hs b/src/Mcmc/Mcmc.hs --- a/src/Mcmc/Mcmc.hs +++ b/src/Mcmc/Mcmc.hs @@ -21,6 +21,7 @@ ) where +import Control.Exception import Control.Monad import Control.Monad.IO.Class import Control.Monad.Trans.Reader @@ -28,13 +29,13 @@ import Mcmc.Environment import Mcmc.Logger import Mcmc.Settings +import System.Exit import System.IO -import Text.Show.Pretty import Prelude hiding (cycle) -- The MCMC algorithm has read access to an environment and uses an algorithm -- transforming the state @a@. -type MCMC a = ReaderT (Environment Settings) IO a +type MCMC = ReaderT (Environment Settings) IO mcmcExecute :: Algorithm a => a -> MCMC a mcmcExecute a = do @@ -52,6 +53,30 @@ logDebugB "Reset acceptance rates." return $ aResetAcceptance a +mcmcExceptionHandler :: Algorithm a => Environment Settings -> a -> AsyncException -> IO b +mcmcExceptionHandler e a UserInterrupt = do + putStrLn "" + putStrLn "USER INTERRUPT!" + putStrLn "Try to terminate gracefully and save chain for continuation." + putStrLn "Press CTRL-C again to terminate now." + putStrLn "Close output files." + _ <- aCloseMonitors a + closeEnvironment e + putStrLn "Try to save settings." + let s = settings e + settingsSave s + putStrLn "Try to save compressed MCMC analysis." + putStrLn "For long traces, or complex objects, this may take a while." + let nm = sAnalysisName s + aSave nm a + putStrLn "Markov chain saved. Analysis can be continued." + putStrLn "Terminate gracefully." + exitWith $ ExitFailure 1 +mcmcExceptionHandler _ _ e = throw e + +-- XXX: Exception handling. Is it enough to mask execution of monitors and catch +-- UserInterrupt during iterations? + mcmcExecuteMonitors :: Algorithm a => a -> MCMC () mcmcExecuteMonitors a = do e <- ask @@ -59,7 +84,8 @@ vb = sVerbosity s t0 = startingTime e iTotal = burnInIterations (sBurnIn s) + fromIterations (sIterations s) - mStdLog <- liftIO (aExecuteMonitors vb t0 iTotal a) + -- NOTE: Mask asynchronous exceptions when writing monitor files. + mStdLog <- liftIO $ mask_ $ aExecuteMonitors vb t0 iTotal a forM_ mStdLog (logOutB " ") mcmcIterate :: Algorithm a => Int -> a -> MCMC a @@ -67,8 +93,10 @@ | n < 0 = error "mcmcIterate: Number of iterations is negative." | n == 0 = return a | otherwise = do + e <- ask p <- sParallelizationMode . settings <$> ask - a' <- liftIO $ aIterate p a + -- NOTE: User interrupt is handled during iterations. + a' <- liftIO $ catch (aIterate p a) (mcmcExceptionHandler e a) mcmcExecuteMonitors a' mcmcIterate (n -1) a' @@ -90,11 +118,16 @@ mcmcContinueRun :: Algorithm a => a -> MCMC a mcmcContinueRun a = do s <- reader settings - let iTotal = fromIterations (sIterations s) + burnInIterations (sBurnIn s) + let iBurnIn = burnInIterations (sBurnIn s) + iNormal = fromIterations (sIterations s) + iTotal = iBurnIn + iNormal logInfoB "Continuation of MCMC sampler." let iCurrent = aIteration a - logInfoS $ "Current iteration: " ++ show iCurrent ++ "." + logInfoS $ "Burn in iterations: " ++ show iBurnIn ++ "." + logInfoS $ "Normal iterations: " ++ show iNormal ++ "." logInfoS $ "Total iterations: " ++ show iTotal ++ "." + logInfoS $ "Current iteration: " ++ show iCurrent ++ "." + when (iCurrent < iBurnIn) $ error "mcmcContinueRun: Can not continue burn in." let di = iTotal - iCurrent logInfoB $ aSummarizeCycle a logInfoS $ "Run chain for " ++ show di ++ " iterations." @@ -122,8 +155,8 @@ 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. + -- Don't add another auto tune period if r == 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." @@ -137,23 +170,23 @@ return a' -- Auto tune the proposals. -mcmcAutotune :: Algorithm a => a -> MCMC a -mcmcAutotune a = do +mcmcAutotune :: Algorithm a => Int -> a -> MCMC a +mcmcAutotune n a = do logDebugB "Auto tune." - return $ aAutoTune a + liftIO $ aAutoTune n 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' + a'' <- mcmcAutotune x a' logInfoB $ aSummarizeCycle a'' logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations." mcmcResetAcceptance a'' -mcmcBurnInWithAutoTuning (x:xs) a = do +mcmcBurnInWithAutoTuning (x : xs) a = do a' <- mcmcIterate x a - a'' <- mcmcAutotune a' + a'' <- mcmcAutotune x a' logDebugB $ aSummarizeCycle a'' logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations." logDebugB $ aStdMonitorHeader a'' @@ -200,11 +233,9 @@ -- Initialize the run, execute the run, and close the run. mcmcRun :: Algorithm a => a -> MCMC a mcmcRun a = do + -- Header. logInfoHeader - - -- Debug settings. - logDebugB "The MCMC settings are:" - reader settings >>= logDebugS . ppShow + reader settings >>= logInfoB . settingsPrettyPrint -- Initialize. a' <- mcmcInitialize a 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 @@ -20,8 +20,7 @@ import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL -import Data.Time.Clock -import Data.Time.Format +import Data.Time import Mcmc.Internal.ByteString -- | Adapted from System.ProgressBar.renderDuration of package diff --git a/src/Mcmc/Posterior.hs b/src/Mcmc/Posterior.hs --- a/src/Mcmc/Posterior.hs +++ b/src/Mcmc/Posterior.hs @@ -11,10 +11,22 @@ -- Creation date: Fri May 28 12:26:35 2021. module Mcmc.Posterior ( Posterior, + PosteriorG, + PosteriorFunction, + PosteriorFunctionG, ) where import Numeric.Log -- | Posterior values are stored in log domain. -type Posterior = Log Double +type Posterior = PosteriorG Double + +-- | Generalized posterior. +type PosteriorG a = Log a + +-- | Posterior function. +type PosteriorFunction a = PosteriorFunctionG a Double + +-- | Generalized posterior function. +type PosteriorFunctionG a b = a -> PosteriorG b diff --git a/src/Mcmc/Prior.hs b/src/Mcmc/Prior.hs --- a/src/Mcmc/Prior.hs +++ b/src/Mcmc/Prior.hs @@ -13,13 +13,15 @@ -- Creation date: Thu Jul 23 13:26:14 2020. module Mcmc.Prior ( Prior, + PriorG, PriorFunction, - noPrior, + PriorFunctionG, -- * Improper priors - largerThan, + noPrior, + greaterThan, positive, - lowerThan, + lessThan, negative, -- * Continuous priors @@ -42,67 +44,88 @@ import Control.Monad import Data.Maybe (fromMaybe) +import Mcmc.Internal.Gamma import Mcmc.Statistics.Types import Numeric.Log import qualified Statistics.Distribution as S -import qualified Statistics.Distribution.Exponential as S -import qualified Statistics.Distribution.Gamma as S -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 +type Prior = PriorG Double +-- | Generalized prior. +type PriorG a = Log a + -- | Prior function. -type PriorFunction a = a -> Prior +type PriorFunction a = PriorFunctionG a Double +-- | Generalized prior function. +type PriorFunctionG a b = a -> PriorG b + -- | Flat prior function. Useful for testing and debugging. -noPrior :: PriorFunction a +noPrior :: RealFloat b => PriorFunctionG a b noPrior = const 1.0 +{-# SPECIALIZE noPrior :: PriorFunction Double #-} --- | Improper uniform prior; strictly larger than a given value. -largerThan :: LowerBoundary -> PriorFunction Double -largerThan a x - | x <= a = 0 - | otherwise = 1 +-- | Improper uniform prior; strictly greater than a given value. +greaterThan :: RealFloat a => LowerBoundary a -> PriorFunctionG a a +greaterThan a x + | x > a = 1.0 + | otherwise = 0.0 +{-# SPECIALIZE greaterThan :: Double -> PriorFunction Double #-} --- | Improper uniform prior; strictly larger than zero. -positive :: PriorFunction Double -positive = largerThan 0 +-- | Improper uniform prior; strictly greater than zero. +positive :: RealFloat a => PriorFunctionG a a +positive = greaterThan 0 +{-# SPECIALIZE positive :: PriorFunction Double #-} --- | Improper uniform prior; strictly lower than a given value. -lowerThan :: UpperBoundary -> PriorFunction Double -lowerThan b x - | x >= b = 0 - | otherwise = 1 +-- | Improper uniform prior; strictly less than a given value. +lessThan :: RealFloat a => UpperBoundary a -> PriorFunctionG a a +lessThan a x + | x < a = 1.0 + | otherwise = 0.0 +{-# SPECIALIZE lessThan :: Double -> PriorFunction Double #-} --- | Improper uniform prior; strictly lower than zero. -negative :: PriorFunction Double -negative = lowerThan 0 +-- | Improper uniform prior; strictly less than zero. +negative :: RealFloat a => PriorFunctionG a a +negative = lessThan 0 +{-# SPECIALIZE negative :: PriorFunction Double #-} -- | Exponential distributed prior. -exponential :: Rate -> PriorFunction Double -exponential l = Exp . S.logDensity d +-- +-- Call 'error' if the rate is zero or negative. +exponential :: RealFloat a => Rate a -> PriorFunctionG a a +exponential l x + | l <= 0 = error "exponential: Rate is zero or negative." + | x < 0 = error "exponential: Negative value." + | otherwise = ll * Exp (negate l * x) where - d = S.exponential l + ll = Exp $ log l +{-# SPECIALIZE exponential :: Double -> PriorFunction Double #-} -- | Gamma distributed prior. -gamma :: Shape -> Scale -> PriorFunction Double -gamma k t = Exp . S.logDensity d - where - d = S.gammaDistr k t +-- +-- Call 'error' if the shape or scale are zero or negative. +gamma :: RealFloat a => Shape a -> Scale a -> PriorFunctionG a a +gamma k t x + | k <= 0 = error "gamma: Shape is zero or negative." + | t <= 0 = error "gamma: Scale is zero or negative." + | x < 0 = error "gamma: Negative value." + | x == 0 = 0.0 + | otherwise = Exp $ log x * (k - 1) - (x / t) - logGammaG k - log t * k +{-# SPECIALIZE gamma :: Double -> Double -> PriorFunction Double #-} -- | 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 +gammaMeanVariance :: RealFloat a => Mean a -> Variance a -> PriorFunctionG a a +gammaMeanVariance m v = gamma k t + where + (k, t) = gammaMeanVarianceToShapeScale m v +{-# SPECIALIZE gammaMeanVariance :: Double -> Double -> PriorFunction Double #-} -- | 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) +gammaMeanOne :: RealFloat a => Shape a -> PriorFunctionG a a +gammaMeanOne k = gamma k (recip k) +{-# SPECIALIZE gammaMeanOne :: Double -> PriorFunction Double #-} -- The mean and variance of the gamma distribution are -- @@ -118,30 +141,51 @@ -- | Calculate mean and variance of the gamma distribution given the shape and -- the scale. -gammaShapeScaleToMeanVariance :: Shape -> Scale -> (Mean, Variance) +gammaShapeScaleToMeanVariance :: Num a => Shape a -> Scale a -> (Mean a, Variance a) gammaShapeScaleToMeanVariance k t = let m = k * t in (m, m * t) +{-# SPECIALIZE gammaShapeScaleToMeanVariance :: Double -> Double -> (Double, Double) #-} -- | Calculate shape and scale of the gamma distribution given the mean and -- the variance. -gammaMeanVarianceToShapeScale :: Mean -> Variance -> (Shape, Scale) +gammaMeanVarianceToShapeScale :: Fractional a => Mean a -> Variance a -> (Shape a, Scale a) gammaMeanVarianceToShapeScale m v = (m * m / v, v / m) +{-# SPECIALIZE gammaMeanVarianceToShapeScale :: Double -> Double -> (Double, Double) #-} +mLnSqrt2Pi :: RealFloat a => a +mLnSqrt2Pi = 0.9189385332046727417803297364056176398613974736377834128171 +{-# INLINE mLnSqrt2Pi #-} + -- | Normal distributed prior. -normal :: Mean -> StandardDeviation -> PriorFunction Double -normal m s = Exp . S.logDensity d +-- +-- Call 'error' if the standard deviation is zero or negative. +normal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a +normal m s x + | s <= 0 = error "normal: Standard deviation is zero or negative." + | otherwise = Exp $ (- xm * xm / (2 * s * s)) - denom where - d = S.normalDistr m s + xm = x - m + denom = mLnSqrt2Pi + log s +{-# SPECIALIZE normal :: Double -> Double -> PriorFunction Double #-} -- | Uniform prior on [a, b]. -uniform :: LowerBoundary -> UpperBoundary -> PriorFunction Double +-- +-- Call 'error' if the lower boundary is greather than the upper boundary. +uniform :: RealFloat a => LowerBoundary a -> UpperBoundary a -> PriorFunctionG a a uniform a b x - | x <= a = 0 - | x >= b = 0 - | otherwise = Exp 0 + | a > b = error "uniform: Lower boundary is greater than upper boundary." + | x < a = 0.0 + | x > b = 0.0 + | otherwise = 1.0 +{-# SPECIALIZE uniform :: Double -> Double -> PriorFunction Double #-} -- | Poisson distributed prior. -poisson :: Rate -> PriorFunction Int -poisson l = Exp . S.logProbability d +-- +-- Call 'error' if the rate is zero or negative. +poisson :: Rate Double -> PriorFunction Int +poisson l n + | l <= 0 = error "poisson: Rate is zero or negative." + | n < 0 = error "poisson: Negative value." + | otherwise = Exp $ S.logProbability d n where d = S.poisson l @@ -149,9 +193,11 @@ -- -- Use with care because the elements are checked for positiveness, and this can -- take some time if the list is long and does not contain any zeroes. -product' :: [Log Double] -> Log Double +product' :: RealFloat a => [Log a] -> Log a product' = fromMaybe 0 . prodM +{-# SPECIALIZE product' :: [Log Double] -> Log Double #-} -- The type could be generalized to any MonadPlus Integer -prodM :: [Log Double] -> Maybe (Log Double) +prodM :: RealFloat a => [Log a] -> Maybe (Log a) prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0)) 1 +{-# SPECIALIZE prodM :: [Log Double] -> Maybe (Log Double) #-} diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs --- a/src/Mcmc/Proposal.hs +++ b/src/Mcmc/Proposal.hs @@ -1,4 +1,3 @@ -{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DerivingVia #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE RankNTypes #-} @@ -31,53 +30,31 @@ ProposalSimple, Tuner (..), Tune (..), + defaultTuningFunction, createProposal, TuningParameter, + AuxiliaryTuningParameters, tuningParameterMin, tuningParameterMax, - tune, + tuneWithTuningParameters, + tuneWithChainParameters, getOptimalRate, - -- * Cycles - Order (..), - Cycle (ccProposals), - cycleFromList, - setOrder, - prepareProposals, - tuneCycle, - autoTuneCycle, - - -- * Acceptance rates - Acceptance (fromAcceptance), - emptyA, - pushA, - resetA, - transformKeysA, - acceptanceRate, - acceptanceRates, - -- * Output proposalHeader, - proposalHLine, summarizeProposal, - summarizeCycle, ) where -import Control.DeepSeq -import Data.Aeson -import Data.Bifunctor import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL -import Data.Default import qualified Data.Double.Conversion.ByteString as BC import Data.Function -import Data.List -import qualified Data.Map.Strict as M -import Data.Maybe +import qualified Data.Vector as VB import Lens.Micro +import Lens.Micro.Extras +import Mcmc.Acceptance import Mcmc.Internal.ByteString -import Mcmc.Internal.Shuffle import Numeric.Log hiding (sum) import System.Random.MWC @@ -97,8 +74,9 @@ -- | Check if the weight is positive. pWeight :: Int -> PWeight -pWeight n | n <= 0 = error "pWeight: Proposal weight is zero or negative." - | otherwise = PWeight n +pWeight n + | n <= 0 = error "pWeight: Proposal weight is zero or negative." + | otherwise = PWeight n -- | Proposal dimension. -- @@ -108,8 +86,8 @@ -- high dimensional ones. -- -- Optimal acceptance rates are still subject to controversies. As far as I --- know, research has focused on random walk proposal with a multivariate normal --- distribution of dimension @d@. In this case, the following acceptance rates +-- know, research has focused on random walk proposals with multivariate normal +-- distributions of dimension @d@. In this case, the following acceptance rates -- are desired: -- -- - one dimension: 0.44 (numerical results); @@ -123,19 +101,27 @@ -- Of course, many proposals may not be classical random walk proposals. For -- example, the beta proposal on a simplex ('Mcmc.Proposal.Simplex.beta') -- samples one new variable of the simplex from a beta distribution while --- rescaling all other variables. What is the dimension of this proposal? I --- don't know, but I set the dimension to 2. The reason is that if the dimension --- of the simplex is 2, two variables are changed. If the dimension of the --- simplex is high, one variable is changed substantially, while all others are --- changed marginally. +-- rescaling all other variables. What is the dimension of this proposal? Here, +-- the dimension is set to 2. The reason is that if the dimension of the simplex +-- is 2, two variables are changed. If the dimension of the simplex is high, one +-- variable is changed substantially, while all others are changed marginally. -- -- Further, if a proposal changes a number of variables in the same way (and not --- independently like in a random walk proposal), I still set the dimension of --- the proposal to the number of variables changed. +-- independently like in a random walk proposal), the dimension of the proposal +-- is set to the number of variables changed. -- --- Finally, I assume that proposals of unknown dimension have high dimension, --- and use the optimal acceptance rate 0.234. -data PDimension = PDimension Int | PDimensionUnknown +-- Moreover, proposals of unknown dimension are assumed to have high dimension, +-- and the optimal acceptance rate 0.234 is used. +-- +-- Finally, special proposals may have completely different desired acceptance +-- rates. For example. the Hamiltonian Monte Carlo proposal (see +-- Mcmc.Proposal.Hamiltonian.hmc) has a desired acceptance rate of 0.65. +-- Specific acceptance rates can be set with 'PSpecial'. +data PDimension + = PDimension Int + | PDimensionUnknown + | -- | Provide dimension ('Int') and desired acceptance rate ('Double'). + PSpecial Int Double -- | A 'Proposal' is an instruction about how the Markov chain will traverse the -- state space @a@. Essentially, it is a probability mass or probability density @@ -192,20 +178,21 @@ -- | Lift a proposal from one data type to another. -- --- Assume the Jacobian is 1.0 (see also 'liftProposal' and 'liftProposalWith'). +-- Assume the Jacobian is 1.0. -- -- For example: -- -- @ -- scaleFirstEntryOfTuple = _1 @~ scale -- @ +-- +-- See also 'liftProposal' and 'liftProposalWith'. infixl 7 @~ + (@~) :: Lens' b a -> Proposal a -> Proposal b (@~) = liftProposal --- | Lift a proposal from one data type to another. --- --- Assume the Jacobian is 1.0 (see also '(@~)' and 'liftProposalWith'). +-- | See '(@~)'. liftProposal :: Lens' b a -> Proposal a -> Proposal b liftProposal = liftProposalWith (const 1.0) @@ -217,7 +204,8 @@ -- 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) +liftProposalWith jf l (Proposal n r d w s t) = + Proposal n r d w (liftProposalSimple jf l s) (liftTuner jf l <$> t) -- | Simple proposal without tuning information. -- @@ -239,8 +227,9 @@ -- determinant of the Jacobian matrix differs from 1.0. type ProposalSimple a = a -> GenIO -> IO (a, KernelRatio, Jacobian) -convertProposalSimple :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b -convertProposalSimple jf l s = s' +-- Lift a simple proposal from one data type to another. +liftProposalSimple :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b +liftProposalSimple jf l s = s' where s' y g = do (x', r, j) <- s (y ^. l) g @@ -250,31 +239,72 @@ j' = j * jyx / jxy return (y', r, j') --- | Tune the acceptance rate of a 'Proposal'; see 'tune', or 'autoTuneCycle'. +-- | Required information to tune 'Proposal's. data Tuner a = Tuner - { tParam :: TuningParameter, - tFunc :: TuningParameter -> ProposalSimple a + { tGetTuningParameter :: TuningParameter, + -- | Instruction about how to compute new tuning parameter from a given + -- acceptance rate and the old tuning parameter. + tComputeTuningParameter :: AcceptanceRate -> TuningParameter -> TuningParameter, + tGetAuxiliaryTuningParameters :: AuxiliaryTuningParameters, + -- | Instruction about how to compute new auxiliary tuning parameters from a + -- given trace and the old auxiliary tuning parameters. + tComputeAuxiliaryTuningParameters :: + VB.Vector a -> + AuxiliaryTuningParameters -> + AuxiliaryTuningParameters, + -- | Given the tuning parameter, and the auxiliary tuning parameters, get + -- the tuned simple proposal. + -- + -- Should return 'Left' if the vector of auxiliary tuning parameters is + -- invalid. + tGetSimpleProposal :: + TuningParameter -> + AuxiliaryTuningParameters -> + Either String (ProposalSimple a) } -convertTuner :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b -convertTuner jf l (Tuner p f) = Tuner p f' +-- Lift tuner from one data type to another. +liftTuner :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b +liftTuner jf l (Tuner p fP ps fPs g) = Tuner p fP ps fPs' g' where - f' x = convertProposalSimple jf l $ f x + fPs' = fPs . VB.map (view l) + g' x xs = liftProposalSimple jf l <$> g x xs --- | Tune the proposal? +-- | Tune proposal? data Tune = Tune | NoTune deriving (Show, Eq) -- | Tuning parameter. +-- +-- The larger the tuning parameter, the larger the proposal and the lower the +-- expected acceptance rate; and vice versa. type TuningParameter = Double --- | Create a tuneable proposal. +-- | Auxiliary tuning parameters; vector may be empty. +type AuxiliaryTuningParameters = VB.Vector TuningParameter + +-- | Default tuning function. +-- +-- Subject to change. +defaultTuningFunction :: + -- Optimal acceptance rate. + PDimension -> + AcceptanceRate -> + TuningParameter -> + TuningParameter +defaultTuningFunction d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t + +noAuxiliaryTuningFunction :: VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters +noAuxiliaryTuningFunction _ ts = ts + +-- | Create a proposal with a single tuning parameter. +-- +-- Proposals with arbitrary tuning parameters have to be created manually. See +-- 'Tuner' for more information, and 'Mcmc.Proposal.Hamiltonian' for an example. createProposal :: -- | Description of the proposal type and parameters. PDescription -> - -- | Function creating a simple proposal for a given tuning parameter. The - -- larger the tuning parameter, the larger the proposal and the lower the - -- expected acceptance rate; and vice versa. + -- | Function creating a simple proposal for a given tuning parameter. (TuningParameter -> ProposalSimple a) -> -- | Dimension. PDimension -> @@ -285,17 +315,30 @@ -- | Activate tuning? Tune -> Proposal a -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 +createProposal r f d n w Tune = + Proposal n r d w (f 1.0) (Just tuner) + where + fT = defaultTuningFunction d + fTs = noAuxiliaryTuningFunction + g t _ = Right $ f t + tuner = Tuner 1.0 fT VB.empty fTs g +createProposal r f d n w NoTune = + Proposal n r d w (f 1.0) Nothing --- | Minimal tuning parameter; @1e-12@, subject to change. +-- IDEA: Per proposal type tuning parameter boundaries. For example, a sliding +-- proposal with a large tuning parameter is not a problem. But then, if the +-- tuning parameters are very different from one, a different base proposal +-- should be chosen. + +-- | Minimal tuning parameter; @1e-5@, subject to change. -- -- >>> tuningParameterMin -- 1e-5 tuningParameterMin :: TuningParameter tuningParameterMin = 1e-5 --- | Maximal tuning parameter; @1e12@, subject to change. +-- | Maximal tuning parameter; @1e3@, subject to change. +-- -- >>> tuningParameterMax -- 1e3 tuningParameterMax :: TuningParameter @@ -303,20 +346,42 @@ -- | Tune a 'Proposal'. -- --- The size of the proposal is proportional to the tuning parameter which has a --- positive lower bound of 'tuningParameterMin'. +-- The size of the proposal is proportional to the tuning parameter which has +-- positive lower and upper boundaries of 'tuningParameterMin' and +-- 'tuningParameterMax', respectively. -- --- The tuning function maps the current tuning parameter to a new one. +-- Auxiliary tuning parameters may also be used by the 'Tuner' of the proposal. -- --- Return 'Nothing' if 'Proposal' is not tuneable. -tune :: (TuningParameter -> TuningParameter) -> Proposal a -> Maybe (Proposal a) -tune f m = do - (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} +-- Return 'Left' if: +-- +-- - the 'Proposal' is not tuneable; +-- +-- - the auxiliary tuning parameters are invalid. +tuneWithTuningParameters :: + TuningParameter -> + AuxiliaryTuningParameters -> + Proposal a -> + Either String (Proposal a) +tuneWithTuningParameters t ts p = case prTuner p of + Nothing -> Left "tuneWithTuningParameters: Proposal is not tunable." + Just (Tuner _ fT _ fTs g) -> + -- Ensure that the tuning parameter is strictly positive and well bounded. + let t' = max tuningParameterMin t + t'' = min tuningParameterMax t' + psE = g t'' ts + in case psE of + Left err -> Left $ "tune: " <> err + Right ps -> Right $ p {prSimple = ps, prTuner = Just $ Tuner t'' fT ts fTs g} +-- | See 'tuneWithTuningParameters' and 'Tuner'. +tuneWithChainParameters :: AcceptanceRate -> VB.Vector a -> Proposal a -> Either String (Proposal a) +tuneWithChainParameters ar xs p = case prTuner p of + Nothing -> Left "tuneWithChainParameters: Proposal is not tunable." + Just (Tuner t fT ts fTs _) -> + let t' = fT ar t + ts' = fTs xs ts + in tuneWithTuningParameters t' ts' p + -- | See 'PDimension'. getOptimalRate :: PDimension -> Double getOptimalRate (PDimension n) @@ -329,6 +394,7 @@ | n >= 5 = 0.234 | otherwise = error "getOptimalRate: Proposal dimension is not an integer?" getOptimalRate PDimensionUnknown = 0.234 +getOptimalRate (PSpecial _ r) = r -- Warn if acceptance rate is lower. rateMin :: Double @@ -338,122 +404,6 @@ rateMax :: Double rateMax = 0.9 --- | Define the order in which 'Proposal's are executed in a 'Cycle'. The total --- number of 'Proposal's per 'Cycle' may differ between 'Order's (e.g., compare --- 'RandomO' and 'RandomReversibleO'). -data Order - = -- | Shuffle the 'Proposal's in the 'Cycle'. The 'Proposal's are replicated - -- according to their weights and executed in random order. If a 'Proposal' has - -- weight @w@, it is executed exactly @w@ times per iteration. - RandomO - | -- | The 'Proposal's are executed sequentially, in the order they appear in the - -- 'Cycle'. 'Proposal's with weight @w>1@ are repeated immediately @w@ times - -- (and not appended to the end of the list). - SequentialO - | -- | Similar to 'RandomO'. However, a reversed copy of the list of - -- shuffled 'Proposal's is appended such that the resulting Markov chain is - -- reversible. - -- Note: the total number of 'Proposal's executed per cycle is twice the number - -- of 'RandomO'. - RandomReversibleO - | -- | Similar to 'SequentialO'. However, a reversed copy of the list of - -- sequentially ordered 'Proposal's is appended such that the resulting Markov - -- chain is reversible. - SequentialReversibleO - deriving (Eq, Show) - -instance Default Order where def = RandomO - --- Describe the order. -describeOrder :: Order -> BL.ByteString -describeOrder RandomO = "The proposals are executed in random order." -describeOrder SequentialO = "The proposals are executed sequentially." -describeOrder RandomReversibleO = - BL.intercalate - "\n" - [ describeOrder RandomO, - "A reversed copy of the shuffled proposals is appended to ensure reversibility." - ] -describeOrder SequentialReversibleO = - BL.intercalate - "\n" - [ describeOrder SequentialO, - "A reversed copy of the sequential proposals is appended to ensure reversibility." - ] - --- | In brief, a 'Cycle' is a list of proposals. --- --- The state of the Markov chain will be logged only after all 'Proposal's in --- the 'Cycle' have been completed, and the iteration counter will be increased --- by one. The order in which the 'Proposal's are executed is specified by --- 'Order'. The default is 'RandomO'. --- --- No proposals with the same name and description are allowed in a 'Cycle', so --- that they can be uniquely identified. -data Cycle a = Cycle - { ccProposals :: [Proposal a], - ccOrder :: Order - } - --- | Create a 'Cycle' from a list of 'Proposal's. -cycleFromList :: [Proposal a] -> Cycle a -cycleFromList [] = - error "cycleFromList: Received an empty list but cannot create an empty Cycle." -cycleFromList xs = - if length (nub xs) == length xs - then Cycle xs def - else error "cycleFromList: Proposals are not unique." - --- | Set the order of 'Proposal's in a 'Cycle'. -setOrder :: Order -> Cycle a -> Cycle a -setOrder o c = c {ccOrder = o} - --- | Replicate 'Proposal's according to their weights and possibly shuffle them. -prepareProposals :: Cycle a -> GenIO -> IO [Proposal a] -prepareProposals (Cycle xs o) g = case o of - RandomO -> shuffle ps g - SequentialO -> return ps - RandomReversibleO -> do - psR <- shuffle ps g - return $ psR ++ reverse psR - SequentialReversibleO -> return $ ps ++ reverse ps - where - !ps = concat [replicate (fromPWeight $ prWeight p) p | p <- xs] - --- The number of proposals depends on the order. -getNProposalsPerCycle :: Cycle a -> Int -getNProposalsPerCycle (Cycle xs o) = case o of - RandomO -> once - SequentialO -> once - RandomReversibleO -> 2 * once - SequentialReversibleO -> 2 * once - where - once = sum $ map (fromPWeight . prWeight) xs - --- | Tune 'Proposal's in the 'Cycle'. See 'tune'. -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} - else error "tuneCycle: Propoals in map and cycle do not match." - where - ps = ccProposals c - tuneF p = case m M.!? p of - Nothing -> p - Just f -> fromMaybe p (tune f p) - --- | Calculate acceptance rates and auto tune the 'Proposal's in the 'Cycle'. For --- now, a 'Proposal' is enlarged when the acceptance rate is above 0.44, and --- shrunk otherwise. Do not change 'Proposal's that are not tuneable. -autoTuneCycle :: Acceptance (Proposal a) -> Cycle a -> Cycle a -autoTuneCycle a = tuneCycle (M.mapWithKey tuningF $ acceptanceRates a) - where - 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 -> BL.ByteString -> @@ -477,67 +427,6 @@ 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 = @@ -552,10 +441,6 @@ "Tuning parameter" "Consider manual adjustment" --- | Horizontal line of proposal summaries. -proposalHLine :: BL.ByteString -proposalHLine = BL.replicate (BL.length proposalHeader) '-' - -- | Proposal summary. summarizeProposal :: PName -> @@ -582,7 +467,7 @@ 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) tuningParameter + tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 4) tuningParameter checkRate rate | rate < rateMin = Just "rate too low" | rate > rateMax = Just "rate too high" @@ -599,32 +484,3 @@ (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 -summarizeCycle a c = - BL.intercalate "\n" $ - [ "Summary of proposal(s) in cycle.", - nProposalsFullStr, - describeOrder (ccOrder c), - proposalHeader, - proposalHLine - ] - ++ [ summarizeProposal - (prName p) - (prDescription p) - (prWeight p) - (tParam <$> prTuner p) - (prDimension p) - (ar p) - | p <- ps - ] - ++ [proposalHLine] - where - ps = ccProposals c - nProposals = getNProposalsPerCycle c - nProposalsStr = BB.toLazyByteString $ BB.intDec nProposals - nProposalsFullStr = case nProposals of - 1 -> nProposalsStr <> " proposal is performed per iteration." - _ -> nProposalsStr <> " proposals are performed per iterations." - ar m = acceptanceRate m a 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 @@ -31,7 +31,7 @@ genBactrian :: SpikeParameter -> - StandardDeviation -> + StandardDeviation Double -> GenIO -> IO Double genBactrian m s g = do @@ -42,7 +42,7 @@ b <- bernoulli 0.5 g return $ if b then x else - x -logDensityBactrian :: SpikeParameter -> StandardDeviation -> Double -> Log Double +logDensityBactrian :: SpikeParameter -> StandardDeviation Double -> Double -> Log Double logDensityBactrian m s x = Exp $ log $ kernel1 + kernel2 where mn = m * s @@ -54,7 +54,7 @@ bactrianAdditive :: SpikeParameter -> - StandardDeviation -> + StandardDeviation Double -> ProposalSimple Double bactrianAdditive m s x g = do dx <- genBactrian m s g @@ -63,8 +63,8 @@ -- bactrianSimple lens spike stdDev tune forwardOp backwardOp bactrianAdditiveSimple :: SpikeParameter -> - StandardDeviation -> - TuningParameter -> + StandardDeviation Double -> + TuningParameter -> ProposalSimple Double bactrianAdditiveSimple m s t | m < 0 = error "bactrianAdditiveSimple: Spike parameter negative." @@ -85,7 +85,7 @@ -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. slideBactrian :: SpikeParameter -> - StandardDeviation -> + StandardDeviation Double -> PName -> PWeight -> Tune -> @@ -105,7 +105,7 @@ bactrianMult :: SpikeParameter -> - StandardDeviation -> + StandardDeviation Double -> ProposalSimple Double bactrianMult m s x g = do du <- genBactrian m s g @@ -115,7 +115,11 @@ jac = Exp $ log $ recip u return (x * u, qYX / qXY, jac) -bactrianMultSimple :: SpikeParameter -> StandardDeviation -> TuningParameter -> ProposalSimple Double +bactrianMultSimple :: + SpikeParameter -> + StandardDeviation Double -> + TuningParameter -> + ProposalSimple Double bactrianMultSimple m s t | m < 0 = error "bactrianMultSimple: Spike parameter negative." | m >= 1 = error "bactrianMultSimple: Spike parameter 1.0 or larger." @@ -123,10 +127,12 @@ | otherwise = bactrianMult m (t * s) -- | Multiplicative proposal with kernel similar to the silhouette of a Bactrian --- camel. See 'slideBactrian'. +-- camel. +-- +-- See 'Mcmc.Proposal.Scale.scale', and 'slideBactrian'. scaleBactrian :: SpikeParameter -> - StandardDeviation -> + StandardDeviation Double -> PName -> PWeight -> Tune -> 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 @@ -19,17 +19,36 @@ import Numeric.Log import Statistics.Distribution --- | Generic function to create proposals for continuous parameters ('Double'). +-- | Generic function to create proposals for continuous parameters (e.g., +-- 'Double'). +-- +-- The procedure is as follows: Let \(\mathbb{X}\) be the state space and \(x\) +-- be the current state. +-- +-- 1. Let \(D\) be a continuous probability distribution on \(\mathbb{D}\); +-- sample an auxiliary variable \(epsilon \sim D\). +-- +-- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). Propose a +-- new state \(x' = x \odot \epsilon\). +-- +-- 3. If the proposal is unbiased, the Metropolis-Hastings-Green ratio can +-- directly be calculated using the posterior function. +-- +-- 4. However, if the proposal is biased: Suppose \(g : \mathbb{D} \to +-- \mathbb{D}\) inverses the auxiliary variable \(\epsilon\) such that \(x = +-- x' \odot g(\epsilon)\). Calculate the Metropolis-Hastings-Green ratio +-- using the posterior function, \(g\), \(D\), \(\epsilon\), and possibly a +-- Jacobian function. genericContinuous :: (ContDistr d, ContGen d) => -- | Probability distribution d -> - -- | Forward operator. + -- | Forward operator \(\odot\). -- -- For example, for a multiplicative proposal on one variable the forward -- operator is @(*)@, so that @x * u = y@. (a -> Double -> a) -> - -- | Inverse operator. + -- | Inverse operator \(g\) of the auxiliary variable. -- -- For example, 'recip' for a multiplicative proposal on one variable, since -- @y * (recip u) = x * u * (recip u) = x@. @@ -46,8 +65,8 @@ -- That is, the determinant of the Jacobian matrix of multiplication is just -- the reciprocal value of @u@ (with conversion to log domain). -- - -- Required for proposals for which absolute value of the determinant of the - -- Jacobian differs from 1.0. + -- Required for proposals for which the absolute value of the determinant of + -- the Jacobian differs from 1.0. -- -- Conversion to log domain is necessary, because some determinants of -- Jacobians are very small (or large). @@ -67,15 +86,22 @@ return (x `f` u, r, j) {-# INLINEABLE genericContinuous #-} --- | Generic function to create proposals for discrete parameters ('Int'). +-- | Generic function to create proposals for discrete parameters (e.g., 'Int'). +-- +-- See 'genericContinuous'. genericDiscrete :: (DiscreteDistr d, DiscreteGen d) => -- | Probability distribution. d -> - -- | Forward operator, e.g. (+), so that x + dx = x'. + -- | Forward operator. + -- + -- For example, (+), so that x + dx = x'. (a -> Int -> a) -> - -- | Inverse operator, e.g., 'negate', so that x' + (negate dx) = x. Only - -- required for biased proposals. + -- | Inverse operator \(g\) of the auxiliary variable. + -- + -- For example, 'negate', so that x' + (negate dx) = x. + -- + -- Only required for biased proposals. Maybe (Int -> Int) -> ProposalSimple a genericDiscrete d f mfInv x g = do diff --git a/src/Mcmc/Proposal/Hamiltonian.hs b/src/Mcmc/Proposal/Hamiltonian.hs new file mode 100644 --- /dev/null +++ b/src/Mcmc/Proposal/Hamiltonian.hs @@ -0,0 +1,448 @@ +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE UndecidableInstances #-} + +-- | +-- Module : Mcmc.Proposal.Hamiltonian +-- Description : Hamiltonian Monte Carlo proposal +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Mon Jul 5 12:59:42 2021. +-- +-- The Hamiltonian Monte Carlo (HMC) proposal. +-- +-- For references, see: +-- +-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using +-- Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.), +-- Handbook of Markov Chain Monte Carlo (2011), CRC press. +-- +-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data +-- analysis (2014), CRC Press. +-- +-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual +-- introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017). +-- +-- NOTE on implementation: +-- +-- - The implementation assumes the existence of the gradient. Like so, the user +-- can use automatic or manual differentiation, depending on the problem at +-- hand. +-- +-- - The state needs to be list like or 'Traversable' so that the structure of +-- the state space is available. A 'Traversable' constraint on the data type +-- is nice because it is more general than, for example, a list, and +-- user-defined data structures can be used. +-- +-- - The state needs to have a zip-like 'Applicative' instance so that +-- - matrix/vector operations can be performed. + +module Mcmc.Proposal.Hamiltonian + ( Gradient, + Masses, + LeapfrogTrajectoryLength, + LeapfrogScalingFactor, + HTune (..), + HSettings (..), + hamiltonian, + ) +where + +import Data.Foldable +import qualified Data.Matrix as M +import Data.Maybe +import Data.Traversable +import qualified Data.Vector as VB +import Mcmc.Prior +import Mcmc.Proposal +import Numeric.Log +import Statistics.Distribution +import Statistics.Distribution.Normal +import qualified Statistics.Function as S +import qualified Statistics.Sample as S +import System.Random.MWC + +-- TODO: At the moment, the HMC proposal is agnostic of the prior and +-- likelihood, that is, the posterior function. This means, that it cannot know +-- when it reaches a point with zero posterior probability. This also affects +-- restricted or constrained parameters. See Gelman p. 303. + +-- TODO: No-U-turn sampler. + +-- TODO: Riemannian adaptation. + +-- | Gradient of the log posterior function. +type Gradient f = f Double -> f Double + +-- | Function validating the state. +-- +-- Useful if parameters are constrained. +type Validate f = f Double -> Bool + +-- | Masses of parameters. +-- +-- NOTE: Full specification of a mass matrix including off-diagonal elements is +-- not supported. +-- +-- NOTE: Parameters without masses ('Nothing') are not changed by the +-- Hamiltonian proposal. +-- +-- The masses roughly describe how reluctant the particle moves through the +-- state space. If a parameter has higher mass, the momentum in this direction +-- will be changed less by the provided gradient, than when the same parameter +-- has lower mass. +-- +-- The proposal is more efficient if masses are assigned according to the +-- inverse (co)-variance structure of the posterior function. That is, +-- parameters changing on larger scales should have lower masses than parameters +-- changing on lower scales. In particular, and for a diagonal mass matrix, the +-- optimal masses are the inverted variances of the parameters distributed +-- according to the posterior function. +-- +-- Of course, the scales of the parameters of the posterior function are usually +-- unknown. Often, it is sufficient to +-- +-- - set the masses to identical values roughly scaled with the inverted +-- estimated average variance of the posterior function; or even to +-- +-- - set all masses to 1.0, and trust the tuning algorithm (see +-- 'HTuneMassesAndLeapfrog') to find the correct values. +type Masses f = f (Maybe Double) + +-- | Mean leapfrog trajectory length \(L\). +-- +-- Number of leapfrog steps per proposal. +-- +-- To avoid problems with ergodicity, the actual number of leapfrog steps is +-- sampled proposal from a discrete uniform distribution over the interval +-- \([\text{floor}(0.8L),\text{ceiling}(1.2L)]\). +-- +-- For a discussion of ergodicity and reasons why randomization is important, +-- see [1] p. 15; also mentioned in [2] p. 304. +-- +-- NOTE: To avoid errors, the left bound has an additional hard minimum of 1, +-- and the right bound is required to be larger equal than the left bound. +-- +-- Usually set to 10, but larger values may be desirable. +type LeapfrogTrajectoryLength = Int + +-- | Mean of leapfrog scaling factor \(\epsilon\). +-- +-- Determines the size of each leapfrog step. +-- +-- To avoid problems with ergodicity, the actual leapfrog scaling factor is +-- sampled per proposal from a continuous uniform distribution over the interval +-- \((0.8\epsilon,1.2\epsilon]\). +-- +-- For a discussion of ergodicity and reasons why randomization is important, +-- see [1] p. 15; also mentioned in [2] p. 304. +-- +-- Usually set such that \( L \epsilon = 1.0 \), but smaller values may be +-- required if acceptance rates are low. +type LeapfrogScalingFactor = Double + +-- Target state containing parameters. +type Positions f = f Double + +-- Momenta of the parameters. +type Momenta f = f (Maybe Double) + +-- | Tuning settings. +-- +-- +-- Tuning of leapfrog parameters: +-- +-- We expect that the larger the leapfrog step size the larger the proposal step +-- size and the lower the acceptance ratio. Consequently, if the acceptance rate +-- is too low, the leapfrog step size is decreased and vice versa. Further, the +-- leapfrog trajectory length is scaled such that the product of the leapfrog +-- step size and trajectory length stays constant. +-- +-- Tuning of masses: +-- +-- The variances of all parameters of the posterior distribution obtained over +-- the last auto tuning interval is calculated and the masses are amended using +-- the old masses and the inverted variances. If, for a specific coordinate, the +-- sample size is too low, or if the calculated variance is out of predefined +-- bounds, the mass of the affected position is not changed. +data HTune + = -- | Tune masses and leapfrog parameters. + HTuneMassesAndLeapfrog + | -- | Tune leapfrog parameters only. + HTuneLeapfrogOnly + | -- | Do not tune at all. + HNoTune + deriving (Eq, Show) + +-- | Specifications for Hamilton Monte Carlo proposal. +data HSettings f = HSettings + { hGradient :: Gradient f, + hMaybeValidate :: Maybe (Validate f), + hMasses :: Masses f, + hLeapfrogTrajectoryLength :: LeapfrogTrajectoryLength, + hLeapfrogScalingFactor :: LeapfrogScalingFactor, + hTune :: HTune + } + +checkHSettings :: Foldable f => HSettings f -> Maybe String +checkHSettings (HSettings _ _ masses l eps _) + | any f masses = Just "checkHSettings: One or more masses are zero or negative." + | l < 1 = Just "checkHSettings: Leapfrog trajectory length is zero or negative." + | eps <= 0 = Just "checkHSettings: Leapfrog scaling factor is zero or negative." + | otherwise = Nothing + where + f (Just m) = m <= 0 + f Nothing = False + +generateMomenta :: + Traversable f => + Masses f -> + GenIO -> + IO (Momenta f) +generateMomenta masses gen = traverse (generateWith gen) masses + where + generateWith g (Just m) = let d = normalDistr 0 (sqrt m) in Just <$> genContVar d g + generateWith _ Nothing = pure Nothing + +priorMomenta :: + (Applicative f, Foldable f) => + Masses f -> + Momenta f -> + Prior +priorMomenta masses phi = foldl' (*) 1.0 $ f <$> masses <*> phi + where + f (Just m) (Just p) = let d = normalDistr 0 (sqrt m) in Exp $ logDensity d p + f Nothing Nothing = 1.0 + f _ _ = error "priorMomenta: Got just a mass and no momentum or the other way around." + +leapfrog :: + Applicative f => + Gradient f -> + Maybe (Validate f) -> + Masses f -> + LeapfrogTrajectoryLength -> + LeapfrogScalingFactor -> + Positions f -> + Momenta f -> + -- Maybe (Positions', Momenta'). + Maybe (Positions f, Momenta f) +leapfrog grad mVal masses l eps theta phi = do + let -- The first half step of the momenta. + phiHalf = leapfrogStepMomenta 0.5 eps grad theta phi + -- L-1 full steps. This gives the positions theta_{L-1}, and the momenta + -- phi_{L-1/2}. + (thetaLM1, phiLM1Half) <- go (l - 1) (Just (theta, phiHalf)) + -- The last full step of the positions. + thetaL <- valF $ leapfrogStepPositions eps masses thetaLM1 phiLM1Half + let -- The last half step of the momenta. + phiL = leapfrogStepMomenta 0.5 eps grad thetaL phiLM1Half + return (thetaL, phiL) + where + valF x = case mVal of + Nothing -> Just x + Just f -> if f x then Just x else Nothing + go _ Nothing = Nothing + go 0 (Just (t, p)) = Just (t, p) + go n (Just (t, p)) = + let t' = leapfrogStepPositions eps masses t p + p' = leapfrogStepMomenta 1.0 eps grad t' p + r = (,p') <$> valF t' + in go (n - 1) r + +leapfrogStepMomenta :: + Applicative f => + -- Size of step (half or full step). + Double -> + LeapfrogScalingFactor -> + Gradient f -> + -- Current positions. + Positions f -> + -- Current momenta. + Momenta f -> + -- New momenta. + Momenta f +leapfrogStepMomenta xi eps grad theta phi = phi <+. ((xi * eps) .* grad theta) + where + (<+.) :: Applicative f => f (Maybe Double) -> f Double -> f (Maybe Double) + (<+.) xs ys = f <$> xs <*> ys + f Nothing _ = Nothing + f (Just x) y = Just $ x + y + +leapfrogStepPositions :: + Applicative f => + LeapfrogScalingFactor -> + Masses f -> + -- Current positions. + Positions f -> + -- Current momenta. + Momenta f -> + Positions f +-- The arguments are flipped to encounter the maybe momentum. +leapfrogStepPositions eps masses theta phi = theta <+. (mScaledReversed .*> phi) + where + (<+.) :: Applicative f => f Double -> f (Maybe Double) -> f Double + (<+.) xs ys = f <$> xs <*> ys + f x Nothing = x + f x (Just y) = x + y + mScaledReversed = (fmap . fmap) ((* eps) . (** (-1))) masses + (.*>) :: Applicative f => f (Maybe Double) -> f (Maybe Double) -> f (Maybe Double) + (.*>) xs ys = g <$> xs <*> ys + g (Just x) (Just y) = Just $ x * y + g Nothing Nothing = Nothing + g _ _ = error "leapfrogStepPositions: Got just a mass and no momentum or the other way around." + +-- Scalar-vector multiplication. +(.*) :: Applicative f => Double -> f Double -> f Double +(.*) x ys = (* x) <$> ys + +-- NOTE: Fixed parameters without mass have a tuning parameter of NaN. +massesToTuningParameters :: Foldable f => Masses f -> AuxiliaryTuningParameters +massesToTuningParameters = VB.fromList . map (fromMaybe nan) . toList + where + nan = 0 / 0 + +-- We need the structure in order to fill it with the given parameters. +tuningParametersToMasses :: + Traversable f => + AuxiliaryTuningParameters -> + Masses f -> + Either String (Masses f) +tuningParametersToMasses xs ms = + if null xs' + then sequenceA msE + else Left "tuningParametersToMasses: Too many values." + where + (xs', msE) = mapAccumL setValue (VB.toList xs) ms + setValue [] _ = ([], Left "tuningParametersToMasses: Too few values.") + -- NOTE: Recover fixed parameters and unset their mass. + setValue (y : ys) _ = let y' = if isNaN y then Nothing else Just y in (ys, Right y') + +hTuningParametersToSettings :: + Traversable f => + TuningParameter -> + AuxiliaryTuningParameters -> + HSettings f -> + Either String (HSettings f) +hTuningParametersToSettings t ts (HSettings g v m l e tn) = + if tn == HTuneMassesAndLeapfrog + then case tuningParametersToMasses ts m of + Left err -> Left err + Right m' -> Right $ HSettings g v m' lTuned eTuned tn + else Right $ HSettings g v m lTuned eTuned tn + where + -- The larger epsilon, the larger the proposal step size and the lower the + -- expected acceptance ratio. + -- + -- Further, we roughly keep \( L * \epsilon = 1.0 \). The equation is not + -- correct, because we pull L closer to the original value to keep the + -- runtime somewhat acceptable. + lTuned = ceiling $ fromIntegral l / (t ** 0.9) :: Int + eTuned = t * e + +hamiltonianSimpleWithTuningParameters :: + (Applicative f, Traversable f) => + HSettings f -> + TuningParameter -> + AuxiliaryTuningParameters -> + Either String (ProposalSimple (Positions f)) +hamiltonianSimpleWithTuningParameters s t ts = case hTuningParametersToSettings t ts s of + Left err -> Left err + Right s' -> Right $ hamiltonianSimple s' + +hamiltonianSimple :: + (Applicative f, Traversable f) => + HSettings f -> + ProposalSimple (Positions f) +hamiltonianSimple (HSettings gradient mVal masses l e _) theta g = do + phi <- generateMomenta masses g + lRan <- uniformR (lL, lR) g + eRan <- uniformR (eL, eR) g + case leapfrog gradient mVal masses lRan eRan theta phi of + Nothing -> return (theta, 0.0, 1.0) + Just (theta', phi') -> + let prPhi = priorMomenta masses phi + -- NOTE: Neal page 12: In order for the proposal to be in detailed + -- balance, the momenta have to be negated before proposing the new value. + -- This is not required here since the prior involves normal distributions + -- centered around 0. However, if the multivariate normal distribution is + -- used, it makes a difference. + prPhi' = priorMomenta masses phi' + kernelR = prPhi' / prPhi + in return (theta', kernelR, 1.0) + where + lL = maximum [1 :: Int, floor $ (0.8 :: Double) * fromIntegral l] + lR = maximum [lL, ceiling $ (1.2 :: Double) * fromIntegral l] + eL = 0.8 * e + eR = 1.2 * e + +minVariance :: Double +minVariance = 1e-6 + +maxVariance :: Double +maxVariance = 1e6 + +minSamples :: Int +minSamples = 60 + +computeAuxiliaryTuningParameters :: + Foldable f => + VB.Vector (Positions f) -> + AuxiliaryTuningParameters -> + AuxiliaryTuningParameters +computeAuxiliaryTuningParameters xss ts = + VB.zipWith (\t -> rescueWith t . calcSamplesAndVariance) ts xssT + where + -- TODO: Improve matrix transposition. + xssT = VB.fromList $ M.toColumns $ M.fromLists $ VB.toList $ VB.map toList xss + calcSamplesAndVariance xs = (VB.length $ VB.uniq $ S.gsort xs, S.variance xs) + rescueWith t (sampleSize, var) = + if var < minVariance || maxVariance < var || sampleSize < minSamples + then -- then traceShow ("Rescue with " <> show t) t + t + else + let t' = sqrt (t * recip var) + in -- in traceShow ("Old mass " <> show t <> " new mass " <> show t') t' + t' + +-- | Hamiltonian Monte Carlo proposal. +-- +-- The 'Applicative' and 'Traversable' instances are used for element-wise +-- operations. +-- +-- Assume a zip-like 'Applicative' instance so that cardinality remains +-- constant. +-- +-- NOTE: The desired acceptance rate is 0.65, although the dimension of the +-- proposal is high. +-- +-- NOTE: The speed of this proposal can change drastically when tuned because +-- the leapfrog trajectory length is changed. +hamiltonian :: + (Applicative f, Traversable f) => + -- | The sample state is used to calculate the dimension of the proposal. + f Double -> + HSettings f -> + PName -> + PWeight -> + Proposal (f Double) +hamiltonian x s n w = case checkHSettings s of + Just err -> error err + Nothing -> + let desc = PDescription "Hamiltonian Monte Carlo (HMC)" + dim = PSpecial (length x) 0.65 + ts = massesToTuningParameters (hMasses s) + ps = hamiltonianSimple s + p' = Proposal n desc dim w ps + fT = defaultTuningFunction dim + tS = hTune s + fTs = + if tS == HTuneMassesAndLeapfrog + then computeAuxiliaryTuningParameters + else \_ xs -> xs + in case tS of + HNoTune -> p' Nothing + _ -> p' $ Just $ Tuner 1.0 fT ts fTs (hamiltonianSimpleWithTuningParameters s) 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 @@ -26,7 +26,7 @@ -- The actual proposal with tuning parameter. The tuning parameter does not -- change the mean. -scaleSimple :: Shape -> Scale -> TuningParameter -> ProposalSimple Double +scaleSimple :: Shape Double -> Scale Double -> TuningParameter -> ProposalSimple Double scaleSimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) @@ -36,10 +36,18 @@ where jac _ = Exp . log . recip --- | Multiplicative proposal with gamma distributed kernel. +-- | Multiplicative proposal. +-- +-- The gamma distribution is used to sample the multiplier. Therefore, this and +-- all derived proposals are log-additive in that they do not change the sign of +-- the state. Further, the value zero is never proposed when having a strictly +-- positive value. +-- +-- Consider using 'Mcmc.Proposal.Slide.slide' to allow proposition of values +-- having opposite sign. scale :: - Shape -> - Scale -> + Shape Double -> + Scale Double -> PName -> PWeight -> Tune -> @@ -48,12 +56,12 @@ where description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th --- | Multiplicative proposal with gamma distributed kernel. +-- | See 'scale'. -- -- The scale of the gamma distribution is set to (shape)^{-1}, so that the mean -- of the gamma distribution is 1.0. scaleUnbiased :: - Shape -> + Shape Double -> PName -> PWeight -> Tune -> @@ -62,7 +70,11 @@ where description = PDescription $ "Scale unbiased; shape: " ++ show k -scaleContrarilySimple :: Shape -> Scale -> TuningParameter -> ProposalSimple (Double, Double) +scaleContrarilySimple :: + Shape Double -> + Scale Double -> + TuningParameter -> + ProposalSimple (Double, Double) scaleContrarilySimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) @@ -73,13 +85,13 @@ contra (x, y) u = (x * u, y / u) jac _ u = Exp $ log $ recip $ u * u --- | Multiplicative proposal with gamma distributed kernel. +-- | See 'scale'. -- -- The two values are scaled contrarily so that their product stays constant. -- Contrary proposals are useful when parameters are confounded. scaleContrarily :: - Shape -> - Scale -> + Shape Double -> + Scale Double -> PName -> PWeight -> Tune -> 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 @@ -26,14 +26,16 @@ import Statistics.Distribution.Uniform -- The actual proposal with tuning parameter. -slideSimple :: Mean -> StandardDeviation -> TuningParameter -> ProposalSimple Double +slideSimple :: Mean Double -> StandardDeviation Double -> TuningParameter -> ProposalSimple Double slideSimple m s t = genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing --- | Additive proposal with normally distributed kernel. +-- | Additive proposal. +-- +-- A normal distribution is used to sample the addend. slide :: - Mean -> - StandardDeviation -> + Mean Double -> + StandardDeviation Double -> PName -> PWeight -> Tune -> @@ -43,15 +45,17 @@ description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s -- The actual proposal with tuning parameter. -slideSymmetricSimple :: StandardDeviation -> TuningParameter -> ProposalSimple Double +slideSymmetricSimple :: StandardDeviation Double -> TuningParameter -> ProposalSimple Double slideSymmetricSimple s t = genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing --- | Additive proposal with normally distributed kernel with mean zero. This --- proposal is very fast, because the Metropolis-Hastings-Green ratio does not --- include calculation of the forwards and backwards kernels. +-- | See 'slide'. +-- +-- Use a normal distribution with mean zero. This proposal is fast, because the +-- Metropolis-Hastings-Green ratio does not include calculation of the forwards +-- and backwards kernels. slideSymmetric :: - StandardDeviation -> + StandardDeviation Double -> PName -> PWeight -> Tune -> @@ -65,9 +69,11 @@ slideUniformSimple d t = genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing --- | Additive proposal with uniformly distributed kernel with mean zero. This --- proposal is very fast, because the Metropolis-Hastings-Green ratio does not --- include calculation of the forwards and backwards kernels. +-- | See 'slide'. +-- +-- Use a uniformly distributed kernel with mean zero. This proposal is fast, +-- because the Metropolis-Hastings-Green ratio does not include calculation of +-- the forwards and backwards kernels. slideUniformSymmetric :: Size -> PName -> @@ -82,20 +88,22 @@ contra (x, y) u = (x + u, y - u) slideContrarilySimple :: - Mean -> - StandardDeviation -> + Mean Double -> + StandardDeviation Double -> TuningParameter -> ProposalSimple (Double, Double) slideContrarilySimple m s t = genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing --- | Additive proposal with normally distributed kernel. +-- | See 'slide'. -- +-- Use a normally distributed kernel. +-- -- The two values are slid contrarily so that their sum stays constant. Contrary -- proposals are useful when parameters are confounded. slideContrarily :: - Mean -> - StandardDeviation -> + Mean Double -> + StandardDeviation Double -> PName -> PWeight -> Tune -> diff --git a/src/Mcmc/Settings.hs b/src/Mcmc/Settings.hs --- a/src/Mcmc/Settings.hs +++ b/src/Mcmc/Settings.hs @@ -17,7 +17,7 @@ ( -- * Data types AnalysisName (..), HasAnalysisName (..), - BurnInSpecification (..), + BurnInSettings (..), burnInIterations, Iterations (..), TraceLength (..), @@ -34,16 +34,21 @@ settingsSave, settingsLoad, settingsCheck, + settingsPrettyPrint, ) where import Data.Aeson import Data.Aeson.TH +import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import Mcmc.Logger import System.Directory import System.IO +bsInt :: Int -> BL.ByteString +bsInt = BB.toLazyByteString . BB.intDec + -- | Analysis name of the MCMC sampler. newtype AnalysisName = AnalysisName {fromAnalysisName :: String} deriving (Eq, Read, Show) @@ -56,7 +61,7 @@ getAnalysisName :: s -> AnalysisName -- | Burn in specification. -data BurnInSpecification +data BurnInSettings = -- | No burn in. NoBurnIn | -- | Burn in for a given number of iterations. @@ -75,17 +80,27 @@ BurnInWithCustomAutoTuning [Int] deriving (Eq, Read, Show) -$(deriveJSON defaultOptions ''BurnInSpecification) +$(deriveJSON defaultOptions ''BurnInSettings) --- Check if the burn in specification is valid. -burnInValid :: BurnInSpecification -> Bool +burnInPrettyPrint :: BurnInSettings -> BL.ByteString +burnInPrettyPrint NoBurnIn = + "None." +burnInPrettyPrint (BurnInWithoutAutoTuning x) = + bsInt x <> " iterations; no auto tune." +burnInPrettyPrint (BurnInWithAutoTuning x y) = + bsInt x <> " iterations; auto tune with a period of " <> bsInt y <> "." +burnInPrettyPrint (BurnInWithCustomAutoTuning xs) = + bsInt (sum xs) <> " iterations; custom auto tune periods." + +-- Check if the burn in settings are valid. +burnInValid :: BurnInSettings -> 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 :: BurnInSettings -> Int burnInIterations NoBurnIn = 0 burnInIterations (BurnInWithoutAutoTuning n) = n burnInIterations (BurnInWithAutoTuning n _) = n @@ -103,13 +118,28 @@ -- -- 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. + = -- | Automatically determine the minimum length of the trace. The value is + -- the maximum of used + -- + -- - 'Mcmc.Monitor.MonitorBatch' sizes + -- + -- - auto tune intervals during burn in TraceAuto | -- | Store a given minimum number of iterations of the chain. Store more -- iterations if required (see 'TraceAuto'). TraceMinimum Int + deriving (Eq, Show) +$(deriveJSON defaultOptions ''TraceLength) + +traceLengthPrettyPrint :: TraceLength -> BL.ByteString +traceLengthPrettyPrint TraceAuto = "Determined automatically." +traceLengthPrettyPrint (TraceMinimum x) = "Minimum length of " <> bsInt x <> "." + +validTraceLength :: TraceLength -> Bool +validTraceLength (TraceMinimum n) = n > 0 +validTraceLength _ = True + -- | Execution mode. data ExecutionMode = -- | Perform new run. @@ -132,6 +162,11 @@ class HasExecutionMode s where getExecutionMode :: s -> ExecutionMode +executionModePrettyPrint :: ExecutionMode -> BL.ByteString +executionModePrettyPrint Fail = "Fail if output files exist." +executionModePrettyPrint Overwrite = "Overwrite existing output files." +executionModePrettyPrint Continue = "Expect output files exist." + -- | Open a file honoring the execution mode. -- -- Call 'error' if execution mode is @@ -194,11 +229,16 @@ $(deriveJSON defaultOptions ''SaveMode) +saveModePrettyPrint :: SaveMode -> BL.ByteString +saveModePrettyPrint NoSave = "Do not save analysis." +saveModePrettyPrint Save = "Save analysis." + -- | Settings of an MCMC sampler. data Settings = Settings { sAnalysisName :: AnalysisName, - sBurnIn :: BurnInSpecification, + sBurnIn :: BurnInSettings, sIterations :: Iterations, + sTraceLength :: TraceLength, sExecutionMode :: ExecutionMode, sParallelizationMode :: ParallelizationMode, sSaveMode :: SaveMode, @@ -270,13 +310,14 @@ -- | Current iteration. Int -> IO () -settingsCheck s@(Settings nm bi i em _ _ _ _) iCurrent +settingsCheck s@(Settings nm bi i tl em _ _ _ _) iCurrent | null (fromAnalysisName nm) = serr "Analysis name is the empty string." | burnInIterations bi < 0 = serr "Number of burn in iterations is negative." - | not $ burnInValid bi = serr $ "Burn in setting invalid: " <> show bi + | 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." + | not $ validTraceLength tl = serr $ "Trace length invalid: " <> show tl <> "." | iCurrent /= 0 && em /= Continue = serr "Current iteration is non-zero but execution mode is not 'Continue'." | iCurrent == 0 && em == Continue = @@ -284,3 +325,24 @@ | otherwise = return () where serr = settingsError s iCurrent + +logModePrettyPrint :: LogMode -> BL.ByteString +logModePrettyPrint LogStdOutAndFile = "Log to standard output and file." +logModePrettyPrint LogStdOutOnly = "Log to standard output only." +logModePrettyPrint LogFileOnly = "Log to file only." + +-- | Pretty print settings. +settingsPrettyPrint :: Settings -> BL.ByteString +settingsPrettyPrint (Settings nm bi is tl em pm sm lm vb) = + BL.unlines + [ "The MCMC settings are:", + " Analysis name: " <> BL.pack (fromAnalysisName nm) <> ".", + " Burn in: " <> burnInPrettyPrint bi, + " Iterations: " <> bsInt (fromIterations is) <> " iterations.", + " Trace length: " <> traceLengthPrettyPrint tl, + " Execution mode: " <> executionModePrettyPrint em, + " Parallelization mode: " <> BL.pack (show pm) <> ".", + " Save mode: " <> saveModePrettyPrint sm, + " Log mode: " <> logModePrettyPrint lm, + " Verbosity: " <> BL.pack (show vb) <> "." + ] diff --git a/src/Mcmc/Statistics/Types.hs b/src/Mcmc/Statistics/Types.hs --- a/src/Mcmc/Statistics/Types.hs +++ b/src/Mcmc/Statistics/Types.hs @@ -24,22 +24,22 @@ where -- | Mean of a distribution. -type Mean = Double +type Mean a = a -- | Standard deviation of a distribution. -type StandardDeviation = Double +type StandardDeviation a = a -- | Variance of a distribution. -type Variance = Double +type Variance a = a -- | Shape of a distribution. -type Shape = Double +type Shape a = a -- | Scale of a distribution. -type Scale = Double +type Scale a = a -- | Rate of a distribution. -type Rate = Double +type Rate a = a -- | Dimension of a distribution. type Dimension = Int @@ -50,7 +50,7 @@ type Size = Double -- | Lower boundary of a distribution. -type LowerBoundary = Double +type LowerBoundary a = a -- | Upper boundary of a distribution. -type UpperBoundary = Double +type UpperBoundary a = a diff --git a/test/Mcmc/ProposalSpec.hs b/test/Mcmc/ProposalSpec.hs --- a/test/Mcmc/ProposalSpec.hs +++ b/test/Mcmc/ProposalSpec.hs @@ -14,6 +14,7 @@ ) where +import Mcmc.Cycle import Mcmc.Proposal import Mcmc.Proposal.Slide import System.Random.MWC diff --git a/test/Mcmc/SaveSpec.hs b/test/Mcmc/SaveSpec.hs --- a/test/Mcmc/SaveSpec.hs +++ b/test/Mcmc/SaveSpec.hs @@ -47,23 +47,27 @@ mon :: Monitor Double mon = Monitor monStd [] [] +settings :: Settings +settings = + Settings + (AnalysisName "SaveSpec") + (BurnInWithAutoTuning 20 10) + (Iterations 100) + TraceAuto + Overwrite + Sequential + NoSave + LogStdOutOnly + Quiet + spec :: Spec spec = do describe "save and load" $ it "doesn't change the MCMC chain" $ do gen <- R.create - let s = - Settings - (AnalysisName "SaveSpec") - (BurnInWithAutoTuning 20 10) - (Iterations 200) - Overwrite - Sequential - NoSave - LogStdOutOnly - Quiet - c <- fromMHG <$> mhg noPrior lh proposals mon TraceAuto 0 gen + a <- mhg settings noPrior lh proposals mon 0 gen + c <- fromMHG <$> mcmc settings a savedChain <- toSavedChain c c' <- fromSavedChain noPrior lh proposals mon savedChain putStrLn "@load . save@ should be @id@." @@ -76,8 +80,8 @@ -- g1' <- R.save $ generator c' -- g1 `shouldBe` g1' putStrLn "Sampling from the chains should be the same." - r <- fromMHG <$> mcmc s (MHG c) - r' <- fromMHG <$> mcmc s (MHG c') + r <- fromMHG <$> mcmcContinue (Iterations 100) settings (MHG c) + r' <- fromMHG <$> mcmcContinue (Iterations 100) settings (MHG c') link r `shouldBe` link r' iteration r `shouldBe` iteration r' frozenT2 <- freezeT (trace c) @@ -87,20 +91,22 @@ g2' <- R.save $ generator r' g2 `shouldBe` g2' --- -- TODO: 'mhContinue'. --- describe "mhContinue" --- $ it "mh 200 + mhContinue 200 == mh 400" --- $ do --- gen1 <- create --- let s1 = chain "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 400 gen1 --- r1 <- mh s1 --- gen2 <- create --- let s2 = chain "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 200 gen2 --- r2' <- mh s2 --- r2 <- mhContinue 200 r2' --- link r1 `shouldBe` link r2 --- iteration r1 `shouldBe` iteration r2 --- trace r1 `shouldBe` trace r2 --- g <- save $ generator r1 --- g' <- save $ generator r2 --- g `shouldBe` g' + describe "mhContinue" $ + it "mcmc 50 + mcmcContinue 50 == mcmc 100" $ + do + gen1 <- R.create + a1 <- mhg settings noPrior lh proposals mon 0 gen1 + r1 <- fromMHG <$> mcmc settings a1 + gen2 <- R.create + let settings' = settings {sIterations = Iterations 50} + a2 <- mhg settings' noPrior lh proposals mon 0 gen2 + r2' <- mcmc settings' a2 + r2 <- fromMHG <$> mcmcContinue (Iterations 50) settings' r2' + link r1 `shouldBe` link r2 + iteration r1 `shouldBe` iteration r2 + frozenT1 <- freezeT (trace r1) + frozenT2 <- freezeT (trace r2) + frozenT1 `shouldBe` frozenT2 + g <- R.save $ generator r1 + g' <- R.save $ generator r2 + g `shouldBe` g'