mcmc 0.5.0.0 → 0.6.0.0

raw patch · 32 files changed

+1765/−672 lines, 32 filesdep +math-functionsdep +matricesdep ~aesondep ~bytestringdep ~circularPVP ok

version bump matches the API change (PVP)

Dependencies added: math-functions, matrices

Dependency ranges changed: aeson, bytestring, circular, containers, criterion, data-default, deepseq, directory, dirichlet, double-conversion, hspec, log-domain, microlens, monad-parallel, mwc-random, pretty-show, primitive, statistics, time, transformers, vector, zlib

API changes (from Hackage documentation)

Files

• ChangeLog.md +12/−0
• README.md +9/−3
• bench/Normal.hs +8/−4
• bench/Poisson.hs +3/−2
• mcmc.cabal +126/−111
• src/Mcmc.hs +7/−1
• src/Mcmc/Acceptance.hs +93/−0
• src/Mcmc/Algorithm.hs +5/−2
• src/Mcmc/Algorithm/MC3.hs +48/−45
• src/Mcmc/Algorithm/MHG.hs +26/−9
• src/Mcmc/Chain/Chain.hs +37/−10
• src/Mcmc/Chain/Save.hs +15/−13
• src/Mcmc/Cycle.hs +187/−0
• src/Mcmc/Environment.hs +1/−1
• src/Mcmc/Internal/Gamma.hs +195/−0
• src/Mcmc/Likelihood.hs +11/−2
• src/Mcmc/Logger.hs +4/−6
• src/Mcmc/MarginalLikelihood.hs +8/−6
• src/Mcmc/Mcmc.hs +49/−18
• src/Mcmc/Monitor/Time.hs +1/−2
• src/Mcmc/Posterior.hs +13/−1
• src/Mcmc/Prior.hs +99/−53
• src/Mcmc/Proposal.hs +144/−288
• src/Mcmc/Proposal/Bactrian.hs +16/−10
• src/Mcmc/Proposal/Generic.hs +35/−9
• src/Mcmc/Proposal/Hamiltonian.hs +448/−0
• src/Mcmc/Proposal/Scale.hs +22/−10
• src/Mcmc/Proposal/Slide.hs +25/−17
• src/Mcmc/Settings.hs +73/−11
• src/Mcmc/Statistics/Types.hs +8/−8
• test/Mcmc/ProposalSpec.hs +1/−0
• test/Mcmc/SaveSpec.hs +36/−30

@@ -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: <https://github.com/dschrempf/mcmc-date>.++ ## 0.5.0.0  -   Marginal likelihood estimation using thermodynamic integration or stepping

@@ -10,6 +10,7 @@ -   Metropolis-Hastings-Green <sup><a id="fnr.1" class="footref" href="#fn.1">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>.   ## 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  <sup><a id="fn.1" href="#fnr.1">1</a></sup> Geyer, C. J., Introduction to Markov chain Monte Carlo, In Handbook of@@ -39,8 +42,11 @@  <sup><a id="fn.2" href="#fnr.2">2</a></sup> Geyer, C. J., Markov chain monte carlo maximum likelihood, Computing Science and Statistics, Proceedings of the 23rd Symposium on the Interface,-(), (1991).+(1991).  <sup><a id="fn.3" href="#fnr.3">3</a></sup> 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).++<sup><a id="fn.4" href="#fnr.4">4</a></sup> 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.

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

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

@@ -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 <https://github.com/dschrempf/mcmc#readme>+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 <https://github.com/dschrempf/mcmc#readme>+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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@@ -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) #-}

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

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

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

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

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

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

@@ -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) <> "."+    ]

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

@@ -14,6 +14,7 @@   ) where +import Mcmc.Cycle import Mcmc.Proposal import Mcmc.Proposal.Slide import System.Random.MWC

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