mcmc 0.2.2 → 0.2.3

raw patch · 21 files changed

+376/−174 lines, 21 filesdep −asyncPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: async

API changes (from Hackage documentation)

Files

• ChangeLog.md +14/−0
• bench/Poisson.hs +2/−2
• mcmc.cabal +2/−2
• src/Mcmc.hs +5/−1
• src/Mcmc/Item.hs +8/−9
• src/Mcmc/Mcmc.hs +0/−8
• src/Mcmc/Metropolis.hs +31/−14
• src/Mcmc/Monitor.hs +39/−36
• src/Mcmc/Monitor/Parameter.hs +4/−5
• src/Mcmc/Monitor/ParameterBatch.hs +11/−14
• src/Mcmc/Monitor/Time.hs +1/−1
• src/Mcmc/Prior.hs +108/−0
• src/Mcmc/Proposal.hs +45/−26
• src/Mcmc/Proposal/Bactrian.hs +1/−1
• src/Mcmc/Proposal/Generic.hs +6/−11
• src/Mcmc/Proposal/Scale.hs +30/−4
• src/Mcmc/Proposal/Slide.hs +28/−3
• src/Mcmc/Save.hs +2/−2
• src/Mcmc/Tools/Shuffle.hs +5/−4
• test/Mcmc/ProposalSpec.hs +11/−11
• test/Mcmc/SaveSpec.hs +23/−20

@@ -5,6 +5,20 @@ ## Unreleased changes  +## 0.2.3++-   Contrary proposals.+-   Change how monitors are lifted (use normal function, not a lens).+-   Priors.+-   Remove concurrent monitors (was slow).+-   Improve MCMC sampler output.+++## 0.2.2++-   Move away from hpack.++ ## 0.2.1  -   Consistently use ByteString instead of Text.

@@ -60,10 +60,10 @@ initial = (0, 0)  monAlpha :: MonitorParameter I-monAlpha = _1 @. monitorDouble "alpha"+monAlpha = fst @. monitorDouble "alpha"  monBeta :: MonitorParameter I-monBeta = _2 @. monitorDouble "beta"+monBeta = snd @. monitorDouble "beta"  monStd :: MonitorStdOut I monStd = monitorStdOut [monAlpha, monBeta] 150

@@ -1,6 +1,6 @@ cabal-version:  2.2 name:           mcmc-version:        0.2.2+version:        0.2.3 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@@ -32,6 +32,7 @@       Mcmc.Monitor.Parameter       Mcmc.Monitor.ParameterBatch       Mcmc.Monitor.Time+      Mcmc.Prior       Mcmc.Proposal       Mcmc.Proposal.Bactrian       Mcmc.Proposal.Generic@@ -51,7 +52,6 @@   ghc-options: -Wall   build-depends:       aeson-    , async     , base >=4.7 && <5     , bytestring     , containers

@@ -71,7 +71,7 @@     --     -- @     -- slideSimple :: Lens' a Double -> Double -> Double -> Double -> ProposalSimple a-    -- slideSimple l m s t = proposalGenericContinuous l (normalDistr m (s * t)) (+) (-)+    -- slideSimple l m s t = genericContinuous l (normalDistr m (s * t)) (+) (-)     -- @     --     -- This specification is more involved. Especially since we need to know the@@ -133,6 +133,9 @@     module Mcmc.Monitor.Parameter,     module Mcmc.Monitor.ParameterBatch, +    -- * Prior distributions+    module Mcmc.Prior,+     -- * Algorithms      -- | At the moment, the library is tailored to the Metropolis-Hastings@@ -150,6 +153,7 @@ import Mcmc.Monitor import Mcmc.Monitor.Parameter import Mcmc.Monitor.ParameterBatch+import Mcmc.Prior import Mcmc.Proposal import Mcmc.Proposal.Bactrian import Mcmc.Proposal.Scale

@@ -30,15 +30,14 @@ -- | An 'Item', or link of the Markov chain. For reasons of computational -- efficiency, each state is associated with the corresponding prior and -- likelihood.-data Item a-  = Item-      { -- | The current state in the state space @a@.-        state :: a,-        -- | The current prior.-        prior :: Log Double,-        -- | The current likelihood.-        likelihood :: Log Double-      }+data Item a = Item+  { -- | The current state in the state space @a@.+    state :: a,+    -- | The current prior.+    prior :: Log Double,+    -- | The current likelihood.+    likelihood :: Log Double+  }   deriving (Eq, Ord, Show, Read)  instance ToJSON a => ToJSON (Item a) where

@@ -25,7 +25,6 @@     mcmcResetA,     mcmcSummarizeCycle,     mcmcReport,-    mcmcMonitorStdOutHeader,     mcmcMonitorExec,     mcmcRun,   )@@ -183,14 +182,7 @@     Nothing -> return ()   mcmcInfoS $ "Run chain for " <> show n <> " iterations."   mcmcInfoT "Initial state."-  mcmcMonitorStdOutHeader   mcmcMonitorExec---- | Print header line of standard output monitor.-mcmcMonitorStdOutHeader :: Mcmc a ()-mcmcMonitorStdOutHeader = do-  m <- gets monitor-  mcmcInfoA $ mcmcOutT $ mHeader m  -- Save the status of an MCMC run. See 'saveStatus'. mcmcSave :: ToJSON a => Mcmc a ()

@@ -32,10 +32,20 @@ import System.Random.MWC import Prelude hiding (cycle) --- For non-symmetric proposals.+-- The Metropolis-Hastings ratio. ----- q = qYX / qXY+-- 'Infinity' if fX is zero. In this case, the proposal is always accepted.+--+-- 'NaN' if (fY or q) and fX are zero. In this case, the proposal is always+-- rejected.++-- There is a discrepancy between authors saying that one should (a) always+-- accept the new state when the current posterior is zero (Chapter 4 of the+-- 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 {-# INLINE mhRatio #-} @@ -82,7 +92,7 @@   mcmcDebugS $ "Run " <> show n <> " iterations."   c <- gets cycle   g <- gets generator-  cycles <- liftIO $ getNCycles c n g+  cycles <- liftIO $ getNIterations c n g   forM_ cycles mhIter  -- Burn in and auto tune.@@ -91,14 +101,12 @@   | t <= 0 = error "mhBurnInN: Auto tuning period smaller equal 0."   | b > t = do     mcmcResetA-    mcmcMonitorStdOutHeader     mhNIter t     mcmcSummarizeCycle >>= mcmcDebugT     mcmcAutotune     mhBurnInN (b - t) (Just t)   | otherwise = do     mcmcResetA-    mcmcMonitorStdOutHeader     mhNIter b     mcmcSummarizeCycle >>= mcmcInfoT     mcmcInfoS $ "Acceptance ratios calculated over the last " <> show b <> " iterations."@@ -120,20 +128,21 @@ mhRun n = do   mcmcResetA   mcmcInfoS $ "Run chain for " <> show n <> " iterations."-  let (m, r) = n `quotRem` 100-  -- Print header to standard output every 100 iterations.-  replicateM_ m $ do-    mcmcMonitorStdOutHeader-    mhNIter 100-  when (r > 0) $ do-    mcmcMonitorStdOutHeader-    mhNIter r+  -- let (m, r) = n `quotRem` 100+  -- -- Print header to standard output every 100 iterations.+  -- replicateM_ m $ do+  --   mcmcMonitorStdOutHeader+  --   mhNIter 100+  -- when (r > 0) $ do+  --   mcmcMonitorStdOutHeader+  --   mhNIter r+  mhNIter n  mhT :: ToJSON a => Mcmc a () mhT = do   mcmcInfoT "Metropolis-Hastings sampler."-  mcmcReport   mcmcSummarizeCycle >>= mcmcInfoT+  mcmcReport   s <- get   let b = fromMaybe 0 (burnInIterations s)   mhBurnIn b (autoTuningPeriod s)@@ -147,6 +156,14 @@   mhRun dn  -- | Continue a Markov chain for a given number of Metropolis-Hastings steps.+--+-- At the moment, when an MCMC run is continued, the old @.mcmc@ file is+-- deleted. This behavior may change in the future.+--+-- This means that an interrupted continuation also breaks previous runs. This+-- step is necessary because, otherwise, incomplete monitor files are left on+-- disk, if a continuation is canceled. Subsequent continuations would append to+-- the incomplete monitor files and produce garbage. mhContinue ::   ToJSON a =>   -- | Additional number of Metropolis-Hastings steps.

@@ -21,16 +21,14 @@     MonitorBatch,     monitorBatch, -    -- * Use monitor+    -- * Use monitors     mOpen,     mAppend,-    mHeader,     mExec,     mClose,   ) where -import Control.Concurrent.Async import Control.Monad import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL@@ -101,6 +99,33 @@     sep = "   " <> BL.replicate (BL.length row - 3) '-'     nms = [BL.pack $ mpName p | p <- msParams m] +msDataLine ::+  Int ->+  Item a ->+  Int ->+  UTCTime ->+  Int ->+  MonitorStdOut a ->+  IO BL.ByteString+msDataLine i (Item x p l) ss st j m = do+  ct <- getCurrentTime+  let dt = ct `diffUTCTime` st+      -- Careful, don't evaluate this when i == ss.+      timePerIter = dt / fromIntegral (i - ss)+      -- -- Always 0; doesn't make much sense.+      -- tpi = if (i - ss) < 10+      --       then ""+      --       else renderDurationS timePerIter+      eta =+        if (i - ss) < 10+          then ""+          else renderDuration $ timePerIter * fromIntegral (j - i)+  return $+    msRenderRow $+      [BL.pack (show i), renderLog p, renderLog l, renderLog (p * l)]+        ++ [BB.toLazyByteString $ mpFunc mp x | mp <- msParams m]+        ++ [renderDuration dt, eta]+ msExec ::   Int ->   Item a ->@@ -109,27 +134,12 @@   Int ->   MonitorStdOut a ->   IO (Maybe BL.ByteString)-msExec i (Item x p l) ss st j m+msExec i it ss st j m   | i `mod` msPeriod m /= 0 = return Nothing-  | otherwise = do-    ct <- getCurrentTime-    let dt = ct `diffUTCTime` st-        -- Careful, don't evaluate this when i == ss.-        timePerIter = dt / fromIntegral (i - ss)-        -- -- Always 0; doesn't make much sense.-        -- tpi = if (i - ss) < 10-        --       then ""-        --       else renderDurationS timePerIter-        eta =-          if (i - ss) < 10-            then ""-            else renderDuration $ timePerIter * fromIntegral (j - i)-    return $-      Just $-        msRenderRow $-          [BL.pack (show i), renderLog p, renderLog l, renderLog (p * l)]-            ++ [BB.toLazyByteString $ mpFunc mp x | mp <- msParams m]-            ++ [renderDuration dt, eta]+  | i `mod` (msPeriod m * 100) == 0 = do+    l <- msDataLine i it ss st j m+    return $ Just $ msHeader m <> "\n" <> l+  | otherwise = Just <$> msDataLine i it ss st j m  -- | Monitor to a file; constructed with 'monitorFile'. data MonitorFile a = MonitorFile@@ -331,6 +341,7 @@   mapM_ mfHeader fs'   bs' <- mapM (mbOpen n frc) bs   mapM_ mbHeader bs'+  hSetBuffering stdout LineBuffering   return $ Monitor s fs' bs'  -- | Open the files associated with the 'Monitor' in append mode.@@ -340,10 +351,6 @@   bs' <- mapM (mbAppend n) bs   return $ Monitor s fs' bs' --- | Get header line of 'MonitorStdOut'.-mHeader :: Monitor a -> BL.ByteString-mHeader (Monitor s _ _) = msHeader s- -- | Execute monitors; print status information to files and return text to be -- printed to standard output and log file. mExec ::@@ -363,15 +370,11 @@   Monitor a ->   IO (Maybe BL.ByteString) mExec v i ss st xs j (Monitor s fs bs) = do-  -- XXX: I am not sure if this concurrency is necessary.-  mf <- async $ mapConcurrently_ (mfExec i $ headT xs) fs-  mb <- async $ mapConcurrently_ (mbExec i xs) bs-  ms <- async $ if v == Quiet-                then return Nothing-                else msExec i (headT xs) ss st j s-  wait mf-  wait mb-  wait ms+  mapM_ (mfExec i $ headT xs) fs+  mapM_ (mbExec i xs) bs+  if v == Quiet+    then return Nothing+    else msExec i (headT xs) ss st j s  -- | Close the files associated with the 'Monitor'. mClose :: Monitor a -> IO (Monitor a)

@@ -24,7 +24,6 @@  import qualified Data.ByteString.Builder as BB import qualified Data.Double.Conversion.ByteString as BC-import Lens.Micro  -- | Instruction about a parameter to monitor. data MonitorParameter a = MonitorParameter@@ -34,15 +33,15 @@     mpFunc :: a -> BB.Builder   } --- | Convert a parameter monitor from one data type to another using a lens.+-- | Convert a parameter monitor from one data type to another. -- -- For example, to monitor a 'Double' value being the first entry of a tuple: -- -- @--- mon = _1 @@ monitorDouble+-- mon = fst @. monitorDouble -- @-(@.) :: Lens' b a -> MonitorParameter a -> MonitorParameter b-(@.) l (MonitorParameter n f) = MonitorParameter n (\x -> f $ x ^. l)+(@.) :: (b -> a) -> MonitorParameter a -> MonitorParameter b+(@.) f (MonitorParameter n m) = MonitorParameter n (m . f)  -- | Monitor 'Int'. monitorInt ::

@@ -26,7 +26,6 @@  import qualified Data.ByteString.Builder as BB import qualified Data.Double.Conversion.ByteString as BC-import Lens.Micro  -- | Instruction about a parameter to monitor via batch means. Usually, the -- monitored parameter is average over the batch size. However, arbitrary@@ -35,24 +34,22 @@ -- -- XXX: Batch monitors are slow at the moment because the monitored parameter -- has to be extracted from the state for each iteration.-data MonitorParameterBatch a-  = MonitorParameterBatch-      { -- | Name of batch monitored parameter.-        mbpName :: String,-        -- | Instruction about how to extract the batch mean from the trace.-        mbpFunc :: [a] -> BB.Builder-      }+data MonitorParameterBatch a = MonitorParameterBatch+  { -- | Name of batch monitored parameter.+    mbpName :: String,+    -- | Instruction about how to extract the batch mean from the trace.+    mbpFunc :: [a] -> BB.Builder+  } --- | Convert a batch parameter monitor from one data type to another using a--- lens.+-- | Convert a batch parameter monitor from one data type to another. ----- For example, to batch monitor a real float value being the first entry of a tuple:+-- For example, to batch monitor the mean of the first entry of a tuple: -- -- @--- mon = _1 @# monitorBatchMeanRealFloat+-- mon = fst @# monitorBatchMean -- @-(@#) :: Lens' b a -> MonitorParameterBatch a -> MonitorParameterBatch b-(@#) l (MonitorParameterBatch n f) = MonitorParameterBatch n (f . map (^. l))+(@#) :: (b -> a) -> MonitorParameterBatch a -> MonitorParameterBatch b+(@#) f (MonitorParameterBatch n m) = MonitorParameterBatch n (m . map f)  mean :: Real a => [a] -> Double mean xs = realToFrac (sum xs) / fromIntegral (length xs)

@@ -19,8 +19,8 @@  import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL-import Mcmc.Internal.ByteString import Data.Time.Clock+import Mcmc.Internal.ByteString  -- | Adapted from System.ProgressBar.renderDuration of package -- [terminal-progressbar-0.4.1](https://hackage.haskell.org/package/terminal-progress-bar-0.4.1).

@@ -0,0 +1,108 @@+{-# LANGUAGE BangPatterns #-}++-- |+-- Module      :  Prior+-- Description :  Convenience functions to compute priors+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jul 23 13:26:14 2020.+module Mcmc.Prior+  ( -- * Continuous priors+    positive,+    negative,+    uniform,+    normal,+    exponential,+    gamma,++    -- * Discrete priors++    -- No discrete priors are available yet.++    -- * Auxiliary functions+    product',+  )+where++import Control.Monad+import Data.Maybe (fromMaybe)+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++-- | Improper uniform prior; larger than 0.+positive :: Double -> Log Double+positive x+  | x <= 0 = 0+  | otherwise = 1++-- | Improper uniform prior; lower than 0.+negative :: Double -> Log Double+negative x+  | x >= 0 = 0+  | otherwise = 1++-- | Uniform prior on [a, b].+uniform ::+  -- | Lower bound a.+  Double ->+  -- | Upper bound b.+  Double ->+  Double ->+  Log Double+uniform a b x+  | x <= a = 0+  | x >= b = 0+  | otherwise = Exp 0++-- | Normal distributed prior.+normal ::+  -- | Mean.+  Double ->+  -- | Standard deviation.+  Double ->+  Double ->+  Log Double+normal m s x = Exp $ S.logDensity d x+  where+    d = S.normalDistr m s++-- | Exponential distributed prior.+exponential ::+  -- | Rate.+  Double ->+  Double ->+  Log Double+exponential l x = Exp $ S.logDensity d x+  where+    d = S.exponential l++-- | Gamma distributed prior.+gamma ::+  -- | Shape.+  Double ->+  -- | Scale.+  Double ->+  Double ->+  Log Double+gamma k t x = Exp $ S.logDensity d x+  where+    d = S.gammaDistr k t++-- | Intelligent product that stops when encountering a zero.+--+-- Use with care because the elements have to be 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' = fromMaybe 0 . prodM++-- The type could be generalized to any MonadPlus Integer+prodM :: [Log Double] -> Maybe (Log Double)+prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0)) 1

@@ -4,11 +4,6 @@  -- TODO: Proposals on simplices: SimplexElementScale (?). --- TODO: Proposals on tree topologies.--- - NNI--- - Narrow (what is this, see RevBayes)--- - FNPR (dito)- -- | -- Module      :  Mcmc.Proposal -- Description :  Proposals and cycles@@ -34,7 +29,7 @@     Cycle (ccProposals),     fromList,     setOrder,-    getNCycles,+    getNIterations,     tuneCycle,     autotuneCycle,     summarizeCycle,@@ -77,7 +72,7 @@     -- | The weight determines how often a 'Proposal' is executed per iteration of     -- the Markov chain.     pWeight :: Int,-    -- | Simple proposal without tuning information.+    -- | Simple proposal without name, weight, and tuning information.     pSimple :: ProposalSimple a,     -- | Tuning is disabled if set to 'Nothing'.     pTuner :: Maybe (Tuner a)@@ -92,9 +87,6 @@ instance Ord (Proposal a) where   compare = compare `on` pName -convertP :: Lens' b a -> Proposal a -> Proposal b-convertP l (Proposal n w s t) = Proposal n w (convertS l s) (convertT l <$> t)- -- | Convert a proposal from one data type to another using a lens. -- -- For example:@@ -103,7 +95,7 @@ -- scaleFirstEntryOfTuple = scale >>> _1 -- @ (@~) :: Lens' b a -> Proposal a -> Proposal b-(@~) = convertP+(@~) l (Proposal n w s t) = Proposal n w (convertS l s) (convertT l <$> t)  -- | Simple proposal without tuning information. --@@ -167,10 +159,9 @@     let t' = max tuningParamMin (t * dt)     return $ m {pSimple = f t', pTuner = Just $ Tuner t' f} --- XXX: The desired acceptance ratio 0.44 is optimal for one-dimensional--- proposals; one could also store the affected number of dimensions with the--- proposal and tune towards an acceptance ratio accounting for the number of--- dimensions.+-- The desired acceptance ratio 0.44 is optimal for one-dimensional proposals;+-- one could also store the affected number of dimensions with the proposal and+-- tune towards an acceptance ratio accounting for the number of dimensions. -- -- The optimal ratios seem to be: -- - One dimension: 0.44 (numerical result).@@ -180,6 +171,14 @@ ratioOpt :: Double ratioOpt = 0.44 +-- Warn if acceptance ratio is lower.+ratioMin :: Double+ratioMin = 0.1++-- Warn if acceptance ratio is larger.+ratioMax :: Double+ratioMax = 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').@@ -206,11 +205,13 @@  instance Default Order where def = RandomO --- | 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'.+-- | 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'.+-- -- __Proposals must have unique names__, so that they can be identified. data Cycle a = Cycle   { ccProposals :: [Proposal a],@@ -233,8 +234,8 @@ setOrder o c = c {ccOrder = o}  -- | Replicate 'Proposal's according to their weights and possibly shuffle them.-getNCycles :: Cycle a -> Int -> GenIO -> IO [[Proposal a]]-getNCycles (Cycle xs o) n g = case o of+getNIterations :: Cycle a -> Int -> GenIO -> IO [[Proposal a]]+getNIterations (Cycle xs o) n g = case o of   RandomO -> shuffleN ps n g   SequentialO -> return $ replicate n ps   RandomReversibleO -> do@@ -269,8 +270,9 @@   BL.ByteString ->   BL.ByteString ->   BL.ByteString ->+  BL.ByteString ->   BL.ByteString-renderRow name weight nAccept nReject acceptRatio tuneParam = "   " <> nm <> wt <> na <> nr <> ra <> tp+renderRow name weight nAccept nReject acceptRatio tuneParam manualAdjustment = "   " <> nm <> wt <> na <> nr <> ra <> tp <> mt   where     nm = alignLeft 30 name     wt = alignRight 8 weight@@ -278,13 +280,22 @@     nr = alignRight 15 nReject     ra = alignRight 15 acceptRatio     tp = alignRight 20 tuneParam+    mt = alignRight 30 manualAdjustment  proposalHeader :: BL.ByteString proposalHeader =-  renderRow "Proposal" "Weight" "Accepted" "Rejected" "Ratio" "Tuning parameter"+  renderRow "Proposal" "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) weight nAccept nReject acceptRatio tuneParamStr+summarizeProposal m r =+  renderRow+    (BL.pack name)+    weight+    nAccept+    nReject+    acceptRatio+    tuneParamStr+    manualAdjustmentStr   where     name = pName m     weight = BB.toLazyByteString $ BB.intDec $ pWeight m@@ -292,17 +303,25 @@     nReject = BB.toLazyByteString $ maybe "" (BB.intDec . (^. _2)) r     acceptRatio = BL.fromStrict $ maybe "" (BC.toFixed 3 . (^. _3)) r     tuneParamStr = BL.fromStrict $ maybe "" (BC.toFixed 3) (tParam <$> pTuner m)+    check v+      | v < ratioMin = "ratio too low"+      | v > ratioMax = "ratio too high"+      | otherwise = ""+    manualAdjustmentStr = BL.fromStrict $ maybe "" (check . (^. _3)) r +hLine :: BL.ByteString -> BL.ByteString+hLine s = "   " <> BL.replicate (BL.length s - 3) '-'+ -- | Summarize the 'Proposal's in the 'Cycle'. Also report acceptance ratios. summarizeCycle :: Acceptance (Proposal a) -> Cycle a -> BL.ByteString summarizeCycle a c =   BL.intercalate "\n" $     [ "Summary of proposal(s) in cycle. " <> mpi <> " proposal(s) per iteration.",       proposalHeader,-      "   " <> BL.replicate (BL.length proposalHeader - 3) '-'+      hLine proposalHeader     ]       ++ [summarizeProposal m (ar m) | m <- ps]-      ++ ["   " <> BL.replicate (BL.length proposalHeader - 3) '-']+      ++ [hLine proposalHeader]   where     ps = ccProposals c     mpi = BB.toLazyByteString $ BB.intDec $ sum $ map pWeight ps

@@ -99,7 +99,7 @@ -- Hence, -- dx' = 1/(1-dx) - 1. fInv :: Double -> Double-fInv dx = recip (1-dx) - 1+fInv dx = recip (1 - dx) - 1  bactrianMult ::   Double ->

@@ -1,8 +1,3 @@--- Technically, only a Getter is needed when calculating the kernel of the proposal--- ('kernelCont', and similar functions). I tried splitting the lens into a getter--- and a setter. However, speed improvements were marginal, and some times not--- even measurable. Using a 'Lens'' is just easier, and has no real drawbacks.- -- | -- Module      :  Mcmc.Proposal.Generic -- Description :  Generic interface to create proposals@@ -15,8 +10,8 @@ -- -- Creation date: Thu May 14 20:26:27 2020. module Mcmc.Proposal.Generic-  ( proposalGenericContinuous,-    proposalGenericDiscrete,+  ( genericContinuous,+    genericDiscrete,   ) where @@ -45,7 +40,7 @@ {-# INLINEABLE sampleCont #-}  -- | Generic function to create proposals for continuous parameters ('Double').-proposalGenericContinuous ::+genericContinuous ::   (ContDistr d, ContGen d) =>   -- | Probability distribution   d ->@@ -55,7 +50,7 @@   -- required for biased proposals.   Maybe (Double -> Double) ->   ProposalSimple a-proposalGenericContinuous d f fInv = ProposalSimple $ sampleCont d f fInv+genericContinuous d f fInv = ProposalSimple $ sampleCont d f fInv  sampleDiscrete ::   (DiscreteDistr d, DiscreteGen d) =>@@ -77,7 +72,7 @@ {-# INLINEABLE sampleDiscrete #-}  -- | Generic function to create proposals for discrete parameters ('Int').-proposalGenericDiscrete ::+genericDiscrete ::   (DiscreteDistr d, DiscreteGen d) =>   -- | Probability distribution.   d ->@@ -87,4 +82,4 @@   -- required for biased proposals.   Maybe (Int -> Int) ->   ProposalSimple a-proposalGenericDiscrete fd f fInv = ProposalSimple $ sampleDiscrete fd f fInv+genericDiscrete fd f fInv = ProposalSimple $ sampleDiscrete fd f fInv

@@ -14,6 +14,7 @@ module Mcmc.Proposal.Scale   ( scale,     scaleUnbiased,+    scaleContrarily,   ) where @@ -24,7 +25,7 @@ -- The actual proposal with tuning parameter. The tuning parameter does not -- change the mean. scaleSimple :: Double -> Double -> Double -> ProposalSimple Double-scaleSimple k th t = proposalGenericContinuous (gammaDistr (k / t) (th * t)) (*) (Just recip)+scaleSimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) (*) (Just recip)  -- | Multiplicative proposal with Gamma distributed kernel. scale ::@@ -41,9 +42,10 @@   Proposal Double scale n w k th = createProposal n w (scaleSimple k th) --- | Multiplicative proposal with Gamma distributed kernel. The scale of the Gamma--- distributions is set to (shape)^{-1}, so that the mean of the Gamma--- distribution is 1.0.+-- | Multiplicative proposal with Gamma distributed kernel.+--+-- The scale of the Gamma distributions is set to (shape)^{-1}, so that the mean+-- of the Gamma distribution is 1.0. scaleUnbiased ::   -- | Name.   String ->@@ -55,3 +57,27 @@   Bool ->   Proposal Double scaleUnbiased n w k = createProposal n w (scaleSimple k (1 / k))++contra :: (Double, Double) -> Double -> (Double, Double)+contra (x, y) z = (x*z, y/z)++scaleContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)+scaleContrarilySimple k th t = genericContinuous (gammaDistr (k / t) (th * t)) contra (Just recip)++-- | Multiplicative proposal with Gamma distributed kernel.+--+-- The two values are scaled contrarily so that their product stays constant.+-- Contrary proposals are useful when parameters are confounded.+scaleContrarily ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Shape.+  Double ->+  -- | Scale.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal (Double, Double)+scaleContrarily n w k th = createProposal n w (scaleContrarilySimple k th)

@@ -15,6 +15,7 @@   ( slide,     slideSymmetric,     slideUniform,+    slideContrarily,   ) where @@ -25,7 +26,7 @@  -- The actual proposal with tuning parameter. slideSimple :: Double -> Double -> Double -> ProposalSimple Double-slideSimple m s t = proposalGenericContinuous (normalDistr m (s * t)) (+) (Just negate)+slideSimple m s t = genericContinuous (normalDistr m (s * t)) (+) (Just negate)  -- | Additive proposal with normally distributed kernel. slide ::@@ -44,7 +45,7 @@  -- The actual proposal with tuning parameter. slideSymmetricSimple :: Double -> Double -> ProposalSimple Double-slideSymmetricSimple s t = proposalGenericContinuous (normalDistr 0.0 (s * t)) (+) Nothing+slideSymmetricSimple s t = genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing  -- | Additive proposal with normally distributed kernel with mean zero. This -- proposal is very fast, because the Metropolis-Hastings ratio does not include@@ -64,7 +65,7 @@ -- The actual proposal with tuning parameter. slideUniformSimple :: Double -> Double -> ProposalSimple Double slideUniformSimple d t =-  proposalGenericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing+  genericContinuous (uniformDistr (- t * d) (t * d)) (+) Nothing  -- | Additive proposal with uniformly distributed kernel. This proposal is very fast, -- because the Metropolis-Hastings ratio does not include calculation of the@@ -80,3 +81,27 @@   Bool ->   Proposal Double slideUniform n w d = createProposal n w (slideUniformSimple d)++contra :: (Double, Double) -> Double -> (Double, Double)+contra (x, y) d = (x + d, y - d)++slideContrarilySimple :: Double -> Double -> Double -> ProposalSimple (Double, Double)+slideContrarilySimple m s t = genericContinuous (normalDistr m (s * t)) contra (Just negate)++-- | Additive proposal with normally distributed kernel.+--+-- The two values are slid contrarily so that their sum stays constant. Contrary+-- proposals are useful when parameters are confounded.+slideContrarily ::+  -- | Name.+  String ->+  -- | Weight.+  Int ->+  -- | Mean.+  Double ->+  -- | Standard deviation.+  Double ->+  -- | Enable tuning.+  Bool ->+  Proposal (Double, Double)+slideContrarily n w m s = createProposal n w (slideContrarilySimple m s)

@@ -26,10 +26,10 @@ import Control.Monad import Data.Aeson import Data.Aeson.TH-import Data.Maybe import qualified Data.ByteString.Lazy.Char8 as BL import Data.List hiding (cycle) import qualified Data.Map as M+import Data.Maybe import Data.Vector.Unboxed (Vector) import Data.Word -- TODO: Splitmix. Reproposal as soon as split mix is used and is available with the@@ -153,7 +153,7 @@ -- - cycle -- - monitor ----- To avoid incomplete continued runs, the 'mcmc' file is removed after load.+-- To avoid incomplete continued runs, the @.mcmc@ file is removed after load. loadStatus ::   FromJSON a =>   (a -> Log Double) ->

@@ -20,8 +20,8 @@  import Control.Monad import Control.Monad.ST-import qualified Data.Vector as V import Data.Vector (Vector)+import qualified Data.Vector as V import qualified Data.Vector.Mutable as M import System.Random.MWC   ( GenIO,@@ -49,9 +49,10 @@ -- 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)+  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

@@ -30,14 +30,14 @@  spec :: Spec spec =-  describe "getNCycles"-    $ it "returns the correct number of proposals in a cycle"-    $ do-      g <- create-      l1 <- length . head <$> getNCycles c 1 g-      l1 `shouldBe` 4-      l2 <- length . head <$> getNCycles (setOrder RandomReversibleO c) 1 g-      l2 `shouldBe` 8-      o3 <- head <$> getNCycles (setOrder SequentialReversibleO c) 1 g-      length o3 `shouldBe` 8-      o3 `shouldBe` [p1, p2, p2, p2, p2, p2, p2, p1]+  describe "getNIterations" $+    it "returns the correct number of proposals in a cycle" $+      do+        g <- create+        l1 <- length . head <$> getNIterations c 1 g+        l1 `shouldBe` 4+        l2 <- length . head <$> getNIterations (setOrder RandomReversibleO c) 1 g+        l2 `shouldBe` 8+        o3 <- head <$> getNIterations (setOrder SequentialReversibleO c) 1 g+        length o3 `shouldBe` 8+        o3 `shouldBe` [p1, p2, p2, p2, p2, p2, p2, p1]

@@ -61,27 +61,30 @@  spec :: Spec spec =-  describe "saveStatus and loadStatus"-    $ it "doesn't change the MCMC chain"-    $ do-      gen <- create-      let s =-            force $ quiet $ 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"-      r <- mh s-      r' <- mh s'-      -- Done during 'loadStatus'.-      -- removeFile "SaveSpec.json"-      item r `shouldBe` item r'-      iteration r `shouldBe` iteration r'-      trace r `shouldBe` trace r'-      g <- save $ generator r-      g' <- save $ generator r'-      g `shouldBe` g'+  describe "saveStatus and loadStatus" $+    it "doesn't change the MCMC chain" $+      do+        gen <- create+        let s =+              force $+                quiet $+                  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"+        r <- mh s+        r' <- mh s'+        -- Done during 'loadStatus'.+        -- removeFile "SaveSpec.json"+        item r `shouldBe` item r'+        iteration r `shouldBe` iteration r'+        trace r `shouldBe` trace r'+        g <- save $ generator r+        g' <- save $ generator r'+        g `shouldBe` g'+ -- -- TODO: Splitmix. This will only work with a splittable generator--- -- because getNCycles changes the generator.+-- -- because getNIterations changes the generator. -- describe "mhContinue" --   $ it "mh 200 + mhContinue 200 == mh 400" --   $ do