mcmc 0.2.4 → 0.3.0

raw patch · 21 files changed

+710/−297 lines, 21 filesdep +dirichletdep +math-functionsdep +primitivePVP ok

version bump matches the API change (PVP)

Dependencies added: dirichlet, math-functions, primitive

API changes (from Hackage documentation)

Files

• ChangeLog.md +26/−0
• bench/Normal.hs +2/−5
• bench/Poisson.hs +4/−4
• mcmc.cabal +9/−5
• src/Mcmc.hs +9/−1
• src/Mcmc/Internal/Shuffle.hs +64/−0
• src/Mcmc/Mcmc.hs +65/−28
• src/Mcmc/Metropolis.hs +14/−12
• src/Mcmc/Monitor/Parameter.hs +1/−1
• src/Mcmc/Prior.hs +17/−7
• src/Mcmc/Proposal.hs +68/−32
• src/Mcmc/Proposal/Bactrian.hs +26/−20
• src/Mcmc/Proposal/Generic.hs +46/−41
• src/Mcmc/Proposal/Scale.hs +43/−20
• src/Mcmc/Proposal/Simplex.hs +203/−0
• src/Mcmc/Proposal/Slide.hs +41/−30
• src/Mcmc/Save.hs +19/−11
• src/Mcmc/Status.hs +46/−9
• src/Mcmc/Tools/Shuffle.hs +0/−64
• test/Mcmc/ProposalSpec.hs +2/−2
• test/Mcmc/SaveSpec.hs +5/−5

@@ -5,6 +5,32 @@ ## Unreleased changes  +## 0.3.0++-   New shorter example/test for dating trees.+-   `noData` allows running a chain without likelihood function.+-   Give proposal parameters `PName`, `PDescription`, and `PWeight` newtype+    wrappers.+-   Give `Tune` a data type.+-   Allow periodical cleansing of state (`Cleaner`).+-   Add description string to proposals, so that they can be identified in an+    easier way.+-   Add simplices and proposals on simplices.+-   `slideUniform` renamed to `slideUniformSymmetric`.+-   Merge tools into internal.+-   Do not export internal modules.+++## 0.2.4++-   **Change order of arguments for proposals**.+-   'slideStem' was renamed to 'slideBranch'.+-   Change ProposalSimple from newtype to type.+-   Contravariant instances of parameter and batch monitors. Use `(>$<)` instead+    of `(@.)` and `(@#)`.+-   Add `gammaDirichlet` prior for partitioned dating analyses.++ ## 0.2.3  -   Contrary proposals.

@@ -37,7 +37,7 @@ proposals :: Cycle Double proposals =   fromList-    [slideSymmetric 1.0 "medium" 1 True]+    [slideSymmetric 1.0 (PName "Medium") (PWeight 1) Tune]  mons :: [MonitorParameter Double] mons = [monitorDouble "mu"]@@ -45,9 +45,6 @@ monStd :: MonitorStdOut Double monStd = monitorStdOut mons 200 --- monFile :: MonitorFile Double--- monFile = monitorFile "Mu" mons 200- mon :: Monitor Double mon = Monitor monStd [] [] @@ -68,7 +65,7 @@ proposalsBactrian :: Cycle Double proposalsBactrian =   fromList-    [slideBactrian 0.5 1.0 "bactrian" 1 True]+    [slideBactrian 0.5 1.0 (PName "Bactrian") (PWeight 1) Tune]  normalBactrianBench :: GenIO -> IO () normalBactrianBench g = do

@@ -39,19 +39,19 @@     m = sum ys / fromIntegral (length ys)  f :: Int -> Double -> I -> Log Double-f ft yr (alpha, beta) = Exp $ logProbability (poisson l) (fromIntegral ft)+f ft yr (a, b) = Exp $ logProbability (poisson l) (fromIntegral ft)   where-    l = exp $ alpha + beta * yr+    l = exp $ a + b * yr  lh :: I -> Log Double lh x =   product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]  proposalAlpha :: Proposal I-proposalAlpha = _1 @~ slideSymmetric 0.2 "alpha" 2 False+proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (PWeight 1) NoTune  proposalBeta :: Proposal I-proposalBeta = _2 @~ slideSymmetric 0.2 "beta" 1 False+proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (PWeight 1) NoTune  proposals :: Cycle I proposals = fromList [proposalAlpha, proposalBeta]

@@ -1,6 +1,6 @@ cabal-version:  2.2 name:           mcmc-version:        0.2.4+version:        0.3.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@@ -23,7 +23,6 @@ library   exposed-modules:       Mcmc-      Mcmc.Internal.ByteString       Mcmc.Item       Mcmc.Mcmc       Mcmc.Metropolis@@ -38,12 +37,14 @@       Mcmc.Proposal.Generic       Mcmc.Proposal.Scale       Mcmc.Proposal.Slide+      Mcmc.Proposal.Simplex       Mcmc.Save       Mcmc.Status-      Mcmc.Tools.Shuffle       Mcmc.Trace       Mcmc.Verbosity   other-modules:+      Mcmc.Internal.ByteString+      Mcmc.Internal.Shuffle       Paths_mcmc   autogen-modules:       Paths_mcmc@@ -57,10 +58,13 @@     , containers     , data-default     , directory+    , dirichlet     , double-conversion     , log-domain+    , math-functions     , microlens     , mwc-random+    , primitive     , statistics     , time     , transformers@@ -77,7 +81,7 @@       Paths_mcmc   hs-source-dirs:       test-  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N+  ghc-options: -Wall   build-depends:       QuickCheck     , base >=4.7 && <5@@ -100,7 +104,7 @@       Paths_mcmc   hs-source-dirs:       bench-  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N+  ghc-options: -Wall   build-depends:       base >=4.7 && <5     , criterion

@@ -82,17 +82,21 @@     -- B. R., Huelsenbeck, J. P., …, Revbayes: bayesian phylogenetic inference using     -- graphical models and an interactive model-specification language, Systematic     -- Biology, 65(4), 726–736 (2016). http://dx.doi.org/10.1093/sysbio/syw021+    PName (..),+    PWeight (..),     Proposal,     (@~),+    Tune (..),     scale,     scaleUnbiased,     scaleContrarily,     scaleBactrian,     slide,     slideSymmetric,-    slideUniform,+    slideUniformSymmetric,     slideContrarily,     slideBactrian,+    module Mcmc.Proposal.Simplex,     Cycle,     fromList,     Order (..),@@ -113,10 +117,13 @@     -- space (see above) and to monitor the MCMC run, as well as some auxiliary     -- information.     status,+    Cleaner (..),+    cleanWith,     saveWith,     force,     quiet,     debug,+    noData,      -- * Monitor @@ -160,6 +167,7 @@ import Mcmc.Proposal import Mcmc.Proposal.Bactrian import Mcmc.Proposal.Scale+import Mcmc.Proposal.Simplex import Mcmc.Proposal.Slide import Mcmc.Save import Mcmc.Status

@@ -0,0 +1,64 @@+-- |+-- Module      :  Mcmc.Internal.Shuffle+-- Description :  Shuffle a list+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 20 14:37:09 2020.+--+-- From https://wiki.haskell.org/Random_shuffle.+module Mcmc.Internal.Shuffle+  ( shuffle,+    shuffleN,+    grabble,+  )+where++import Control.Monad+import Control.Monad.ST+import Data.Vector (Vector)+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as M+import System.Random.MWC+  ( GenIO,+    uniformR,+  )++-- | Shuffle a list.+shuffle :: [a] -> GenIO -> IO [a]+shuffle xs g = head <$> grabble xs 1 (length xs) g++-- | Shuffle a list @n@ times.+shuffleN :: [a] -> Int -> GenIO -> IO [[a]]+shuffleN xs n = grabble xs n (length xs)++-- -- Using System.Random.Shuffle. Speed is the same, so stay without additional dependency.+-- -- | Shuffle a list @n@ times.+-- shuffleN :: [a] -> Int -> GenIO -> IO [[a]]+-- shuffleN xs n g = replicateM n $ fmap (shuffle xs) (rseqM (length xs - 1) g)+--   where+--     rseqM :: Int -> GenIO -> IO [Int]+--     rseqM 0 _ = return []+--     rseqM i gen = liftM2 (:) (uniformR (0, i) gen) (rseqM (i - 1) gen)++-- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@+-- elements from @xs@, without replacement, and that @m@ times.+grabble :: [a] -> Int -> Int -> GenIO -> IO [[a]]+grabble xs m n gen = do+  swapss <- replicateM m $+    forM [0 .. min (l - 1) n] $ \i -> do+      j <- uniformR (i, l) gen+      return (i, j)+  return $ map (V.toList . V.take n . swapElems (V.fromList xs)) swapss+  where+    l = length xs - 1++swapElems :: Vector a -> [(Int, Int)] -> Vector a+swapElems xs swaps = runST $ do+  mxs <- V.unsafeThaw xs+  mapM_ (uncurry $ M.unsafeSwap mxs) swaps+  V.unsafeFreeze mxs

@@ -15,15 +15,16 @@ -- Functions to work with the 'Mcmc' state transformer. module Mcmc.Mcmc   ( Mcmc,-    mcmcOutT,+    mcmcOutB,     mcmcOutS,-    mcmcWarnT,+    mcmcWarnB,     mcmcWarnS,-    mcmcInfoT,+    mcmcInfoB,     mcmcInfoS,-    mcmcDebugT,+    mcmcDebugB,     mcmcDebugS,     mcmcAutotune,+    mcmcClean,     mcmcResetA,     mcmcSummarizeCycle,     mcmcReport,@@ -40,12 +41,14 @@ import Data.Maybe import Data.Time.Clock import Data.Time.Format+import Mcmc.Item import Mcmc.Monitor import Mcmc.Monitor.Time import Mcmc.Proposal import Mcmc.Save import Mcmc.Status hiding (debug) import Mcmc.Verbosity+import Numeric.Log import System.Directory import System.IO import Prelude hiding (cycle)@@ -58,66 +61,96 @@ msgPrepare c t = BL.cons c $ ": " <> t  -- | Write to standard output and log file.-mcmcOutT :: BL.ByteString -> Mcmc a ()-mcmcOutT msg = do+mcmcOutB :: BL.ByteString -> Mcmc a ()+mcmcOutB msg = do   h <- fromMaybe (error "mcmcOut: Log handle is missing.") <$> gets logHandle   liftIO $ BL.putStrLn msg >> BL.hPutStrLn h msg  -- | Write to standard output and log file. mcmcOutS :: String -> Mcmc a ()-mcmcOutS = mcmcOutT . BL.pack+mcmcOutS = mcmcOutB . BL.pack  -- Perform warning action. mcmcWarnA :: Mcmc a () -> Mcmc a () mcmcWarnA a = gets verbosity >>= \v -> info v a  -- | Print warning message.-mcmcWarnT :: BL.ByteString -> Mcmc a ()-mcmcWarnT = mcmcWarnA . mcmcOutT . msgPrepare 'W'+mcmcWarnB :: BL.ByteString -> Mcmc a ()+mcmcWarnB = mcmcWarnA . mcmcOutB . msgPrepare 'W'  -- | Print warning message. mcmcWarnS :: String -> Mcmc a ()-mcmcWarnS = mcmcWarnT . BL.pack+mcmcWarnS = mcmcWarnB . BL.pack  -- Perform info action. mcmcInfoA :: Mcmc a () -> Mcmc a () mcmcInfoA a = gets verbosity >>= \v -> info v a  -- | Print info message.-mcmcInfoT :: BL.ByteString -> Mcmc a ()-mcmcInfoT = mcmcInfoA . mcmcOutT . msgPrepare 'I'+mcmcInfoB :: BL.ByteString -> Mcmc a ()+mcmcInfoB = mcmcInfoA . mcmcOutB . msgPrepare 'I'  -- | Print info message. mcmcInfoS :: String -> Mcmc a ()-mcmcInfoS = mcmcInfoT . BL.pack+mcmcInfoS = mcmcInfoB . BL.pack  -- Perform debug action. mcmcDebugA :: Mcmc a () -> Mcmc a () mcmcDebugA a = gets verbosity >>= \v -> debug v a  -- | Print debug message.-mcmcDebugT :: BL.ByteString -> Mcmc a ()-mcmcDebugT = mcmcDebugA . mcmcOutT . msgPrepare 'D'+mcmcDebugB :: BL.ByteString -> Mcmc a ()+mcmcDebugB = mcmcDebugA . mcmcOutB . msgPrepare 'D'  -- | Print debug message. mcmcDebugS :: String -> Mcmc a ()-mcmcDebugS = mcmcDebugT . BL.pack+mcmcDebugS = mcmcDebugB . BL.pack  -- | Auto tune the 'Proposal's in the 'Cycle' of the chain. Reset acceptance counts. -- See 'autotuneCycle'. mcmcAutotune :: Mcmc a () mcmcAutotune = do-  mcmcDebugT "Auto tune."+  mcmcDebugB "Auto tune."   s <- get   let a = acceptance s       c = cycle s       c' = autotuneCycle a c   put $ s {cycle = c'} +-- | Clean the state.+mcmcClean :: Mcmc a ()+mcmcClean = do+  s <- get+  let cl = cleaner s+      i = iteration s+  case cl of+    Just (Cleaner n f) | i `mod` n == 0 -> do+      mcmcDebugB "Clean state."+      let (Item st pr lh) = item s+      mcmcDebugS $+        "Old log prior and log likelihood: " ++ show (ln pr) ++ ", " ++ show (ln lh) ++ "."+      let prF = priorF s+          lhF = likelihoodF s+          st' = f st+          pr' = prF st'+          lh' = lhF st'+      mcmcDebugS $+        "New log prior and log likelihood: " ++ show (ln pr') ++ ", " ++ show (ln lh') ++ "."+      let dLogPr = abs $ ln pr - ln pr'+          dLogLh = abs $ ln lh - ln lh'+      when+        (dLogPr > 0.01)+        (mcmcWarnS $ "Log of old and new prior differ by " ++ show dLogPr ++ ".")+      when+        (dLogPr > 0.01)+        (mcmcWarnS $ "Log of old and new likelihood differ by " ++ show dLogLh ++ ".")+      put $ s {item = Item st' pr' lh'}+    _ -> return ()+ -- | Reset acceptance counts. mcmcResetA :: Mcmc a () mcmcResetA = do-  mcmcDebugT "Reset acceptance ratios."+  mcmcDebugB "Reset acceptance ratios."   s <- get   let a = acceptance s   put $ s {acceptance = resetA a}@@ -150,7 +183,7 @@         (True, _, _) -> openFile lfn AppendMode   put s {logHandle = mh}   mcmcDebugS $ "Log file name: " ++ lfn ++ "."-  mcmcDebugT "Log file opened."+  mcmcDebugB "Log file opened."  -- Set the total number of iterations, the current time and open the 'Monitor's -- of the chain. See 'mOpen'.@@ -176,14 +209,18 @@   let b = burnInIterations s       t = autoTuningPeriod s       n = iterations s+      c = cleaner s   case b of     Just b' -> mcmcInfoS $ "Burn in for " <> show b' <> " iterations."     Nothing -> return ()   case t of     Just t' -> mcmcInfoS $ "Auto tune every " <> show t' <> " iterations (during burn in only)."     Nothing -> return ()+  case c of+    Just (Cleaner c' _) -> mcmcInfoS $ "Clean state every " <> show c' <> " iterations."+    Nothing -> return ()   mcmcInfoS $ "Run chain for " <> show n <> " iterations."-  mcmcInfoT "Initial state."+  mcmcInfoB "Initial state."   mcmcMonitorExec  -- Save the status of an MCMC run. See 'saveStatus'.@@ -192,11 +229,11 @@   s <- get   case save s of     Just n -> do-      mcmcInfoT $ "Save Markov chain with trace of length " <> BL.pack (show n) <> "."-      mcmcInfoT "For long traces, or complex objects, this may take a while."+      mcmcInfoB $ "Save Markov chain with trace of length " <> BL.pack (show n) <> "."+      mcmcInfoB "For long traces, or complex objects, this may take a while."       liftIO $ saveStatus (name s <> ".mcmc") s-      mcmcInfoT "Done saving Markov chain."-    Nothing -> mcmcInfoT "Do not save the Markov chain."+      mcmcInfoB "Done saving Markov chain."+    Nothing -> mcmcInfoB "Do not save the Markov chain."  -- | Execute the 'Monitor's of the chain. See 'mExec'. mcmcMonitorExec :: ToJSON a => Mcmc a ()@@ -209,14 +246,14 @@       tr = trace s       vb = verbosity s   mt <- liftIO $ mExec vb i ss st tr j m-  forM_ mt mcmcOutT+  forM_ mt mcmcOutB  -- Close the 'Monitor's of the chain. See 'mClose'. mcmcClose :: ToJSON a => Mcmc a () mcmcClose = do   s <- get-  mcmcSummarizeCycle >>= mcmcInfoT-  mcmcInfoT "Metropolis-Hastings sampler finished."+  mcmcSummarizeCycle >>= mcmcInfoB+  mcmcInfoB "Metropolis-Hastings sampler finished."   let m = monitor s   m' <- liftIO $ mClose m   put $ s {monitor = m'}@@ -225,7 +262,7 @@   let rt = case start s of         Nothing -> error "mcmcClose: Start time not set."         Just (_, st) -> t `diffUTCTime` st-  mcmcInfoT $ "Wall clock run time: " <> renderDuration rt <> "."+  mcmcInfoB $ "Wall clock run time: " <> renderDuration rt <> "."   mcmcInfoS $ "End time: " <> fTime t   case logHandle s of     Just h -> liftIO $ hClose h

@@ -46,9 +46,10 @@ -- Handbook of Markov Chain Monte Carlo), or (b) almost surely reject the -- proposal when either fY or q are zero (Chapter 1). Since I trust the author -- of Chapter 1 (Charles Geyer) I choose to follow option (b).-mhRatio :: Log Double -> Log Double -> Log Double -> Log Double--- q = qYX / qXY-mhRatio fX fY q = fY * q / fX+mhRatio :: Log Double -> Log Double -> Log Double -> Log Double -> Log Double+-- q = qYX / qXY * jXY; see 'ProposalSimple'.+-- j = Jacobian.+mhRatio fX fY q j = fY / fX * q * j {-# INLINE mhRatio #-}  mhPropose :: Proposal a -> Mcmc a ()@@ -61,11 +62,11 @@       a = acceptance s       g = generator s   -- 1. Sample new state.-  (!y, !q) <- liftIO $ p x g+  (!y, !q, !j) <- liftIO $ p x g   -- 2. Calculate Metropolis-Hastings ratio.   let !pY = pF y       !lY = lF y-      !r = mhRatio (pX * lX) (pY * lY) q+      !r = mhRatio (pX * lX) (pY * lY) q j   -- 3. Accept or reject.   if ln r >= 0.0     then put $ s {item = Item y pY lY, acceptance = pushA m True a}@@ -86,6 +87,7 @@       t = trace s       n = iteration s   put $ s {trace = pushT i t, iteration = succ n}+  mcmcClean   mcmcMonitorExec  -- Run N iterations.@@ -104,13 +106,13 @@   | b > t = do     mcmcResetA     mhNIter t-    mcmcSummarizeCycle >>= mcmcDebugT+    mcmcSummarizeCycle >>= mcmcDebugB     mcmcAutotune     mhBurnInN (b - t) (Just t)   | otherwise = do     mcmcResetA     mhNIter b-    mcmcSummarizeCycle >>= mcmcInfoT+    mcmcSummarizeCycle >>= mcmcInfoB     mcmcInfoS $ "Acceptance ratios calculated over the last " <> show b <> " iterations." mhBurnInN b Nothing = mhNIter b @@ -123,7 +125,7 @@     mcmcInfoS $ "Burn in for " <> show b <> " cycles."     mcmcDebugS $ "Auto tuning period is " <> show t <> "."     mhBurnInN b t-    mcmcInfoT "Burn in finished."+    mcmcInfoB "Burn in finished."  -- Run for given number of iterations. mhRun :: ToJSON a => Int -> Mcmc a ()@@ -142,8 +144,8 @@  mhT :: ToJSON a => Mcmc a () mhT = do-  mcmcInfoT "Metropolis-Hastings sampler."-  mcmcSummarizeCycle >>= mcmcInfoT+  mcmcInfoB "Metropolis-Hastings sampler."+  mcmcSummarizeCycle >>= mcmcInfoB   mcmcReport   s <- get   let b = fromMaybe 0 (burnInIterations s)@@ -152,9 +154,9 @@  mhContinueT :: ToJSON a => Int -> Mcmc a () mhContinueT dn = do-  mcmcInfoT "Continuation of Metropolis-Hastings sampler."+  mcmcInfoB "Continuation of Metropolis-Hastings sampler."   mcmcInfoS $ "Run chain for " <> show dn <> " additional iterations."-  mcmcSummarizeCycle >>= mcmcInfoT+  mcmcSummarizeCycle >>= mcmcInfoB   mhRun dn  -- | Continue a Markov chain for a given number of Metropolis-Hastings steps.

@@ -43,7 +43,7 @@     mpFunc :: a -> BB.Builder   } -instance Contravariant (MonitorParameter) where+instance Contravariant MonitorParameter where   contramap f (MonitorParameter n m) = MonitorParameter n (m . f)  -- | Convert a parameter monitor from one data type to another.

@@ -13,7 +13,9 @@ -- Creation date: Thu Jul 23 13:26:14 2020. module Mcmc.Prior   ( -- * Continuous priors+    largerThan,     positive,+    lowerThan,     negative,     uniform,     normal,@@ -35,17 +37,25 @@ import qualified Statistics.Distribution.Gamma as S import qualified Statistics.Distribution.Normal as S --- | Improper uniform prior; larger than 0.+-- | Improper uniform prior; strictly larger than a given value.+largerThan :: Double -> Double -> Log Double+largerThan a x+  | x <= a = 0+  | otherwise = 1++-- | Improper uniform prior; strictly larger than zero. positive :: Double -> Log Double-positive x-  | x <= 0 = 0+positive = largerThan 0++-- | Improper uniform prior; strictly lower than a given value.+lowerThan :: Double -> Double -> Log Double+lowerThan b x+  | x >= b = 0   | otherwise = 1 --- | Improper uniform prior; lower than 0.+-- | Improper uniform prior; strictly lower than zero. negative :: Double -> Log Double-negative x-  | x >= 0 = 0-  | otherwise = 1+negative = lowerThan 0  -- | Uniform prior on [a, b]. uniform ::

@@ -2,8 +2,6 @@ {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE RankNTypes #-} --- TODO: Proposals on simplices: SimplexElementScale (?).- -- | -- Module      :  Mcmc.Proposal -- Description :  Proposals and cycles@@ -17,10 +15,14 @@ -- Creation date: Wed May 20 13:42:53 2020. module Mcmc.Proposal   ( -- * Proposal+    PName (..),+    PDescription (..),+    PWeight (..),     Proposal (..),     (@~),     ProposalSimple,     Tuner (tParam, tFunc),+    Tune (..),     createProposal,     tune, @@ -57,36 +59,54 @@ import Data.Maybe import Lens.Micro import Mcmc.Internal.ByteString-import Mcmc.Tools.Shuffle+import Mcmc.Internal.Shuffle import Numeric.Log hiding (sum) import System.Random.MWC +-- | Proposal name.+newtype PName = PName {fromPName :: String}+  deriving (Show, Eq, Ord)++-- | Proposal description.+newtype PDescription = PDescription {fromPDescription :: String}+  deriving (Show, Eq, Ord)++-- | The weight determines how often a 'Proposal' is executed per iteration of+-- the Markov chain.+newtype PWeight = PWeight {fromPWeight :: Int}+  deriving (Show, Eq, Ord)+ -- | 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 -- conditioned on the current state (i.e., a kernel). -- -- A 'Proposal' may be tuneable in that it contains information about how to enlarge -- or shrink the step size to tune the acceptance ratio.+--+-- No proposals with the same name and description are allowed in a 'Cycle'. data Proposal a = Proposal-  { -- | Name (no proposals with the same name are allowed in a 'Cycle').-    pName :: String,+  { -- | Name of the affected variable.+    pName :: PName,+    -- | Description of the proposal type and parameters.+    pDescription :: PDescription,     -- | The weight determines how often a 'Proposal' is executed per iteration of     -- the Markov chain.-    pWeight :: Int,+    pWeight :: PWeight,     -- | Simple proposal without name, weight, and tuning information.     pSimple :: ProposalSimple a,     -- | Tuning is disabled if set to 'Nothing'.     pTuner :: Maybe (Tuner a)   } +-- XXX: This should be removed. instance Show (Proposal a) where-  show m = show $ pName m+  show m = fromPName (pName m) <> " " <> fromPDescription (pDescription m) <> ", weight " <> show (fromPWeight $ pWeight m)  instance Eq (Proposal a) where-  m == n = pName m == pName n+  m == n = pName m == pName n && pDescription m == pDescription n  instance Ord (Proposal a) where-  compare = compare `on` pName+  compare = compare `on` (\p -> (pDescription p, pName p, pWeight p))  -- | Convert a proposal from one data type to another using a lens. --@@ -96,7 +116,7 @@ -- scaleFirstEntryOfTuple = _1 @~ scale -- @ (@~) :: Lens' b a -> Proposal a -> Proposal b-(@~) l (Proposal n w s t) = Proposal n w (convertS l s) (convertT l <$> t)+(@~) l (Proposal n d w s t) = Proposal n d w (convertS l s) (convertT l <$> t)  -- | Simple proposal without tuning information. --@@ -105,16 +125,26 @@ -- -- In order to calculate the Metropolis-Hastings ratio, we need to know the -- ratio of the backward to forward kernels (i.e., the probability masses or--- probability densities). For unbiased proposals, this ratio is 1.0. For biased--- proposals, the ratio is ?.-type ProposalSimple a = a -> GenIO -> IO (a, Log Double)+-- probability densities) and the absolute value of the determinant of the+-- Jacobian matrix.+--+-- For unbiased proposals, these values are 1.0 such that+--+-- @+-- proposalSimpleUnbiased x g = return (x', 1.0, 1.0)+-- @+--+-- For biased proposals, the kernel ratio is qYX / qXY, where qXY is the+-- probability density to move from X to Y, and the absolute value of the+-- determinant of the Jacobian matrix differs from 1.0.+type ProposalSimple a = a -> GenIO -> IO (a, Log Double, Log Double)  convertS :: Lens' b a -> ProposalSimple a -> ProposalSimple b convertS l s = s'   where     s' v g = do-      (x', r) <- s (v ^. l) g-      return (set l x' v, r)+      (x', r, j) <- s (v ^. l) g+      return (set l x' v, r, j)  -- | Tune the acceptance ratio of a 'Proposal'; see 'tune', or 'autotuneCycle'. data Tuner a = Tuner@@ -127,21 +157,27 @@   where     f' x = convertS l $ f x +-- | Tune the proposal?+data Tune = Tune | NoTune+  deriving (Show, Eq)+ -- | Create a possibly tuneable proposal. createProposal ::+  -- | Description of the proposal type and parameters.+  PDescription ->   -- | Function creating a simple proposal for a given tuning parameter. The   -- larger the tuning parameter, the larger the proposal (and the lower the   -- expected acceptance ratio), and vice versa.   (Double -> ProposalSimple a) ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Activate tuning?-  Bool ->+  Tune ->   Proposal a-createProposal f n w True = Proposal n w (f 1.0) (Just $ Tuner 1.0 f)-createProposal f n w False = Proposal n w (f 1.0) Nothing+createProposal d f n w Tune = Proposal n d w (f 1.0) (Just $ Tuner 1.0 f)+createProposal d f n w NoTune = Proposal n d w (f 1.0) Nothing  -- Minimal tuning parameter; subject to change. tuningParamMin :: Double@@ -223,11 +259,9 @@ fromList [] =   error "fromList: Received an empty list but cannot create an empty Cycle." fromList xs =-  if length (nub nms) == length nms+  if length (nub xs) == length xs     then Cycle xs def-    else error "fromList: Proposals don't have unique names."-  where-    nms = map pName xs+    else error "fromList: Proposals are not unique."  -- | Set the order of 'Proposal's in a 'Cycle'. setOrder :: Order -> Cycle a -> Cycle a@@ -243,7 +277,7 @@     return [psR ++ reverse psR | psR <- psRs]   SequentialReversibleO -> return $ replicate n $ ps ++ reverse ps   where-    !ps = concat [replicate (pWeight m) m | m <- xs]+    !ps = concat [replicate (fromPWeight $ pWeight m) m | m <- xs]  -- | Tune 'Proposal's in the 'Cycle'. See 'tune'. tuneCycle :: Map (Proposal a) Double -> Cycle a -> Cycle a@@ -271,10 +305,12 @@   BL.ByteString ->   BL.ByteString ->   BL.ByteString ->+  BL.ByteString ->   BL.ByteString-renderRow name weight nAccept nReject acceptRatio tuneParam manualAdjustment = "   " <> nm <> wt <> na <> nr <> ra <> tp <> mt+renderRow name ptype weight nAccept nReject acceptRatio tuneParam manualAdjustment = "   " <> nm <> pt <> wt <> na <> nr <> ra <> tp <> mt   where     nm = alignLeft 30 name+    pt = alignLeft 50 ptype     wt = alignRight 8 weight     na = alignRight 15 nAccept     nr = alignRight 15 nReject@@ -284,12 +320,13 @@  proposalHeader :: BL.ByteString proposalHeader =-  renderRow "Proposal" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter" "Consider manual adjustment"+  renderRow "Name" "Description" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter" "Consider manual adjustment"  summarizeProposal :: Proposal a -> Maybe (Int, Int, Double) -> BL.ByteString summarizeProposal m r =   renderRow-    (BL.pack name)+    (BL.pack $ fromPName $ pName m)+    (BL.pack $ fromPDescription $ pDescription m)     weight     nAccept     nReject@@ -297,8 +334,7 @@     tuneParamStr     manualAdjustmentStr   where-    name = pName m-    weight = BB.toLazyByteString $ BB.intDec $ pWeight m+    weight = BB.toLazyByteString $ BB.intDec $ fromPWeight $ pWeight m     nAccept = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _1)) r     nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) r     acceptRatio = BL.fromStrict $ maybe "" (BC.toFixed 3 . (^. _3)) r@@ -324,7 +360,7 @@       ++ [hLine proposalHeader]   where     ps = ccProposals c-    mpi = BB.toLazyByteString $ BB.intDec $ sum $ map pWeight ps+    mpi = BB.toLazyByteString $ BB.intDec $ sum $ map (fromPWeight . pWeight) ps     ar m = acceptanceRatio m a  -- | For each key @k@, store the number of accepted and rejected proposals.

@@ -52,10 +52,10 @@   Double ->   Double ->   GenIO ->-  IO (Double, Log Double)+  IO (Double, Log Double, Log Double) bactrianAdditive m s x g = do   dx <- genBactrian m s g-  return (x + dx, 1.0)+  return (x + dx, 1.0, 1.0)  -- bactrianSimple lens spike stdDev tune forwardOp backwardOp bactrianAdditiveSimple ::@@ -75,9 +75,9 @@ -- The [Bactrian -- kernel](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/figure/fig01) is -- a mixture of two symmetrically arranged normal distributions. The spike--- parameter loosely determines the standard deviations of the individual humps--- while the second parameter refers to the standard deviation of the complete--- Bactrian kernel.+-- parameter (0, 1) loosely determines the standard deviations (>0.0) of the+-- individual humps while the second parameter refers to the standard deviation+-- of the complete Bactrian kernel. -- -- See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3845170/. slideBactrian ::@@ -86,13 +86,15 @@   -- | Standard deviation.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-slideBactrian m s = createProposal (bactrianAdditiveSimple m s)+slideBactrian m s = createProposal description (bactrianAdditiveSimple m s)+  where+    description = PDescription $ "Slide Bactrian; spike: " ++ show m ++ ", sd: " ++ show s  -- We have: -- x  (1+dx ) = x'@@ -108,12 +110,14 @@   Double ->   Double ->   GenIO ->-  IO (Double, Log Double)+  IO (Double, Log Double, Log Double) bactrianMult m s x g = do-  dx <- genBactrian m s g-  let qXY = logDensityBactrian m s dx-      qYX = logDensityBactrian m s (fInv dx)-  return (x * (1 + dx), qYX / qXY)+  du <- genBactrian m s g+  let qXY = logDensityBactrian m s du+      qYX = logDensityBactrian m s (fInv du)+      u = 1.0 + du+      jac = Exp $ log $ recip u+  return (x * u, qYX / qXY, jac)  bactrianMultSimple :: Double -> Double -> Double -> ProposalSimple Double bactrianMultSimple m s t@@ -130,10 +134,12 @@   -- | Standard deviation.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-scaleBactrian m s = createProposal (bactrianMultSimple m s)+scaleBactrian m s = createProposal description (bactrianMultSimple m s)+  where+    description = PDescription $ "Scale Bactrian; spike: " ++ show m <> ", sd: " <> show s

@@ -18,58 +18,54 @@ import Mcmc.Proposal import Numeric.Log import Statistics.Distribution-import System.Random.MWC -sampleCont ::-  (ContDistr d, ContGen d) =>-  d ->-  (a -> Double -> a) ->-  Maybe (Double -> Double) ->-  a ->-  GenIO ->-  IO (a, Log Double)-sampleCont d f mfInv x g = do-  dx <- genContVar d g-  let r = case mfInv of-        Nothing -> 1.0-        Just fInv ->-          let qXY = Exp $ logDensity d dx-              qYX = Exp $ logDensity d (fInv dx)-           in qYX / qXY-  return (x `f` dx, r)-{-# INLINEABLE sampleCont #-}- -- | Generic function to create proposals for continuous parameters ('Double'). genericContinuous ::   (ContDistr d, ContGen d) =>   -- | Probability distribution   d ->-  -- | Forward operator, e.g. (+), so that x + dx = x'.+  -- | Forward operator.+  --+  -- For example, for a multiplicative proposal on one variable the forward+  -- operator is @(*)@, so that @x * u = y@.   (a -> Double -> a) ->-  -- | Inverse operator, e.g., 'negate', so that x' + (negate dx) = x. Only-  -- required for biased proposals.+  -- | Inverse operator.+  --+  -- For example, 'recip' for a multiplicative proposal on one variable, since+  -- @y * (recip u) = x * u * (recip u) = x@.+  --+  -- Required for biased proposals.   Maybe (Double -> Double) ->+  -- | Function to compute the absolute value of the determinant of the Jacobian+  -- matrix. For example, for a multiplicative proposal on one variable, we have+  --+  -- @+  -- detJacobian _ u = Exp $ log $ recip u+  -- @+  --+  -- 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.+  --+  -- Conversion to log domain is necessary, because some determinants of+  -- Jacobians are very small (or large).+  Maybe (a -> Double -> Log Double) ->   ProposalSimple a-genericContinuous d f fInv = sampleCont d f fInv--sampleDiscrete ::-  (DiscreteDistr d, DiscreteGen d) =>-  d ->-  (a -> Int -> a) ->-  Maybe (Int -> Int) ->-  a ->-  GenIO ->-  IO (a, Log Double)-sampleDiscrete d f mfInv x g = do-  dx <- genDiscreteVar d g-  let r = case mfInv of+genericContinuous d f mInv mJac x g = do+  u <- genContVar d g+  let r = case mInv of         Nothing -> 1.0         Just fInv ->-          let qXY = Exp $ logProbability d dx-              qYX = Exp $ logProbability d (fInv dx)+          let qXY = Exp $ logDensity d u+              qYX = Exp $ logDensity d (fInv u)            in qYX / qXY-  return (x `f` dx, r)-{-# INLINEABLE sampleDiscrete #-}+      j = case mJac of+        Nothing -> 1.0+        Just fJac -> fJac x u+  return (x `f` u, r, j)+{-# INLINEABLE genericContinuous #-}  -- | Generic function to create proposals for discrete parameters ('Int'). genericDiscrete ::@@ -82,4 +78,13 @@   -- required for biased proposals.   Maybe (Int -> Int) ->   ProposalSimple a-genericDiscrete fd f fInv = sampleDiscrete fd f fInv+genericDiscrete d f mfInv x g = do+  u <- genDiscreteVar d g+  let r = case mfInv of+        Nothing -> 1.0+        Just fInv ->+          let qXY = Exp $ logProbability d u+              qYX = Exp $ logProbability d (fInv u)+           in qYX / qXY+  return (x `f` u, r, 1.0)+{-# INLINEABLE genericDiscrete #-}

@@ -20,12 +20,20 @@  import Mcmc.Proposal import Mcmc.Proposal.Generic+import Numeric.Log import Statistics.Distribution.Gamma  -- The actual proposal with tuning parameter. The tuning parameter does not -- change the mean. scaleSimple :: Double -> Double -> Double -> ProposalSimple Double-scaleSimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) (*) (Just recip)+scaleSimple k th t =+  genericContinuous+    (gammaDistr (k / t) (th * t))+    (*)+    (Just recip)+    (Just jac)+  where+    jac _ = Exp . log . recip  -- | Multiplicative proposal with Gamma distributed kernel. scale ::@@ -34,13 +42,15 @@   -- | Scale.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-scale k th = createProposal (scaleSimple k th)+scale k th = createProposal description (scaleSimple k th)+  where+    description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th  -- | Multiplicative proposal with Gamma distributed kernel. --@@ -50,20 +60,31 @@   -- | Shape.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-scaleUnbiased k = createProposal (scaleSimple k (1 / k))--contra :: (Double, Double) -> Double -> (Double, Double)-contra (x, y) z = (x * z, y / z)+scaleUnbiased k = createProposal description (scaleSimple k (1 / k))+  where+    description = PDescription $ "Scale unbiased; shape: " ++ show k  scaleContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)-scaleContrarilySimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) contra (Just recip)+scaleContrarilySimple k th t =+  genericContinuous+    (gammaDistr (k / t) (th * t))+    contra+    (Just recip)+    (Just jac)+  where+    contra (x, y) u = (x * u, y / u)+    jac _ u = Exp $ log $ recip $ u * u +-- -- Determinant of Jacobian matrix.+-- contraJac :: (Double, Double) -> Double+-- contraJac (x, y) = x * y+ -- | Multiplicative proposal with Gamma distributed kernel. -- -- The two values are scaled contrarily so that their product stays constant.@@ -74,10 +95,12 @@   -- | Scale.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal (Double, Double)-scaleContrarily k th = createProposal (scaleContrarilySimple k th)+scaleContrarily k th = createProposal description (scaleContrarilySimple k th)+  where+    description = PDescription $ "Scale contrariliy; shape: " ++ show k ++ ", scale: " ++ show th

@@ -0,0 +1,203 @@+{-# LANGUAGE TemplateHaskell #-}++-- |+-- Module      :  Mcmc.Proposal.Simplex+-- Description :  Proposals on simplices+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon Oct 19 15:32:31 2020.+module Mcmc.Proposal.Simplex+  ( -- * Elements of simplices+    Simplex (toVector),+    simplexUniform,+    simplexFromVector,++    -- * Proposals on simplices+    dirichlet,+    beta,+  )+where++import Data.Aeson+import Data.Aeson.TH+import qualified Data.Vector.Unboxed as V+import Mcmc.Proposal+import Numeric.Log+import Statistics.Distribution+import Statistics.Distribution.Beta+import Statistics.Distribution.Dirichlet++-- import Debug.Trace++-- | An element of a simplex.+--+-- A vector of non-negative values summing to one.+--+-- The nomenclature is not very consistent, because a K-dimensional simplex is+-- usually considered to be the set containing all @K@-dimensional vectors with+-- non-negative elements that sum to 1.0. However, I couldn't come up with a+-- better name. Maybe @SimplexElement@, but that was too long.+newtype Simplex = SimplexUnsafe {toVector :: V.Vector Double}+  deriving (Eq, Show)++$(deriveJSON defaultOptions ''Simplex)++-- Tolerance.+eps :: Double+eps = 1e-14++-- Check if vector is normalized with tolerance 'eps'.+isNormalized :: V.Vector Double -> Bool+isNormalized v+  | abs (V.sum v - 1.0) > eps = False+  | otherwise = True++-- Check if vector contains negative elements.+isNegative :: V.Vector Double -> Bool+isNegative = V.any (< 0)++-- | Ensure that the value vector is an element of a simplex.+--+-- Return 'Left' if:+-- - The value vector is empty.+-- - The value vector contains negative elements.+-- - The value vector is not normalized.+simplexFromVector :: V.Vector Double -> Either String Simplex+simplexFromVector v+  | V.null v = Left "simplexFromVector: Vector is empty."+  | isNegative v = Left "simplexFromVector: Vector contains negative elements."+  | not (isNormalized v) = Left "simplexFromVector: Vector is not normalized."+  | otherwise = Right $ SimplexUnsafe v++-- | Create the uniform element of the K-dimensional simplex.+--+-- Set all values to \(1/D\).+simplexUniform :: Int -> Simplex+simplexUniform k = either error id $ simplexFromVector $ V.replicate k (1.0 / fromIntegral k)++-- Tuning function is inverted (high alpha means small steps).+getTuningFunction :: Double -> (Double -> Double)+getTuningFunction t = (/ t'')+  where+    -- Start with small steps.+    t' = t / 10000+    -- Extremely small tuning parameters lead to numeric overflow. The square+    -- root pulls the tuning parameter closer to 1.0. However, overflow may+    -- still occur (the involved Gamma functions grow faster than the+    -- exponential). I did not observe numeric underflow in my tests.+    t'' = sqrt t'++-- The tuning parameter is proportional to the inverted mean of the shape+-- parameter values.+--+-- The values determining the proposal size have been set using an example+-- analysis. They are good values for this analysis, but may fail for other+-- analyses.+dirichletSimple :: Double -> ProposalSimple Simplex+dirichletSimple t (SimplexUnsafe xs) g = do+  -- If @t@ is high and above 1.0, the parameter vector will be low, and the+  -- variance will be high. If @t@ is low and below 1.0, the parameter vector+  -- will be high, and the Dirichlet distribution will be very concentrated with+  -- low variance.+  let ddXs = either error id $ dirichletDistribution $ V.map tf xs+  -- traceShowM $ V.map tf xs+  ys <- dirichletSample ddXs g+  -- traceShowM ys+  -- Have to check if parameters are valid (because zeroes do occur).+  let eitherDdYs = dirichletDistribution $ V.map tf ys+  let r = case eitherDdYs of+        -- Set ratio to 0; so that the proposal will not be accepted.+        Left _ -> 0+        Right ddYs -> dirichletDensity ddYs xs / dirichletDensity ddXs ys+  -- I do not think a Jacobian is necessary in this case. I do know that if a+  -- subset of states is updated a Jacobian would be necessary.+  --+  -- traceShowM mhRatio+  return (SimplexUnsafe ys, r, 1.0)+  where+    tf = getTuningFunction t++-- | Dirichlet proposal on a simplex.+--+-- For a given element of a K-dimensional simplex, propose a new element of the+-- K-dimensional simplex. The new element is sampled from the multivariate+-- Dirichlet distribution with parameter vector being the old element of the+-- simplex.+--+-- The tuning parameter is used to determine the concentration of the Dirichlet+-- distribution: the lower the tuning parameter, the higher the concentration.+--+-- This proposal may have low acceptance ratios. In this case, please see the+-- coordinate wise 'beta' proposal.+dirichlet :: PName -> PWeight -> Tune -> Proposal Simplex+dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple++-- The tuning parameter is the inverted mean of the shape values.+--+-- The values determining the proposal size have been set using an example+-- analysis. They are good values for this analysis, but may fail for other+-- analyses.+--+-- See also the 'dirichlet' proposal.+betaSimple :: Int -> Double -> ProposalSimple Simplex+betaSimple i t (SimplexUnsafe xs) g = do+  -- Shape parameters of beta distribution. Do not assume that the sum of the+  -- elements of 'xs' is 1.0, because then repeated proposals let the sum of the+  -- vector diverge.+  let aX = xI+      bX = xsSum - xI+      bdXI = betaDistr (tf aX) (tf bX)+  -- New value of element i.+  yI <- genContVar bdXI g+  -- Shape parameters of beta distribution.+  let aY = yI+      bY = 1.0 - yI+      eitherBdYI = betaDistrE (tf aY) (tf bY)+  -- See 'dirichlet', which has the same construct.+  let r = case eitherBdYI of+        Nothing -> 0+        Just bdYI -> Exp $ logDensity bdYI xI - logDensity bdXI yI+      -- The absolute value of the determinant of the Jacobian. Derivation takes+      -- a while...+      ja1 = bY / bX+      jac = Exp $ fromIntegral (V.length xs - 2) * log ja1+  -- Construct new vector.+  let -- Normalization function for other elements.+      -- nf x = x * bY / bX+      --+      -- It turns out, that this factor is also needed to compute the determinant+      -- of the Jacobian above.+      nf x = x * ja1+      ys = V.generate (V.length xs) (\j -> if i == j then yI else nf (xs V.! j))+  return (either error id $ simplexFromVector ys, r, jac)+  where+    xI = xs V.! i+    xsSum = V.sum xs+    tf = getTuningFunction t++-- | Beta proposal on a specific coordinate @i@ on a simplex.+--+-- For a given element of a K-dimensional simplex, propose a new element of the+-- K-dimensional simplex. The coordinate @i@ of the new element is sampled from+-- the beta distribution. The other coordinates are normalized such that the+-- values sum to 1.0. The parameters of the beta distribution are chosen such+-- that the expected value of the beta distribution is the value of the old+-- coordinate.+--+-- The tuning parameter is used to determine the concentration of the beta+-- distribution: the lower the tuning parameter, the higher the concentration.+--+-- See also the 'dirichlet' proposal.+--+-- No "out of bounds" checks are performed during compile time. Run time errors+-- can occur if @i@ is negative, or if @i-1@ is larger than the length of the+-- element vector of the simplex.+beta :: Int -> PName -> PWeight -> Tune -> Proposal Simplex+beta i = createProposal description (betaSimple i)+  where+    description = PDescription $ "Beta; coordinate: " ++ show i

@@ -14,7 +14,7 @@ module Mcmc.Proposal.Slide   ( slide,     slideSymmetric,-    slideUniform,+    slideUniformSymmetric,     slideContrarily,   ) where@@ -26,7 +26,8 @@  -- The actual proposal with tuning parameter. slideSimple :: Double -> Double -> Double -> ProposalSimple Double-slideSimple m s t = genericContinuous (normalDistr m (s * t)) (+) (Just negate)+slideSimple m s t =+  genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing  -- | Additive proposal with normally distributed kernel. slide ::@@ -35,17 +36,20 @@   -- | Standard deviation.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-slide m s = createProposal (slideSimple m s)+slide m s = createProposal description (slideSimple m s)+  where+    description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s  -- The actual proposal with tuning parameter. slideSymmetricSimple :: Double -> Double -> ProposalSimple Double-slideSymmetricSimple s t = genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing+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 ratio does not include@@ -54,39 +58,44 @@   -- | Standard deviation.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-slideSymmetric s = createProposal (slideSymmetricSimple s)+slideSymmetric s = createProposal description (slideSymmetricSimple s)+  where+    description = PDescription $ "Slide symmetric; sd: " ++ show s  -- The actual proposal with tuning parameter. slideUniformSimple :: Double -> Double -> ProposalSimple Double slideUniformSimple d t =-  genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing+  genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing Nothing --- | Additive proposal with uniformly distributed kernel. This proposal is very fast,--- because the Metropolis-Hastings ratio does not include calculation of the--- forwards and backwards kernels.-slideUniform ::+-- | Additive proposal with uniformly distributed kernel with mean zero. This+-- proposal is very fast, because the Metropolis-Hastings ratio does not include+-- calculation of the forwards and backwards kernels.+slideUniformSymmetric ::   -- | Delta.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal Double-slideUniform d = createProposal (slideUniformSimple d)+slideUniformSymmetric d = createProposal description (slideUniformSimple d)+  where+    description = PDescription $ "Slide uniform symmetric; delta: " ++ show d  contra :: (Double, Double) -> Double -> (Double, Double)-contra (x, y) d = (x + d, y - d)+contra (x, y) u = (x + u, y - u)  slideContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)-slideContrarilySimple m s t = genericContinuous (normalDistr m (s * t)) contra (Just negate)+slideContrarilySimple m s t =+  genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing  -- | Additive proposal with normally distributed kernel. --@@ -98,10 +107,12 @@   -- | Standard deviation.   Double ->   -- | Name.-  String ->-  -- | Weight.-  Int ->+  PName ->+  -- | PWeight.+  PWeight ->   -- | Enable tuning.-  Bool ->+  Tune ->   Proposal (Double, Double)-slideContrarily m s = createProposal (slideContrarilySimple m s)+slideContrarily m s = createProposal description (slideContrarilySimple m s)+  where+    description = PDescription $ "Slide contrarily; mean: " ++ show m ++ ", sd: " ++ show s

@@ -68,7 +68,7 @@ $(deriveJSON defaultOptions ''Save)  toSave :: Status a -> Save a-toSave (Status nm it i tr ac br at is f sv vb g _ _ _ _ c _) =+toSave (Status nm it i tr ac br at is f sv vb g _ _ _ _ _ c _) =   Save     nm     it@@ -104,9 +104,10 @@   (a -> Log Double) ->   Cycle a ->   Monitor a ->+  Maybe (Cleaner a) ->   Save a ->   Status a-fromSave p l c m (Save nm it i tr ac' br at is f sv vb g' ts) =+fromSave pr lh cc m cl (Save nm it i tr ac' br at is f sv vb g' ts) =   Status     nm     it@@ -122,16 +123,17 @@     g     Nothing     Nothing-    p-    l-    c'+    pr+    lh+    cl+    cc'     m   where-    ac = transformKeysA [0 ..] (ccProposals c) ac'+    ac = transformKeysA [0 ..] (ccProposals cc) ac'     -- TODO: Splitmix. Remove as soon as split mix is used and is available with     -- the statistics package.     g = unsafePerformIO $ restore $ toSeed g'-    c' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals c) ts) c+    cc' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals cc) ts) cc  -- | Load a 'Status' from file. --@@ -139,23 +141,29 @@ -- a chain is restored: -- - prior function -- - likelihood function+-- - cleaning function -- - cycle -- - monitor -- -- To avoid incomplete continued runs, the @.mcmc@ file is removed after load. loadStatus ::   FromJSON a =>+  -- | Prior function.   (a -> Log Double) ->+  -- | Likelihood function.   (a -> Log Double) ->   Cycle a ->   Monitor a ->+  -- | Cleaner, if needed.+  Maybe (Cleaner a) ->+  -- | Path of status to load.   FilePath ->   IO (Status a)-loadStatus p l c m fn = do+loadStatus pr lh cc mn cl fn = do   res <- eitherDecode . decompress <$> BL.readFile fn   let s = case res of         Left err -> error err-        Right sv -> fromSave p l c m sv+        Right sv -> fromSave pr lh cc mn cl sv   -- Check if prior and likelihood matches.   let Item x svp svl = item s   -- Recompute and check the prior and likelihood for the last state because the@@ -163,10 +171,10 @@   -- function, but having the same prior and likelihood at the last state is   -- already a good indicator.   when-    (p x /= svp)+    (pr x /= svp)     (error "loadStatus: Provided prior function does not match the saved prior.")   when-    (l x /= svl)+    (lh x /= svl)     (error "loadStatus: Provided likelihood function does not match the saved likelihood.")   removeFile fn   return s

@@ -21,12 +21,15 @@ -- -- Creation date: Tue May  5 18:01:15 2020. module Mcmc.Status-  ( Status (..),+  ( Cleaner (..),+    Status (..),     status,+    cleanWith,     saveWith,     force,     quiet,     debug,+    noData,   ) where @@ -42,6 +45,28 @@ import System.Random.MWC hiding (save) import Prelude hiding (cycle) +-- | Clean the state periodically.+--+-- The prior and the likelihood will be updated after the cleaning process.+--+-- For long chains, successive numerical errors can accumulate such that the+-- state diverges from honoring specific required constraints. In these cases, a+-- 'Cleaner' can be used to ensure that the required constraints of the state+-- are honored. For example, the branches of an ultrametric phylogeny may+-- diverge slightly after successful many proposals such that the phylogeny is+-- not anymore ultrametric.+--+-- Please be aware that the Markov chain will not converge to the true posterior+-- distribution if the state is changed substantially! Only apply subtle changes+-- that are absolutely necessary to preserve the required properties of the+-- state such as specific numerical constraints.+data Cleaner a = Cleaner+  { -- | Clean every given number of iterations.+    clEvery :: Int,+    -- | Cleaning function. Executed before monitoring the state.+    clFunction :: a -> a+  }+ -- | The 'Status' contains all information to run an MCMC chain. It is -- constructed using the function 'status'. --@@ -101,6 +126,8 @@     -- | The likelihood function. The un-normalized posterior is the product of     -- the prior and the likelihood.     likelihoodF :: a -> Log Double,+    -- | Clean the state periodically.+    cleaner :: Maybe (Cleaner a),     --     -- Variables related to the algorithm; not saved. @@ -118,8 +145,8 @@   (a -> Log Double) ->   -- | The likelihood function.   (a -> Log Double) ->-  -- | A list of 'Proposal's executed in forward order. The-  -- chain will be logged after each cycle.+  -- | A list of 'Proposal's executed in forward order. The chain will be logged+  -- after each cycle.   Cycle a ->   -- | A 'Monitor' observing the chain.   Monitor a ->@@ -127,14 +154,14 @@   a ->   -- | Number of burn in iterations; deactivate burn in with 'Nothing'.   Maybe Int ->-  -- | Auto tuning period (only during burn in); deactivate-  -- auto tuning with 'Nothing'.+  -- | Auto tuning period (only during burn in); deactivate auto tuning with+  -- 'Nothing'.   Maybe Int ->-  -- | Number of normal iterations excluding burn in. Note-  -- that auto tuning only happens during burn in.+  -- | Number of normal iterations excluding burn in. Note that auto tuning only+  -- happens during burn in.   Int ->-  -- | A source of randomness. For reproducible runs, make-  -- sure to use generators with the same, fixed seed.+  -- | A source of randomness. For reproducible runs, make sure to use+  -- generators with the same, fixed seed.   GenIO ->   Status a status n p l c m x mB mT nI g@@ -157,11 +184,17 @@       Nothing       p       l+      Nothing       c       m   where     i = Item x (p x) (l x) +-- | Clean the state every given number of generations using the given function.+-- See 'Cleaner'.+cleanWith :: Cleaner a -> Status a -> Status a+cleanWith c s = s {cleaner = Just c}+ -- | Save the Markov chain with trace of given length. saveWith :: Int -> Status a -> Status a saveWith n s = s {save = Just n}@@ -179,3 +212,7 @@ -- | Be verbose. debug :: Status a -> Status a debug s = s {verbosity = Debug}++-- | Set the likelihood function to 1.0. Useful for debugging and testing.+noData :: Status a -> Status a+noData s = s {likelihoodF = const 1.0}

@@ -1,64 +0,0 @@--- |--- Module      :  Mcmc.Tools.Shuffle--- Description :  Shuffle a list--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 20 14:37:09 2020.------ From https://wiki.haskell.org/Random_shuffle.-module Mcmc.Tools.Shuffle-  ( shuffle,-    shuffleN,-    grabble,-  )-where--import Control.Monad-import Control.Monad.ST-import Data.Vector (Vector)-import qualified Data.Vector as V-import qualified Data.Vector.Mutable as M-import System.Random.MWC-  ( GenIO,-    uniformR,-  )---- | Shuffle a list.-shuffle :: [a] -> GenIO -> IO [a]-shuffle xs g = head <$> grabble xs 1 (length xs) g---- | Shuffle a list @n@ times.-shuffleN :: [a] -> Int -> GenIO -> IO [[a]]-shuffleN xs n = grabble xs n (length xs)---- -- Using System.Random.Shuffle. Speed is the same, so stay without additional dependency.--- -- | Shuffle a list @n@ times.--- shuffleN :: [a] -> Int -> GenIO -> IO [[a]]--- shuffleN xs n g = replicateM n $ fmap (shuffle xs) (rseqM (length xs - 1) g)---   where---     rseqM :: Int -> GenIO -> IO [Int]---     rseqM 0 _ = return []---     rseqM i gen = liftM2 (:) (uniformR (0, i) gen) (rseqM (i - 1) gen)---- | @grabble xs m n@ is /O(m*n')/, where @n' = min n (length xs)@. Choose @n'@--- elements from @xs@, without replacement, and that @m@ times.-grabble :: [a] -> Int -> Int -> GenIO -> IO [[a]]-grabble xs m n gen = do-  swapss <- replicateM m $-    forM [0 .. min (l - 1) n] $ \i -> do-      j <- uniformR (i, l) gen-      return (i, j)-  return $ map (V.toList . V.take n . swapElems (V.fromList xs)) swapss-  where-    l = length xs - 1--swapElems :: Vector a -> [(Int, Int)] -> Vector a-swapElems xs swaps = runST $ do-  mxs <- V.unsafeThaw xs-  mapM_ (uncurry $ M.unsafeSwap mxs) swaps-  V.unsafeFreeze mxs

@@ -20,10 +20,10 @@ import Test.Hspec  p1 :: Proposal Double-p1 = slideSymmetric 1.0 "test1" 1 True+p1 = slideSymmetric 1.0 (PName "Test 1") (PWeight 1) Tune  p2 :: Proposal Double-p2 = slideSymmetric 1.0 "test2" 3 True+p2 = slideSymmetric 1.0 (PName "Test 2") (PWeight 3) Tune  c :: Cycle Double c = fromList [p1, p2]

@@ -38,10 +38,10 @@ proposals :: Cycle Double proposals =   fromList-    [ slideSymmetric 0.1 "small" 5 True,-      slideSymmetric 1.0 "medium" 2 True,-      slideSymmetric 5.0 "large" 2 True,-      slide 1.0 4.0 "skewed" 1 True+    [ slideSymmetric 0.1 (PName "Small") (PWeight 5) Tune,+      slideSymmetric 1.0 (PName "Medium") (PWeight 2) Tune,+      slideSymmetric 5.0 (PName "Large") (PWeight 2) Tune,+      slide 1.0 4.0 (PName "Skewed") (PWeight 1) Tune     ]  monStd :: MonitorStdOut Double@@ -71,7 +71,7 @@                   saveWith 100 $                     status "SaveSpec" (const 1) lh proposals mon 0 nBurn nAutoTune nIter gen         saveStatus "SaveSpec.json" s-        s' <- loadStatus (const 1) lh proposals mon "SaveSpec.json"+        s' <- loadStatus (const 1) lh proposals mon Nothing "SaveSpec.json"         r <- mh s         r' <- mh s'         -- Done during 'loadStatus'.