mcmc 0.6.2.4 → 0.6.2.5

raw patch · 9 files changed

+135/−42 lines, 9 filesdep +parallelPVP: minor bump suggested

API additions: PVP suggests at least a minor version bump

Dependencies added: parallel

API changes (from Hackage documentation)

Files

• ChangeLog.md +7/−0
• README.md +7/−7
• mcmc.cabal +2/−1
• src/Mcmc.hs +18/−9
• src/Mcmc/Algorithm/MHG.hs +38/−18
• src/Mcmc/Chain/Save.hs +40/−6
• src/Mcmc/Internal/Gamma.hs +2/−0
• src/Mcmc/Prior.hs +20/−0
• src/Mcmc/Settings.hs +1/−1

@@ -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 

@@ -7,22 +7,22 @@  At the moment, the following algorithms are available: --   Metropolis-Hastings-Green <sup><a id="fnr.1" class="footref" href="#fn.1">1</a></sup>;+-   Metropolis-Hastings-Green <sup><a id="fnr.1" class="footref" href="#fn.1" role="doc-backlink">1</a></sup>; -   Metropolis-coupled Markov chain Monte Carlo (also known as parallel-    tempering) <sup><a id="fnr.2" class="footref" href="#fn.2">2</a></sup> <sup>, </sup><sup><a id="fnr.3" class="footref" href="#fn.3">3</a></sup>.--   Hamilton Monte Carlo proposal <sup><a id="fnr.4" class="footref" href="#fn.4">4</a></sup>.+    tempering) <sup><a id="fnr.2" class="footref" href="#fn.2" role="doc-backlink">2</a></sup> <sup>, </sup><sup><a id="fnr.3" class="footref" href="#fn.3" role="doc-backlink">3</a></sup>.+-   Hamilton Monte Carlo proposal <sup><a id="fnr.4" class="footref" href="#fn.4" role="doc-backlink">4</a></sup>.   ## 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

@@ -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 <https://github.com/dschrempf/mcmc#readme>@@ -101,6 +101,7 @@     , math-functions     , microlens     , mwc-random+    , parallel     , pretty-show     , primitive     , statistics

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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. --

@@ -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