mcmc 0.2.1 → 0.2.2

raw patch · 12 files changed

+129/−95 lines, 12 filesdep +asyncPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: async

API changes (from Hackage documentation)

Files

• bench/Normal.hs +2/−2
• bench/Poisson.hs +1/−1
• mcmc.cabal +4/−7
• src/Mcmc.hs +1/−1
• src/Mcmc/Mcmc.hs +5/−4
• src/Mcmc/Metropolis.hs +14/−5
• src/Mcmc/Monitor.hs +11/−6
• src/Mcmc/Proposal.hs +8/−2
• src/Mcmc/Save.hs +9/−2
• src/Mcmc/Status.hs +62/−58
• src/Mcmc/Trace.hs +7/−2
• test/Mcmc/SaveSpec.hs +5/−5

@@ -62,7 +62,7 @@  normalBench :: GenIO -> IO () normalBench g = do-  let s = quiet $ noSave $ status "Normal" (const 1) lh proposals mon 0 nBurn nAutoTune nIter g+  let s = quiet $ status "Normal" (const 1) lh proposals mon 0 nBurn nAutoTune nIter g   void $ mh s  proposalsBactrian :: Cycle Double@@ -72,5 +72,5 @@  normalBactrianBench :: GenIO -> IO () normalBactrianBench g = do-  let s = quiet $ noSave $ status "NormalBactrian" (const 1) lh proposalsBactrian mon 0 nBurn nAutoTune nIter g+  let s = quiet $ status "NormalBactrian" (const 1) lh proposalsBactrian mon 0 nBurn nAutoTune nIter g   void $ mh s

@@ -82,5 +82,5 @@  poissonBench :: GenIO -> IO () poissonBench g = do-  let s = quiet $ noSave $ status "Poisson" (const 1) lh proposals mon initial nBurn nAutoTune nIter g+  let s = quiet $ status "Poisson" (const 1) lh proposals mon initial nBurn nAutoTune nIter g   void $ mh s

@@ -1,11 +1,6 @@-cabal-version: 2.2---- This file has been generated from package.yaml by hpack version 0.34.2.------ see: https://github.com/sol/hpack-+cabal-version:  2.2 name:           mcmc-version:        0.2.1+version:        0.2.2 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@@ -16,6 +11,7 @@ copyright:      Dominik Schrempf (2020) license:        GPL-3.0-or-later build-type:     Simple+ extra-source-files:     README.md     ChangeLog.md@@ -55,6 +51,7 @@   ghc-options: -Wall   build-depends:       aeson+    , async     , base >=4.7 && <5     , bytestring     , containers

@@ -110,7 +110,7 @@     -- space (see above) and to monitor the MCMC run, as well as some auxiliary     -- information.     status,-    noSave,+    saveWith,     force,     quiet,     debug,

@@ -196,12 +196,13 @@ mcmcSave :: ToJSON a => Mcmc a () mcmcSave = do   s <- get-  if save s-    then do-      mcmcInfoT "Save Markov chain. For long chains, this may take a while."+  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."       liftIO $ saveStatus (name s <> ".mcmc") s       mcmcInfoT "Done saving Markov chain."-    else mcmcInfoT "Do not save the Markov chain."+    Nothing -> mcmcInfoT "Do not save the Markov chain."  -- | Execute the 'Monitor's of the chain. See 'mExec'. mcmcMonitorExec :: ToJSON a => Mcmc a ()

@@ -91,12 +91,14 @@   | 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."@@ -110,16 +112,22 @@   | otherwise = do     mcmcInfoS $ "Burn in for " <> show b <> " cycles."     mcmcDebugS $ "Auto tuning period is " <> show t <> "."-    mcmcMonitorStdOutHeader     mhBurnInN b t     mcmcInfoT "Burn in finished."  -- Run for given number of iterations. mhRun :: ToJSON a => Int -> Mcmc a () mhRun n = do+  mcmcResetA   mcmcInfoS $ "Run chain for " <> show n <> " iterations."-  mcmcMonitorStdOutHeader-  mhNIter n+  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  mhT :: ToJSON a => Mcmc a () mhT = do@@ -149,8 +157,9 @@ mhContinue dn s   | dn <= 0 = error "mhContinue: The number of iterations is zero or negative."   | otherwise = mcmcRun (mhContinueT dn) s'-    where n' = iterations s + dn-          s' = s {iterations = n'}+  where+    n' = iterations s + dn+    s' = s {iterations = n'}  -- | Run a Markov chain for a given number of Metropolis-Hastings steps. mh ::

@@ -30,6 +30,7 @@   ) where +import Control.Concurrent.Async import Control.Monad import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL@@ -310,7 +311,7 @@           renderLog mlos :             [BB.toLazyByteString $ mbpFunc mbp (map state t) | mbp <- mbParams m]   where-    t = takeT (mbSize m) t'+    t = takeItems (mbSize m) t'     lps = map prior t     lls = map likelihood t     los = zipWith (*) lps lls@@ -362,11 +363,15 @@   Monitor a ->   IO (Maybe BL.ByteString) mExec v i ss st xs j (Monitor s fs bs) = do-  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+  -- 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  -- | Close the files associated with the 'Monitor'. mClose :: Monitor a -> IO (Monitor a)

@@ -168,9 +168,15 @@     return $ m {pSimple = f t', pTuner = Just $ Tuner t' f}  -- XXX: The desired acceptance ratio 0.44 is optimal for one-dimensional--- 'Proposal's; one could also store the affected number of dimensions with the--- 'Proposal' and tune towards an acceptance ratio accounting for the number of+-- 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).+-- - Five and more dimensions: 0.234 seems to be a good value (numerical result).+-- - Infinite dimensions: 0.234 (theorem for specific target distributions).+-- See Handbook of Markov chain Monte Carlo, chapter 4. ratioOpt :: Double ratioOpt = 0.44 

@@ -26,6 +26,7 @@ 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@@ -40,6 +41,7 @@ import Mcmc.Trace import Mcmc.Verbosity import Numeric.Log+import System.Directory import System.IO.Unsafe (unsafePerformIO) import System.Random.MWC import Prelude hiding (cycle)@@ -56,7 +58,7 @@       (Maybe Int) -- Auto tune.       Int -- Iterations.       Bool -- Force.-      Bool -- Save.+      (Maybe Int) -- Save.       Verbosity       (Vector Word32) -- Current seed. @@ -71,7 +73,7 @@     nm     it     i-    tr+    tr'     ac'     br     at@@ -82,6 +84,7 @@     g'     ts   where+    tr' = takeT (fromMaybe 0 sv) tr     ac' = transformKeysA (ccProposals c) [0 ..] ac     -- TODO: Splitmix. Remove as soon as split mix is used and is available with     -- the statistics package.@@ -142,12 +145,15 @@     c' = tuneCycle (M.mapMaybe id $ M.fromList $ zip (ccProposals c) ts) c  -- | Load a 'Status' from file.+-- -- Important information that cannot be saved and has to be provided again when -- a chain is restored: -- - prior function -- - likelihood function -- - cycle -- - monitor+--+-- To avoid incomplete continued runs, the 'mcmc' file is removed after load. loadStatus ::   FromJSON a =>   (a -> Log Double) ->@@ -177,4 +183,5 @@     ( error         "loadStatus: Provided likelihood function does not match the saved likelihood."     )+  removeFile fn   return s

@@ -28,7 +28,7 @@ module Mcmc.Status   ( Status (..),     status,-    noSave,+    saveWith,     force,     quiet,     debug,@@ -49,64 +49,68 @@  -- | The 'Status' contains all information to run an MCMC chain. It is -- constructed using the function 'status'.-data Status a-  = Status-      { -- MCMC related variables; saved.+data Status a = Status+  { -- MCMC related variables; saved. -        -- | The name of the MCMC chain; used as file prefix.-        name :: String,-        -- | The current 'Item' of the chain combines the current state and the-        -- current likelihood.-        item :: Item a,-        -- | The iteration is the number of completed cycles.-        iteration :: Int,-        -- | The 'Trace' of the Markov chain in reverse order, the most recent-        -- 'Item' is at the head of the list.-        trace :: Trace a,-        -- | For each 'Proposal', store the list of accepted (True) and rejected (False)-        -- proposals; for reasons of efficiency, the list is also stored in reverse-        -- order.-        acceptance :: Acceptance (Proposal a),-        -- | Number of burn in iterations; deactivate burn in with 'Nothing'.-        burnInIterations :: Maybe Int,-        -- | Auto tuning period (only during burn in); deactivate auto tuning with-        -- 'Nothing'.-        autoTuningPeriod :: Maybe Int,-        -- | Number of normal iterations excluding burn in. Note that auto tuning-        -- only happens during burn in.-        iterations :: Int,-        -- Auxiliary variables; saved.+    -- | The name of the MCMC chain; used as file prefix.+    name :: String,+    -- | The current 'Item' of the chain combines the current state and the+    -- current likelihood.+    item :: Item a,+    -- | The iteration is the number of completed cycles.+    iteration :: Int,+    -- | The 'Trace' of the Markov chain in reverse order, the most recent+    -- 'Item' is at the head of the list.+    trace :: Trace a,+    -- | For each 'Proposal', store the list of accepted (True) and rejected (False)+    -- proposals; for reasons of efficiency, the list is also stored in reverse+    -- order.+    acceptance :: Acceptance (Proposal a),+    -- | Number of burn in iterations; deactivate burn in with 'Nothing'.+    burnInIterations :: Maybe Int,+    -- | Auto tuning period (only during burn in); deactivate auto tuning with+    -- 'Nothing'.+    autoTuningPeriod :: Maybe Int,+    -- | Number of normal iterations excluding burn in. Note that auto tuning+    -- only happens during burn in.+    iterations :: Int,+    --+    -- Auxiliary variables; saved. -        -- | Overwrite output files? Default is 'False', change with 'force'.-        forceOverwrite :: Bool,-        -- | Save the chain at the end of the run? Default is 'True', change with 'noSave'.-        save :: Bool,-        -- | Verbosity.-        verbosity :: Verbosity,-        -- | The random number generator.-        generator :: GenIO,-        -- Auxiliary variables; not saved.+    -- | Overwrite output files? Default is 'False', change with 'force'.+    forceOverwrite :: Bool,+    -- | Save the chain with trace of given length at the end of the run?+    -- Default is no save ('Nothing'). Change with 'saveWith'.+    save :: Maybe Int,+    -- | Verbosity.+    verbosity :: Verbosity,+    -- | The random number generator.+    generator :: GenIO,+    --+    -- Auxiliary variables; not saved. -        -- | Starting time and starting iteration of chain; used to calculate-        -- run time and ETA.-        start :: Maybe (Int, UTCTime),-        -- | Handle to log file.-        logHandle :: Maybe Handle,-        -- Auxiliary functions; not saved.+    -- | Starting time and starting iteration of chain; used to calculate+    -- run time and ETA.+    start :: Maybe (Int, UTCTime),+    -- | Handle to log file.+    logHandle :: Maybe Handle,+    --+    -- Auxiliary functions; not saved. -        -- | The prior function. The un-normalized posterior is the product of the-        -- prior and the likelihood.-        priorF :: a -> Log Double,-        -- | The likelihood function. The un-normalized posterior is the product of-        -- the prior and the likelihood.-        likelihoodF :: a -> Log Double,-        -- Variables related to the algorithm; not saved.+    -- | The prior function. The un-normalized posterior is the product of the+    -- prior and the likelihood.+    priorF :: a -> Log Double,+    -- | The likelihood function. The un-normalized posterior is the product of+    -- the prior and the likelihood.+    likelihoodF :: a -> Log Double,+    --+    -- Variables related to the algorithm; not saved. -        -- | A set of 'Proposal's form a 'Cycle'.-        cycle :: Cycle a,-        -- | A 'Monitor' observing the chain.-        monitor :: Monitor a-      }+    -- | A set of 'Proposal's form a 'Cycle'.+    cycle :: Cycle a,+    -- | A 'Monitor' observing the chain.+    monitor :: Monitor a+  }  -- | Initialize the 'Status' of a Markov chain Monte Carlo run. status ::@@ -148,7 +152,7 @@       mT       nI       False-      True+      Nothing       Info       g       Nothing@@ -160,9 +164,9 @@   where     i = Item x (p x) (l x) --- | Do not save the Markov chain at the end.-noSave :: Status a -> Status a-noSave s = s {save = False}+-- | Save the Markov chain with trace of given length.+saveWith :: Int -> Status a -> Status a+saveWith n s = s {save = Just n}  -- | Overwrite existing files; it is not necessary to use 'force', when a chain -- is continued.

@@ -14,6 +14,7 @@     singletonT,     pushT,     headT,+    takeItems,     takeT,   ) where@@ -54,5 +55,9 @@ headT = head . fromTrace  -- | Get the N most recent items of the trace.-takeT :: Int -> Trace a -> [Item a]-takeT n = take n . fromTrace+takeItems :: Int -> Trace a -> [Item a]+takeItems n = take n . fromTrace++-- | Shorten the trace to given length.+takeT :: Int -> Trace a -> Trace a+takeT n = Trace . take n . fromTrace

@@ -23,7 +23,6 @@     stdDev,   ) import Statistics.Distribution.Normal-import System.Directory import System.Random.MWC import Test.Hspec @@ -67,13 +66,14 @@     $ do       gen <- create       let s =-            force $ quiet $ noSave $+            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'-      removeFile "SaveSpec.json"+      -- Done during 'loadStatus'.+      -- removeFile "SaveSpec.json"       item r `shouldBe` item r'       iteration r `shouldBe` iteration r'       trace r `shouldBe` trace r'@@ -86,10 +86,10 @@ --   $ it "mh 200 + mhContinue 200 == mh 400" --   $ do --     gen1 <- create---     let s1 = noSave $ status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 400 gen1+--     let s1 = status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 400 gen1 --     r1 <- mh s1 --     gen2 <- create---     let s2 = noSave $ status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 200 gen2+--     let s2 = status "SaveSpec" (const 1) likelihood proposals mon 0 nBurn nAutoTune 200 gen2 --     r2' <- mh s2 --     r2 <- mhContinue 200 r2' --     item r1 `shouldBe` item r2