diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -4,6 +4,13 @@ ## Unreleased changes +- Provide unsafe loader for MHG algorithm (useful for initializing chains with + different prior or likelihood functions from saves). +- Documentation and readme. +- Log normal prior distribution. +- Parallel computation of prior and likelihood. Speculative parallelization; + this change is not always beneficial, we will see. + ## 0.6.2.4 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -7,22 +7,22 @@ At the moment, the following algorithms are available: -- Metropolis-Hastings-Green 1; +- Metropolis-Hastings-Green 1; - Metropolis-coupled Markov chain Monte Carlo (also known as parallel - tempering) 2 , 3. -- Hamilton Monte Carlo proposal 4. + tempering) 2 , 3. +- Hamilton Monte Carlo proposal 4. ## Documentation -The [source code](https://hackage.haskell.org/package/mcmc) contains detailed documentation about general concepts as well +The [source code](https://hackage.haskell.org/package/mcmc/docs/Mcmc.html) contains detailed documentation about general concepts as well as specific functions. ## Examples -[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. +The Git repository also includes [example MCMC analyses](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples). Build them with +[cabal-install](https://cabal.readthedocs.io/en/latest/cabal-commands.html#) or [Stack](https://docs.haskellstack.org/en/stable/README/). git clone https://github.com/dschrempf/mcmc.git cd mcmc @@ -32,7 +32,7 @@ stack exec archery -For a more involved example, have a look at the [phylogenetic dating project](https://github.com/dschrempf/mcmc-dating). +For a more involved example, have a look at a [phylogenetic dating project](https://github.com/dschrempf/mcmc-dating). # Footnotes diff --git a/mcmc.cabal b/mcmc.cabal --- a/mcmc.cabal +++ b/mcmc.cabal @@ -1,6 +1,6 @@ cabal-version: 3.0 name: mcmc -version: 0.6.2.4 +version: 0.6.2.5 synopsis: Sample from a posterior using Markov chain Monte Carlo description: Please see the README on GitHub at @@ -101,6 +101,7 @@ , math-functions , microlens , mwc-random + , parallel , pretty-show , primitive , statistics diff --git a/src/Mcmc.hs b/src/Mcmc.hs --- a/src/Mcmc.hs +++ b/src/Mcmc.hs @@ -18,22 +18,28 @@ -- Markov Chain Monte Carlo (pp. 45), CRC press. -- -- This library focusses on classical Markov chain Monte Carlo algorithms such --- as the Metropolis-Hastings-Green [1] algorithm, or population methods +-- as the Metropolis-Hastings-Green (MHG) [1] algorithm, or population methods -- involving parallel chains such as the Metropolic-coupled Markov chain Monte -- Carlo [2] algorithm. In particular, sequential Monte Carlo [3] algorithms --- following a moving posterior distribution are not provided. +-- following a moving posterior distribution are not provided. Recently, +-- Hamiltonian Monte Carlo (HMC) proposals have been added [4]. HMC proposals +-- can be used with [automatic +-- differentiation](https://hackage.haskell.org/package/ad). HMC proposals with +-- automatic differentiation are quite slow for complicated prior or likelihood +-- functions, but they are incredibly useful when specialized MHG proposals are +-- not readily available. -- -- An MCMC sampler can be run with 'mcmc', for example using the -- Metropolis-Hastings-Green algorithm 'mhg'. -- -- Usually, it is best to start with an example: -- --- - Basic inference of the [accuracy of an --- archer](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples/Archery/Archery.hs) +-- - Basic inference of the [accuracy of an archer](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples/Archery/Archery.hs) (see the [statement of the problem](https://revbayes.github.io/tutorials/mcmc/archery.html)). -- --- - [More involved --- examples](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples/PhylogeneticLikelihood/PhylogeneticLikelihood.hs) +-- - [Other examples](https://github.com/dschrempf/mcmc/tree/master/mcmc-examples). -- +-- - More involved example performing [phylogenetic dating](https://github.com/dschrempf/mcmc-date). +-- -- __The import of this module alone should cover most use cases.__ -- -- @[1]@ Geyer, C. J. (2011), Introduction to markov chain monte carlo, In @@ -45,6 +51,9 @@ -- -- @[3]@ Sequential monte carlo methods in practice (2001), Editors: Arnaud -- Doucet, Nando de Freitas, and Neil Gordon, Springer New York. +-- +-- @[4]@ Review by Betancourt and notes: Betancourt, M., A conceptual +-- introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017). module Mcmc ( -- * Proposals @@ -172,12 +181,12 @@ module Mcmc.Monitor.Parameter, module Mcmc.Monitor.ParameterBatch, - -- * Prior, likelihood, and posterior values and functions + -- * Priors, likelihoodl, and posteriors module Mcmc.Prior, module Mcmc.Likelihood, module Mcmc.Posterior, - -- * MCMC samplers + -- * Run MCMC samplers mcmc, mcmcContinue, -- | See also 'settingsLoad', 'mhgLoad', and 'mc3Load'. @@ -189,7 +198,7 @@ -- * Marginal likelihood calculation module Mcmc.MarginalLikelihood, - -- * Useful types + -- * Types used in statistics module Mcmc.Statistics.Types, -- * Useful re-exports 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 @@ -22,6 +22,7 @@ mhg, mhgSave, mhgLoad, + mhgLoadUnsafe, MHGRatio, mhgAccept, ) @@ -30,6 +31,7 @@ import Codec.Compression.GZip import Control.Monad import Control.Monad.IO.Class +import Control.Parallel.Strategies import Data.Aeson import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time @@ -124,9 +126,35 @@ Monitor a -> AnalysisName -> IO (MHG a) -mhgLoad pr lh cc mn nm = do +mhgLoad = mhgLoadWith fromSavedChain + +-- | See 'mhgLoad' but do not perform sanity checks. +-- +-- Useful when restarting a run with changed prior function, likelihood function +-- or proposals. Use with care! +mhgLoadUnsafe :: + FromJSON a => + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + AnalysisName -> + IO (MHG a) +mhgLoadUnsafe = mhgLoadWith fromSavedChainUnsafe + +-- Nice type :-). +mhgLoadWith :: + FromJSON a => + (PriorFunction a -> LikelihoodFunction a -> Cycle a -> Monitor a -> SavedChain a -> IO (Chain a)) -> + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + AnalysisName -> + IO (MHG a) +mhgLoadWith f pr lh cc mn nm = do savedChain <- eitherDecode . decompress <$> BL.readFile (mhgFn nm) - chain <- either error (fromSavedChain pr lh cc mn) savedChain + chain <- either error (f pr lh cc mn) savedChain return $ MHG chain -- | MHG ratios are stored in log domain. @@ -178,23 +206,15 @@ -- 1. Sample new state. (!y, !q, !j) <- liftIO $ s x g -- 2. Calculate Metropolis-Hastings-Green ratio. - let !pY = pF y - !lY = lF y - !r = mhgRatio (pX * lX) (pY * lY) q j + -- + -- Most often, parallelization is not helpful, because the prior and + -- likelihood functions are too fast; see + -- https://stackoverflow.com/a/46603680/3536806. + let (pY, lY) = (pF y, lF y) `using` parTuple2 rdeepseq rdeepseq + -- let !pY = pF y + -- !lY = lF y + let !r = mhgRatio (pX * lX) (pY * lY) q j -- 3. Accept or reject. - -- if ln r >= 0.0 - -- then do - -- let !ac' = pushA p True ac - -- return $ MHG $ c {link = Link y pY lY, acceptance = ac'} - -- else do - -- b <- uniform g - -- if b < exp (ln r) - -- then do - -- let !ac' = pushA p True ac - -- return $ MHG $ c {link = Link y pY lY, acceptance = ac'} - -- else do - -- let !ac' = pushA p False ac - -- return $ MHG $ c {acceptance = pushA p False ac'} accept <- mhgAccept r g if accept then do 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 @@ -20,6 +20,7 @@ ( SavedChain (..), toSavedChain, fromSavedChain, + fromSavedChainUnsafe, ) where @@ -78,6 +79,10 @@ -- | Load a saved chain. -- +-- Perform some safety checks: +-- +-- Check that the number of proposals is equal. +-- -- Recompute and check the prior and likelihood for the last state because the -- functions may have changed. Of course, we cannot test for the same function, -- but having the same prior and likelihood at the last state is already a good @@ -89,7 +94,7 @@ Monitor a -> SavedChain a -> IO (Chain a) -fromSavedChain pr lh cc mn (SavedChain ci it i tr ac' g' ts) +fromSavedChain pr lh cc mn sv | pr (state it) /= prior it = let msg = unlines @@ -99,11 +104,40 @@ ] in error msg | lh (state it) /= likelihood it = - error "fromSave: Provided likelihood function does not match the saved likelihood." - | otherwise = do - g <- loadGen g' - tr' <- thawT tr - return $ Chain ci it i tr' ac g i pr lh cc' mn + let msg = + unlines + [ "fromSave: Provided likelihood function does not match the saved likelihood function.", + "fromSave: Current likelihood:" <> show (likelihood it) <> ".", + "fromSave: Given likelihood:" <> show (lh $ state it) <> "." + ] + in error msg + | length (fromAcceptance ac) /= length (ccProposals cc) = + let msg = + unlines + [ "fromSave: The number of proposals does not match.", + "fromSave: Number of saved proposals:" <> show (length $ fromAcceptance ac) <> ".", + "fromSave: Number of given proposals:" <> show (length $ ccProposals cc) <> "." + ] + in error msg + | otherwise = fromSavedChainUnsafe pr lh cc mn sv + where + it = savedLink sv + ac = savedAcceptance sv + +-- | See 'fromSavedChain' but do not perform sanity checks. Useful when +-- restarting a run with changed prior function, likelihood function or +-- proposals. +fromSavedChainUnsafe :: + PriorFunction a -> + LikelihoodFunction a -> + Cycle a -> + Monitor a -> + SavedChain a -> + IO (Chain a) +fromSavedChainUnsafe pr lh cc mn (SavedChain ci it i tr ac' g' ts) = do + g <- loadGen g' + tr' <- thawT tr + return $ Chain ci it i tr' ac g i pr lh cc' mn where ac = transformKeysA [0 ..] (ccProposals cc) ac' tunePs mt p = case mt of diff --git a/src/Mcmc/Internal/Gamma.hs b/src/Mcmc/Internal/Gamma.hs --- a/src/Mcmc/Internal/Gamma.hs +++ b/src/Mcmc/Internal/Gamma.hs @@ -28,6 +28,8 @@ mEulerMascheroniG :: RealFloat a => a mEulerMascheroniG = 0.5772156649015328606065121 +-- | Generalized version of the log gamma distribution. See +-- 'Numeric.SpecFunctions.logGamma'. logGammaG :: (Typeable a, RealFloat a) => a -> a logGammaG z | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z diff --git a/src/Mcmc/Prior.hs b/src/Mcmc/Prior.hs --- a/src/Mcmc/Prior.hs +++ b/src/Mcmc/Prior.hs @@ -31,6 +31,7 @@ gammaMeanOne, gammaShapeScaleToMeanVariance, gammaMeanVarianceToShapeScale, + logNormal, normal, uniform, @@ -155,6 +156,25 @@ mLnSqrt2Pi :: RealFloat a => a mLnSqrt2Pi = 0.9189385332046727417803297364056176398613974736377834128171 {-# INLINE mLnSqrt2Pi #-} + +-- | Log normal distributed prior. +-- +-- NOTE: The log normal distribution is parametrized with the mean \(\mu\) and +-- the standard deviation \(\sigma\) of the underlying normal distribution. The +-- mean and variance of the log normal distribution itself are functions of +-- \(\mu\) and \(\sigma\), but are not the same as \(\mu\) and \(\sigma\)! +-- +-- Call 'error' if the standard deviation is zero or negative. +logNormal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a +logNormal m s x + | s <= 0 = error "logNormal: Standard deviation is zero or negative." + | x <= 0 = 0 + | otherwise = Exp $ t + e + where + t = negate $ mLnSqrt2Pi + log (x * s) + a = recip $ 2 * s * s + b = log x - m + e = negate $ a * b * b -- | Normal distributed prior. -- diff --git a/src/Mcmc/Settings.hs b/src/Mcmc/Settings.hs --- a/src/Mcmc/Settings.hs +++ b/src/Mcmc/Settings.hs @@ -99,7 +99,7 @@ burnInPrettyPrint (BurnInWithAutoTuning x y) = bsInt x <> " iterations; auto tune with a period of " <> bsInt y <> "." burnInPrettyPrint (BurnInWithCustomAutoTuning xs ys) = - bsInt (sum xs) <> " fast," <> bsInt (sum ys) <> " slow iterations; custom auto tune periods." + bsInt (sum xs) <> " fast, " <> bsInt (sum ys) <> " slow iterations; custom auto tune periods." -- Check if the burn in settings are valid. burnInValid :: BurnInSettings -> Bool