mcmc 0.6.2.5 → 0.7.0.0

raw patch · 49 files changed

+2614/−1264 lines, 49 filesdep +addep +fixed-vectordep +randomdep −deepseqPVP ok

version bump matches the API change (PVP)

Dependencies added: ad, fixed-vector, random, splitmix

Dependencies removed: deepseq

API changes (from Hackage documentation)

Files

• ChangeLog.md +18/−0
• bench/Bench.hs +66/−44
• bench/Normal.hs +6/−6
• bench/Poisson.hs +74/−15
• mcmc.cabal +17/−8
• src/Mcmc.hs +12/−7
• src/Mcmc/Acceptance.hs +44/−14
• src/Mcmc/Algorithm.hs +7/−2
• src/Mcmc/Algorithm/MC3.hs +27/−18
• src/Mcmc/Algorithm/MHG.hs +95/−42
• src/Mcmc/Chain/Chain.hs +3/−3
• src/Mcmc/Chain/Link.hs +1/−1
• src/Mcmc/Chain/Save.hs +3/−4
• src/Mcmc/Chain/Trace.hs +8/−1
• src/Mcmc/Cycle.hs +46/−33
• src/Mcmc/Environment.hs +1/−1
• src/Mcmc/Internal/ByteString.hs +1/−1
• src/Mcmc/Internal/Gamma.hs +0/−205
• src/Mcmc/Internal/Random.hs +12/−26
• src/Mcmc/Internal/Shuffle.hs +8/−5
• src/Mcmc/Internal/SpecFunctions.hs +410/−0
• src/Mcmc/Jacobian.hs +32/−0
• src/Mcmc/Likelihood.hs +1/−1
• src/Mcmc/Logger.hs +1/−1
• src/Mcmc/MarginalLikelihood.hs +10/−12
• src/Mcmc/Mcmc.hs +24/−17
• src/Mcmc/Monitor.hs +1/−1
• src/Mcmc/Monitor/Log.hs +1/−1
• src/Mcmc/Monitor/Parameter.hs +1/−1
• src/Mcmc/Monitor/ParameterBatch.hs +1/−1
• src/Mcmc/Monitor/Time.hs +1/−1
• src/Mcmc/Posterior.hs +1/−1
• src/Mcmc/Prior.hs +26/−22
• src/Mcmc/Proposal.hs +108/−88
• src/Mcmc/Proposal/Bactrian.hs +21/−22
• src/Mcmc/Proposal/Generic.hs +6/−6
• src/Mcmc/Proposal/Hamiltonian.hs +0/−601
• src/Mcmc/Proposal/Hamiltonian/Common.hs +215/−0
• src/Mcmc/Proposal/Hamiltonian/Hamiltonian.hs +210/−0
• src/Mcmc/Proposal/Hamiltonian/Internal.hs +413/−0
• src/Mcmc/Proposal/Hamiltonian/Masses.hs +329/−0
• src/Mcmc/Proposal/Hamiltonian/Nuts.hs +305/−0
• src/Mcmc/Proposal/Scale.hs +9/−9
• src/Mcmc/Proposal/Simplex.hs +10/−15
• src/Mcmc/Proposal/Slide.hs +14/−14
• src/Mcmc/Settings.hs +1/−1
• src/Mcmc/Statistics/Types.hs +1/−1
• test/Mcmc/ProposalSpec.hs +3/−3
• test/Mcmc/SaveSpec.hs +10/−9

@@ -4,6 +4,24 @@  ## Unreleased changes ++## 0.7.0.0++-   Use random-1.2. This means, that the random number generator type has changed+    from `GenIO` (`System.Random.MWC`) to `IOGenM StdGen`+    (`System.Random.Stateful`).+-   Always try to save chain, not only on `UserInterupt`.+-   Hamiltonian proposals:+    -   Split `HSettings` into `HParams` and `HStructure` (preparation for NUTS).+    -   Rename `HTune` to `HTuningConf`.+    -   No-U-Turn sampler (NUTS).+-   Significant changes to `Proposal` and `Tuner`.+    -   `PResult`.+    -   `ProposalSimple` -> `PFunction`; Return type is `(PResult, Maybe AcceptanceCounts)`.+++## 0.6.2.5+ -   Provide unsafe loader for MHG algorithm (useful for initializing chains with     different prior or likelihood functions from saves). -   Documentation and readme.

@@ -1,6 +1,6 @@ -- | -- Module      :  Main--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -15,11 +15,11 @@  import Criterion.Main import Data.List-import Mcmc.Internal.Gamma+import Mcmc.Internal.SpecFunctions import Normal import Numeric.SpecFunctions import Poisson-import System.Random.MWC+import System.Random.Stateful  gammaBenchG :: RealFloat a => (a -> a) -> [a] -> a gammaBenchG f = foldl' (\acc x -> acc + (f x)) 0@@ -30,7 +30,7 @@  main :: IO () main = do-  g <- create+  g <- newIOGenM $ mkStdGen 0   defaultMain     [ bgroup         "Normal"@@ -38,8 +38,12 @@           bench "Bactrian" $ nfIO (normalBactrianBench g),           bench "LargeCycle" $ nfIO (normalLargeCycleBench g)         ],-      bench "Poisson" $ nfIO (poissonBench g),       bgroup+        "Poisson"+        [ bench "Random walk" $ nfIO (poissonBench g),+          bench "Hamiltonian" $ nfIO (poissonHamiltonianBench g)+        ],+      bgroup         "MC3"         [ bench "MC3 2" $ nfIO (normalMC3 g 2),           bench "MC3 3" $ nfIO (normalMC3 g 3),@@ -57,58 +61,76 @@     ]  -- benchmarking Normal/Slide--- time                 42.31 ms   (41.88 ms .. 42.60 ms)---                      1.000 R²   (0.999 R² .. 1.000 R²)--- mean                 42.75 ms   (42.52 ms .. 43.29 ms)--- std dev              661.4 μs   (347.4 μs .. 1.074 ms)+-- time                 24.12 ms   (23.52 ms .. 24.76 ms)+--                      0.997 R²   (0.995 R² .. 0.999 R²)+-- mean                 23.70 ms   (23.44 ms .. 24.02 ms)+-- std dev              658.2 μs   (448.5 μs .. 837.2 μs)  -- benchmarking Normal/Bactrian--- time                 45.51 ms   (45.30 ms .. 45.92 ms)---                      1.000 R²   (0.999 R² .. 1.000 R²)--- mean                 45.41 ms   (45.31 ms .. 45.61 ms)--- std dev              276.0 μs   (141.5 μs .. 460.8 μs)+-- time                 26.65 ms   (26.08 ms .. 27.26 ms)+--                      0.998 R²   (0.996 R² .. 1.000 R²)+-- mean                 26.53 ms   (26.32 ms .. 26.91 ms)+-- std dev              626.9 μs   (322.3 μs .. 1.040 ms)  -- benchmarking Normal/LargeCycle--- time                 68.82 ms   (67.18 ms .. 70.81 ms)---                      0.999 R²   (0.997 R² .. 1.000 R²)--- mean                 67.68 ms   (67.26 ms .. 68.59 ms)--- std dev              1.074 ms   (618.2 μs .. 1.602 ms)+-- time                 52.48 ms   (51.51 ms .. 53.83 ms)+--                      0.999 R²   (0.998 R² .. 1.000 R²)+-- mean                 52.54 ms   (51.62 ms .. 54.90 ms)+-- std dev              2.765 ms   (758.0 μs .. 4.788 ms)+-- variance introduced by outliers: 14% (moderately inflated) --- benchmarking Poisson--- time                 72.94 ms   (63.07 ms .. 87.73 ms)---                      0.953 R²   (0.920 R² .. 1.000 R²)--- mean                 64.76 ms   (62.84 ms .. 71.62 ms)--- std dev              5.785 ms   (783.9 μs .. 10.08 ms)--- variance introduced by outliers: 26% (moderately inflated)+-- benchmarking Poisson/Random walk+-- time                 86.57 ms   (86.06 ms .. 86.97 ms)+--                      1.000 R²   (1.000 R² .. 1.000 R²)+-- mean                 86.57 ms   (86.20 ms .. 87.14 ms)+-- std dev              775.9 μs   (415.4 μs .. 1.205 ms) +-- benchmarking Poisson/Hamiltonian+-- time                 2.890 s    (2.313 s .. 3.359 s)+--                      0.993 R²   (0.992 R² .. 1.000 R²)+-- mean                 3.082 s    (2.960 s .. 3.153 s)+-- std dev              120.4 ms   (39.93 ms .. 164.0 ms)+-- variance introduced by outliers: 19% (moderately inflated)+ -- benchmarking MC3/MC3 2--- time                 13.08 ms   (12.73 ms .. 13.44 ms)---                      0.993 R²   (0.986 R² .. 0.997 R²)--- mean                 13.41 ms   (13.16 ms .. 13.72 ms)--- std dev              682.7 μs   (520.3 μs .. 874.0 μs)--- variance introduced by outliers: 22% (moderately inflated)+-- time                 7.771 ms   (7.573 ms .. 8.008 ms)+--                      0.996 R²   (0.993 R² .. 0.999 R²)+-- mean                 7.810 ms   (7.722 ms .. 7.974 ms)+-- std dev              328.3 μs   (223.0 μs .. 517.6 μs)+-- variance introduced by outliers: 19% (moderately inflated)  -- benchmarking MC3/MC3 3--- time                 19.19 ms   (18.86 ms .. 19.59 ms)---                      0.998 R²   (0.996 R² .. 1.000 R²)--- mean                 19.28 ms   (19.11 ms .. 19.51 ms)--- std dev              454.1 μs   (339.1 μs .. 608.4 μs)+-- time                 10.69 ms   (10.57 ms .. 10.87 ms)+--                      0.997 R²   (0.993 R² .. 1.000 R²)+-- mean                 10.73 ms   (10.65 ms .. 10.90 ms)+-- std dev              308.1 μs   (132.2 μs .. 528.9 μs)  -- benchmarking MC3/MC3 4--- time                 25.01 ms   (24.21 ms .. 25.66 ms)---                      0.997 R²   (0.996 R² .. 0.999 R²)--- mean                 24.21 ms   (23.99 ms .. 24.55 ms)--- std dev              606.1 μs   (414.7 μs .. 738.5 μs)+-- time                 13.81 ms   (13.66 ms .. 14.08 ms)+--                      0.999 R²   (0.997 R² .. 1.000 R²)+-- mean                 13.73 ms   (13.67 ms .. 13.85 ms)+-- std dev              213.2 μs   (95.38 μs .. 334.4 μs)  -- benchmarking MC3/MC3 5--- time                 28.39 ms   (26.99 ms .. 29.34 ms)---                      0.995 R²   (0.990 R² .. 0.999 R²)--- mean                 31.13 ms   (30.09 ms .. 33.49 ms)--- std dev              3.009 ms   (984.0 μs .. 4.319 ms)--- variance introduced by outliers: 40% (moderately inflated)+-- time                 16.80 ms   (16.65 ms .. 16.93 ms)+--                      0.999 R²   (0.997 R² .. 1.000 R²)+-- mean                 16.99 ms   (16.87 ms .. 17.23 ms)+-- std dev              376.6 μs   (162.3 μs .. 589.3 μs)  -- benchmarking MC3/MC3 10--- time                 57.25 ms   (56.98 ms .. 57.61 ms)+-- time                 32.59 ms   (32.10 ms .. 33.20 ms)+--                      0.999 R²   (0.998 R² .. 1.000 R²)+-- mean                 32.21 ms   (32.07 ms .. 32.46 ms)+-- std dev              390.9 μs   (219.6 μs .. 640.0 μs)++-- benchmarking probability/gamma function general+-- time                 60.73 ms   (60.37 ms .. 61.26 ms) --                      1.000 R²   (1.000 R² .. 1.000 R²)--- mean                 57.46 ms   (57.34 ms .. 57.56 ms)--- std dev              192.4 μs   (140.1 μs .. 284.0 μs)+-- mean                 60.70 ms   (60.49 ms .. 60.91 ms)+-- std dev              385.5 μs   (287.3 μs .. 513.5 μs)++-- benchmarking probability/gamma function specialized+-- time                 59.00 ms   (58.74 ms .. 59.32 ms)+--                      1.000 R²   (1.000 R² .. 1.000 R²)+-- mean                 58.83 ms   (58.76 ms .. 58.97 ms)+-- std dev              168.3 μs   (105.7 μs .. 245.7 μs)

@@ -1,7 +1,7 @@ -- | -- Module      :  Normal -- Description :  Benchmark Metropolis-Hastings-Green algorithm--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -19,7 +19,7 @@  import Control.Monad import Mcmc-import System.Random.MWC+import System.Random.Stateful  trueMean :: Double trueMean = 5@@ -42,7 +42,7 @@ mon :: Monitor Double mon = Monitor monStd [] [] -normalSlideBench :: GenIO -> IO ()+normalSlideBench :: IOGenM StdGen -> IO () normalSlideBench g = do   let s =         Settings@@ -64,7 +64,7 @@     [slideSymmetric 1.0 (PName $ "Medium " ++ show i) (pWeight 1) Tune | i <- [0 .. 100 :: Int]]  -- Should have the same run time as 'normalSlide'.-normalLargeCycleBench :: GenIO -> IO ()+normalLargeCycleBench :: IOGenM StdGen -> IO () normalLargeCycleBench g = do   let s =         Settings@@ -83,7 +83,7 @@ ccBactrian :: Cycle Double ccBactrian = cycleFromList [slideBactrian 0.5 1.0 (PName "Bactrian") (pWeight 1) Tune] -normalBactrianBench :: GenIO -> IO ()+normalBactrianBench :: IOGenM StdGen -> IO () normalBactrianBench g = do   let s =         Settings@@ -99,7 +99,7 @@   a <- mhg s noPrior lh ccBactrian mon 0 g   void $ mcmc s a -normalMC3 :: GenIO -> Int -> IO ()+normalMC3 :: IOGenM StdGen -> Int -> IO () normalMC3 g n = do   let mcmcS =         Settings

@@ -1,7 +1,14 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_GHC -Wno-orphans #-}+ -- | -- Module      :  Poisson -- Description :  Poisson regression model for airline fatalities--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -13,16 +20,26 @@ -- See https://revbayes.github.io/tutorials/mcmc/poisson.html. module Poisson   ( poissonBench,+    poissonHamiltonianBench,   ) where  import Control.Monad-import Lens.Micro+import Data.Aeson+import Data.Typeable+import qualified Data.Vector.Fixed as VB+import qualified Data.Vector.Fixed.Boxed as VB+import qualified Data.Vector.Storable as VS import Mcmc-import System.Random.MWC+import System.Random.Stateful -type I = (Double, Double)+type IG = VB.Vec2 +type I = IG Double++instance (VB.Arity n, ToJSON a) => ToJSON (VB.Vec n a) where+  toJSON = toJSON . VB.toList+ fatalities :: [Int] fatalities = [24, 25, 31, 31, 22, 21, 26, 20, 16, 22] @@ -32,43 +49,48 @@     ys = [1976.0 .. 1985.0]     m = sum ys / fromIntegral (length ys) -f :: Int -> Double -> I -> Likelihood-f ft yr (a, b) = poisson l (fromIntegral ft)+f :: (RealFloat a, Typeable a) => Int -> Double -> IG a -> LikelihoodG a+f ft yr xs = poisson l ft   where-    l = exp $ a + b * yr+    a = xs VB.! 0+    b = xs VB.! 1+    l = exp $ a + b * realToFrac yr -lh :: LikelihoodFunction I+lh :: (RealFloat a, Typeable a) => LikelihoodFunctionG (IG a) a lh x = product [f ft yr x | (ft, yr) <- zip fatalities normalizedYears]  proposalAlpha :: Proposal I-proposalAlpha = _1 @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune+proposalAlpha = (VB.element 0) @~ slideSymmetric 0.2 (PName "Alpha") (pWeight 1) NoTune  proposalBeta :: Proposal I-proposalBeta = _2 @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune+proposalBeta = (VB.element 1) @~ slideSymmetric 0.2 (PName "Beta") (pWeight 1) NoTune  proposals :: Cycle I proposals = cycleFromList [proposalAlpha, proposalBeta]  initial :: I-initial = (0, 0)+initial = VB.mk2 0 0  monAlpha :: MonitorParameter I-monAlpha = fst >$< monitorDouble "alpha"+monAlpha = (VB.! 0) >$< monitorDouble "alpha"  monBeta :: MonitorParameter I-monBeta = snd >$< monitorDouble "beta"+monBeta = (VB.! 1) >$< monitorDouble "beta"  monStd :: MonitorStdOut I monStd = monitorStdOut [monAlpha, monBeta] 150 +-- monF :: MonitorFile I+-- monF = monitorFile "params" [monAlpha, monBeta] 10+ mon :: Monitor I mon = Monitor monStd [] [] -poissonBench :: GenIO -> IO ()+poissonBench :: IOGenM StdGen -> IO () poissonBench g = do   let s =         Settings-          (AnalysisName "Poisson")+          (AnalysisName "PoissonR")           (BurnInWithAutoTuning 2000 200)           (Iterations 10000)           TraceAuto@@ -78,4 +100,41 @@           LogStdOutOnly           Quiet   a <- mhg s noPrior lh proposals mon initial g+  void $ mcmc s a++toVec :: I -> VS.Vector Double+toVec xs = VS.generate 2 (\i -> xs VB.! i)++fromVec :: I -> VS.Vector Double -> I+fromVec _ xs = VB.mk2 (xs VS.! 0) (xs VS.! 1)++hparams :: HParams+hparams = HParams Nothing Nothing Nothing++htconf :: HTuningConf+htconf = HTuningConf HTuneLeapfrog HTuneAllMasses++hstruct :: HStructure IG+hstruct = HStructure initial toVec fromVec++target :: HTarget IG+target = HTarget Nothing lh Nothing++hmcProposal :: Cycle I+hmcProposal = cycleFromList [hamiltonian hparams htconf hstruct target (PName "Hamiltonian") (pWeight 1)]++poissonHamiltonianBench :: IOGenM StdGen -> IO ()+poissonHamiltonianBench g = do+  let s =+        Settings+          (AnalysisName "PoissonH")+          (BurnInWithAutoTuning 2000 200)+          (Iterations 10000)+          TraceAuto+          Overwrite+          Sequential+          NoSave+          LogStdOutOnly+          Quiet+  a <- mhg s noPrior lh hmcProposal mon initial g   void $ mcmc s a

@@ -1,6 +1,6 @@ cabal-version:      3.0 name:               mcmc-version:            0.6.2.5+version:            0.7.0.0 synopsis:           Sample from a posterior using Markov chain Monte Carlo description:   Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme>@@ -10,7 +10,7 @@ bug-reports:        https://github.com/dschrempf/mcmc/issues author:             Dominik Schrempf maintainer:         dominik.schrempf@gmail.com-copyright:          Dominik Schrempf (2021)+copyright:          2021 Dominik Schrempf license:            GPL-3.0-or-later build-type:         Simple extra-source-files:@@ -53,7 +53,8 @@     Mcmc.Chain.Trace     Mcmc.Cycle     Mcmc.Environment-    Mcmc.Internal.Gamma+    Mcmc.Internal.SpecFunctions+    Mcmc.Jacobian     Mcmc.Likelihood     Mcmc.Logger     Mcmc.MarginalLikelihood@@ -68,7 +69,9 @@     Mcmc.Proposal     Mcmc.Proposal.Bactrian     Mcmc.Proposal.Generic-    Mcmc.Proposal.Hamiltonian+    Mcmc.Proposal.Hamiltonian.Common+    Mcmc.Proposal.Hamiltonian.Hamiltonian+    Mcmc.Proposal.Hamiltonian.Nuts     Mcmc.Proposal.Scale     Mcmc.Proposal.Simplex     Mcmc.Proposal.Slide@@ -79,11 +82,14 @@     Mcmc.Internal.ByteString     Mcmc.Internal.Random     Mcmc.Internal.Shuffle+    Mcmc.Proposal.Hamiltonian.Internal+    Mcmc.Proposal.Hamiltonian.Masses     Paths_mcmc    autogen-modules: Paths_mcmc   hs-source-dirs:  src   build-depends:+    , ad     , aeson     , async     , base               >=4.7 && <5@@ -92,7 +98,6 @@     , containers     , covariance         >=0.2     , data-default-    , deepseq     , directory     , dirichlet     , double-conversion@@ -104,6 +109,8 @@     , parallel     , pretty-show     , primitive+    , random+    , splitmix     , statistics     , time     , transformers@@ -124,7 +131,7 @@     , base        >=4.7 && <5     , hspec     , mcmc-    , mwc-random+    , random     , statistics  benchmark mcmc-bench@@ -139,9 +146,11 @@   autogen-modules: Paths_mcmc   hs-source-dirs:  bench   build-depends:+    , aeson     , base            >=4.7 && <5     , criterion+    , fixed-vector     , math-functions     , mcmc-    , microlens-    , mwc-random+    , random+    , vector

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc -- Description :  Markov chain Monte Carlo samplers, batteries included--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -110,8 +110,8 @@     -- tuning parameter @t@ is implemented as     --     -- @-    -- slideSimple :: Double -> Double -> Double -> ProposalSimple Double-    -- slideSimple m s t =+    -- slide :: Double -> Double -> Double -> PFunction Double+    -- slide m s t =     --   genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing     -- @     --@@ -129,7 +129,6 @@     PWeight (..),     pWeight,     Proposal,-    JacobianFunction,     (@~),     liftProposal,     liftProposalWith,@@ -143,8 +142,10 @@     slideUniformSymmetric,     slideContrarily,     slideBactrian,-    module Mcmc.Proposal.Hamiltonian,     module Mcmc.Proposal.Simplex,+    module Mcmc.Proposal.Hamiltonian.Common,+    module Mcmc.Proposal.Hamiltonian.Hamiltonian,+    module Mcmc.Proposal.Hamiltonian.Nuts,      -- ** Cycles     Cycle,@@ -181,9 +182,10 @@     module Mcmc.Monitor.Parameter,     module Mcmc.Monitor.ParameterBatch, -    -- * Priors, likelihoodl, and posteriors+    -- * Priors, likelihoods, Jacobians, and posteriors     module Mcmc.Prior,     module Mcmc.Likelihood,+    module Mcmc.Jacobian,     module Mcmc.Posterior,      -- * Run MCMC samplers@@ -209,6 +211,7 @@ import Mcmc.Algorithm.MC3 import Mcmc.Algorithm.MHG import Mcmc.Cycle+import Mcmc.Jacobian import Mcmc.Likelihood import Mcmc.MarginalLikelihood import Mcmc.Mcmc@@ -219,7 +222,9 @@ import Mcmc.Prior import Mcmc.Proposal import Mcmc.Proposal.Bactrian-import Mcmc.Proposal.Hamiltonian+import Mcmc.Proposal.Hamiltonian.Common+import Mcmc.Proposal.Hamiltonian.Hamiltonian+import Mcmc.Proposal.Hamiltonian.Nuts import Mcmc.Proposal.Scale import Mcmc.Proposal.Simplex import Mcmc.Proposal.Slide

@@ -1,7 +1,9 @@+{-# LANGUAGE TemplateHaskell #-}+ -- | -- Module      :  Mcmc.Acceptance -- Description :  Handle acceptance rates--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -12,9 +14,12 @@ module Mcmc.Acceptance   ( -- * Acceptance rates     AcceptanceRate,+    AcceptanceCounts (..),     Acceptance (fromAcceptance),     emptyA,-    pushA,+    pushAccept,+    pushReject,+    pushAcceptanceCounts,     resetA,     transformKeysA,     acceptanceRate,@@ -22,18 +27,36 @@   ) where -import Control.DeepSeq import Data.Aeson-import Data.Bifunctor+import Data.Aeson.TH import Data.Foldable import qualified Data.Map.Strict as M  -- | Acceptance rate. type AcceptanceRate = Double +-- | Number of accepted and rejected proposals.+data AcceptanceCounts = AcceptanceCounts+  { nAccepted :: !Int,+    nRejected :: !Int+  }+  deriving (Show, Eq, Ord)++$(deriveJSON defaultOptions ''AcceptanceCounts)++addAccept :: AcceptanceCounts -> AcceptanceCounts+addAccept (AcceptanceCounts a r) = AcceptanceCounts (a + 1) r++addReject :: AcceptanceCounts -> AcceptanceCounts+addReject (AcceptanceCounts a r) = AcceptanceCounts a (r + 1)++addAcceptanceCounts :: AcceptanceCounts -> AcceptanceCounts -> AcceptanceCounts+addAcceptanceCounts (AcceptanceCounts al rl) (AcceptanceCounts ar rr) =+  AcceptanceCounts (al + ar) (rl + rr)+ -- | For each key @k@, store the number of accepted and rejected proposals.-newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k (Int, Int)}-  deriving (Eq, Read, Show)+newtype Acceptance k = Acceptance {fromAcceptance :: M.Map k AcceptanceCounts}+  deriving (Eq, Show)  instance ToJSONKey k => ToJSON (Acceptance k) where   toJSON (Acceptance m) = toJSON m@@ -46,13 +69,20 @@ -- -- Initialize an empty storage of accepted/rejected values. emptyA :: Ord k => [k] -> Acceptance k-emptyA ks = Acceptance $ M.fromList [(k, (0, 0)) | k <- ks]+emptyA ks = Acceptance $ M.fromList [(k, AcceptanceCounts 0 0) | k <- ks] --- | For key @k@, prepend an accepted (True) or rejected (False) proposal.-pushA :: Ord k => k -> Bool -> Acceptance k -> Acceptance k-pushA k True = Acceptance . M.adjust (force . first succ) k . fromAcceptance-pushA k False = Acceptance . M.adjust (force . second succ) k . fromAcceptance+-- | For key @k@, add an accept.+pushAccept :: Ord k => k -> Acceptance k -> Acceptance k+pushAccept k = Acceptance . M.adjust addAccept k . fromAcceptance +-- | For key @k@, add a reject.+pushReject :: Ord k => k -> Acceptance k -> Acceptance k+pushReject k = Acceptance . M.adjust addReject k . fromAcceptance++-- | For key @k@, add acceptance counts.+pushAcceptanceCounts :: Ord k => k -> AcceptanceCounts -> Acceptance k -> Acceptance k+pushAcceptanceCounts k c = Acceptance . M.adjust (addAcceptanceCounts c) k . fromAcceptance+ -- | Reset acceptance storage. resetA :: Ord k => Acceptance k -> Acceptance k resetA = emptyA . M.keys . fromAcceptance@@ -73,8 +103,8 @@ -- zero). acceptanceRate :: Ord k => k -> Acceptance k -> Maybe (Int, Int, AcceptanceRate) acceptanceRate k a = case fromAcceptance a M.!? k of-  Just (0, 0) -> Nothing-  Just (as, rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))+  Just (AcceptanceCounts 0 0) -> Nothing+  Just (AcceptanceCounts as rs) -> Just (as, rs, fromIntegral as / fromIntegral (as + rs))   Nothing -> error "acceptanceRate: Key not found in map."  -- | Acceptance rates for all proposals.@@ -84,7 +114,7 @@ acceptanceRates :: Acceptance k -> M.Map k (Maybe AcceptanceRate) acceptanceRates =   M.map-    ( \(as, rs) ->+    ( \(AcceptanceCounts as rs) ->         if as + rs == 0           then Nothing           else Just $ fromIntegral as / fromIntegral (as + rs)

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Algorithm -- Description :  Algortihms for Markov chain Monte Carlo samplers--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -17,6 +17,7 @@ import qualified Data.ByteString.Lazy.Char8 as BL import Data.Time import Mcmc.Cycle+import Mcmc.Proposal import Mcmc.Settings  -- | Class for algorithms used by MCMC samplers.@@ -38,10 +39,14 @@   --   -- NOTE: Computation in the 'IO' Monad is necessary because the trace is   -- mutable.-  aAutoTune :: Int -> a -> IO a+  aAutoTune :: TuningType -> Int -> a -> IO a    -- | Reset acceptance counts.   aResetAcceptance :: a -> a++  -- | Clean after burn in. In particular, this is used to reduce the length of+  -- the trace, if required.+  aCleanAfterBurnIn :: TraceLength -> a -> IO a    -- | Summarize the cycle.   aSummarizeCycle :: IterationMode -> a -> BL.ByteString

@@ -5,7 +5,7 @@ -- | -- Module      :  Mcmc.Algorithm.MC3 -- Description :  Metropolis-coupled Markov chain Monte Carlo algorithm--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -79,9 +79,7 @@ import Mcmc.Proposal import Mcmc.Settings import Numeric.Log hiding (sum)-import System.Random.MWC---- import Debug.Trace hiding (trace)+import System.Random.Stateful  -- | Total number of parallel chains. --@@ -130,9 +128,9 @@     savedMC3ReciprocalTemperatures :: ReciprocalTemperatures,     savedMC3Iteration :: Int,     savedMC3SwapAcceptance :: Acceptance Int,-    savedMC3Generator :: U.Vector Word32+    savedMC3Generator :: (Word64, Word64)   }-  deriving (Eq, Read, Show)+  deriving (Eq, Show)  $(deriveJSON defaultOptions ''SavedMC3) @@ -175,7 +173,7 @@     mc3Iteration :: Int,     -- | Number of accepted and rejected swaps.     mc3SwapAcceptance :: Acceptance Int,-    mc3Generator :: GenIO+    mc3Generator :: IOGenM StdGen   }  instance ToJSON a => Algorithm (MC3 a) where@@ -185,6 +183,7 @@   aIterate = mc3Iterate   aAutoTune = mc3AutoTune   aResetAcceptance = mc3ResetAcceptance+  aCleanAfterBurnIn = mc3CleanAfterBurnIn   aSummarizeCycle = mc3SummarizeCycle   aOpenMonitors = mc3OpenMonitors   aExecuteMonitors = mc3ExecuteMonitors@@ -266,8 +265,9 @@     l = link c     t = trace c --- TODO: Splitmix. Initialization of the MC3 algorithm is an IO action because--- the generators have to be split. And also because of the mutable trace.+-- TODO (medium): Splitmix. Initialization of the MC3 algorithm is an IO action+-- because the generators have to be split. And also because of the mutable+-- trace.  -- | Initialize an MC3 algorithm with a given number of chains. --@@ -284,7 +284,7 @@   Cycle a ->   Monitor a ->   InitialState a ->-  GenIO ->+  IOGenM StdGen ->   IO (MC3 a) mc3 sMc3 s pr lh cc mn i0 g   | n < 2 = error "mc3: The number of chains must be two or larger."@@ -292,7 +292,8 @@   | sn < 1 || sn > n - 1 = error "mc3: The number of swaps must be in [1, NChains - 1]."   | otherwise = do       -- Split random number generators.-      gs <- V.fromList <$> splitGen n g+      rs <- replicateM n $ splitGenM g+      gs <- V.fromList <$> mapM newIOGenM rs       cs <- V.mapM (mhg s pr lh cc mn i0) gs       hcs <- V.izipWithM (initMHG pr lh) (V.convert bs) cs       return $ MC3 sMc3 hcs bs 0 (emptyA [0 .. n - 2]) g@@ -407,10 +408,10 @@       -- traceIO $ "Log priors (left, right, after swap): " <> show (ln prL') <> " " <> show (ln prR')       -- traceIO $ "Log likelihoods (left, right, before swap): " <> show (ln lhL) <> " " <> show (ln lhR)       -- traceIO $ "Log likelihood (left, right, after swap): " <> show (ln lhL') <> " " <> show (ln lhR')-      let !ac' = pushA i True (mc3SwapAcceptance a)+      let !ac' = pushAccept i (mc3SwapAcceptance a)       return $ a {mc3MHGChains = y, mc3SwapAcceptance = ac'}     else do-      let !ac' = pushA i False (mc3SwapAcceptance a)+      let !ac' = pushReject i (mc3SwapAcceptance a)       return $ a {mc3SwapAcceptance = ac'}   where     g = mc3Generator a@@ -420,8 +421,9 @@   where     mhgs = mc3MHGChains a --- TODO: Splitmix. 'mc3Iterate' is actually not parallel, but concurrent because--- of the IO constraint. Use pure parallel code when we have a pure generator.+-- TODO (medium): Splitmix. 'mc3Iterate' is actually not parallel, but+-- concurrent because of the IO constraint. Use pure parallel code when we have+-- a pure generator. -- -- However, we have to take care of the mutable traces. mc3Iterate ::@@ -475,10 +477,10 @@     rNew = (brOld / blOld) ** xi     brNew = blNew * rNew -mc3AutoTune :: ToJSON a => Int -> MC3 a -> IO (MC3 a)-mc3AutoTune l a = do+mc3AutoTune :: ToJSON a => TuningType -> Int -> MC3 a -> IO (MC3 a)+mc3AutoTune b l a = do   -- 1. Auto tune all chains.-  mhgs' <- V.mapM (aAutoTune l) $ mc3MHGChains a+  mhgs' <- V.mapM (aAutoTune b l) $ mc3MHGChains a   -- 2. Auto tune temperatures.   let optimalRate = getOptimalRate PDimensionUnknown       mCurrentRates = acceptanceRates $ mc3SwapAcceptance a@@ -519,6 +521,13 @@     ac' = resetA $ mc3SwapAcceptance a     --     a' = a {mc3MHGChains = mhgs', mc3SwapAcceptance = ac'}++mc3CleanAfterBurnIn :: ToJSON a => TraceLength -> MC3 a -> IO (MC3 a)+mc3CleanAfterBurnIn tl a = do+  cs' <- V.mapM (aCleanAfterBurnIn tl) cs+  pure $ a {mc3MHGChains = cs'}+  where+    cs = mc3MHGChains a  -- Information in cycle summary: --

@@ -4,7 +4,7 @@ -- | -- Module      :  Mcmc.Algorithm.MHG -- Description :  Metropolis-Hastings-Green algorithm--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -34,6 +34,7 @@ import Control.Parallel.Strategies import Data.Aeson import qualified Data.ByteString.Lazy.Char8 as BL+import Data.Maybe import Data.Time import qualified Data.Vector as VB import Mcmc.Acceptance@@ -50,7 +51,7 @@ import Mcmc.Proposal import Mcmc.Settings import Numeric.Log-import System.Random.MWC+import System.Random.Stateful import Text.Printf import Prelude hiding (cycle) @@ -64,6 +65,7 @@   aIterate = mhgIterate   aAutoTune = mhgAutoTune   aResetAcceptance = mhgResetAcceptance+  aCleanAfterBurnIn = mhgCleanAfterBurnIn   aSummarizeCycle = mhgSummarizeCycle   aOpenMonitors = mhgOpenMonitors   aExecuteMonitors = mhgExecuteMonitors@@ -71,6 +73,23 @@   aCloseMonitors = mhgCloseMonitors   aSave = mhgSave +getTraceLength ::+  Maybe BurnInSettings ->+  TraceLength ->+  Monitor a ->+  Cycle a ->+  Int+getTraceLength burnIn tl mn cc = maximum $ minimumTraceLength : bi : batchMonitorSizes+  where+    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn+    minimumTraceLength = case tl of+      TraceAuto -> 1+      TraceMinimum n -> n+    bi = case (ccRequireTrace cc, burnIn) of+      (True, Just (BurnInWithAutoTuning _ n)) -> n+      (True, Just (BurnInWithCustomAutoTuning ns ms)) -> max (maximum $ 0 : ns) (maximum $ 0 : ms)+      _ -> 0+ -- | Initialize an MHG algorithm. -- -- NOTE: Computation in the 'IO' Monad is necessary because the trace is@@ -82,25 +101,17 @@   Cycle a ->   Monitor a ->   InitialState a ->-  GenIO ->+  IOGenM StdGen ->   IO (MHG a) mhg s pr lh cc mn i0 g = do   -- The trace is a mutable vector and the mutable state needs to be handled by   -- a monad.-  tr <- replicateT traceLength l0+  tr <- replicateT tl l0   return $ MHG $ Chain Nothing l0 0 tr ac g 0 pr lh cc mn   where     l0 = Link i0 (pr i0) (lh i0)     ac = emptyA $ ccProposals cc-    batchMonitorSizes = map getMonitorBatchSize $ mBatches mn-    minimumTraceLength = case sTraceLength s of-      TraceAuto -> 1-      TraceMinimum n -> n-    bi = case sBurnIn s of-      BurnInWithAutoTuning _ n -> n-      BurnInWithCustomAutoTuning ns ms -> max (maximum $ 0 : ns) (maximum $ 0 : ms)-      _ -> 0-    traceLength = maximum $ minimumTraceLength : bi : batchMonitorSizes+    tl = getTraceLength (Just $ sBurnIn s) (sTraceLength s) mn cc  mhgFn :: AnalysisName -> FilePath mhgFn (AnalysisName nm) = nm ++ ".mcmc.mhg"@@ -162,11 +173,14 @@  -- The MHG ratio. This implementation has the following properties: --+-- - The kernel ratio and the Jacobian are checked carefully and should be+-- - strictly positive, finite numbers.+-- -- - The ratio is 'Infinity' if fX is zero. In this case, the proposal is always --   accepted. ----- - The ratio 'NaN' if (fY or q or j) and fX are zero. In this case, the---   proposal is always rejected.+-- - The ratio is 'NaN' if fY and fX are zero. In this case, the proposal is+--   always rejected. -- -- This means that a chain in a state with posterior probability zero (fX=0) can -- only move if a state with non-zero posterior probability is proposed.@@ -182,49 +196,71 @@ -- Since I trust the author of Chapter 1 (Charles Geyer) I choose to follow -- option (b). However, Option (a) is more user-friendly. ----- [1] Handbook of markov chain monte carlo (2011), CRC press.+-- [1] Handbook of Markov chain Monte Carlo (2011), CRC press. ----- [2] Dellaportas, P., & Roberts, G. O., An introduction to mcmc, Lecture Notes+-- [2] Dellaportas, P., & Roberts, G. O., An introduction to MCMC, Lecture Notes -- in Statistics, (), 1–41 (2003). -- http://dx.doi.org/10.1007/978-0-387-21811-3_1. mhgRatio :: Posterior -> Posterior -> KernelRatio -> Jacobian -> MHGRatio--- q = qYX / qXY * jXY; see 'ProposalSimple'.--- j = Jacobian.-mhgRatio fX fY q j = fY / fX * q * j+mhgRatio fX fY q j+  | q == 0.0 = error "mhgRatio: Kernel ratio is negative infinity. Use 'ForceReject'."+  | q == 1.0 / 0.0 = error "mhgRatio: Kernel ratio is infinity. Use 'ForceAccept'."+  | q == 0.0 / 0.0 = error "mhgRatio: Kernel ratio is NaN."+  | j == 0.0 = error "mhgRatio: Jacobian is negative infinity. Use 'ForceReject'."+  | j == 1.0 / 0.0 = error "mhgRatio: Jacobian is infinity. Use 'ForceAccept'."+  | j == 0.0 / 0.0 = error "mhgRatio: Jacobian is NaN."+  | otherwise = fY / fX * q * j {-# INLINE mhgRatio #-}  -- | Accept or reject a proposal with given MHG ratio?-mhgAccept :: MHGRatio -> GenIO -> IO Bool+mhgAccept :: MHGRatio -> IOGenM StdGen -> IO Bool mhgAccept r g   | ln r >= 0.0 = return True   | otherwise = do-      b <- uniform g+      b <- uniformRM (0, 1) g       return $ b < exp (ln r)  mhgPropose :: MHG a -> Proposal a -> IO (MHG a) mhgPropose (MHG c) p = do   -- 1. Sample new state.-  (!y, !q, !j) <- liftIO $ s x g-  -- 2. Calculate Metropolis-Hastings-Green ratio.+  !(pres, mcs) <- liftIO $ s x g+  -- 2. Define new prior and likelihood calculation functions. Avoid actual+  -- calculation of the values.   --   -- Most often, parallelization is not helpful, because the prior and   -- likelihood functions are too fast; see   -- https://stackoverflow.com/a/46603680/3536806.-  let (pY, lY) = (pF y, lF y) `using` parTuple2 rdeepseq rdeepseq-  -- let !pY = pF y-  --     !lY = lF y-  let !r = mhgRatio (pX * lX) (pY * lY) q j+  let calcPrLh y = (pF y, lF y) `using` parTuple2 rdeepseq rdeepseq+      accept y pr lh =+        let !ac' = case mcs of+              Nothing -> pushAccept p ac+              Just cs -> pushAcceptanceCounts p cs ac+         in pure $ MHG $ c {link = Link y pr lh, acceptance = ac'}+      reject =+        let !ac' = case mcs of+              Nothing -> pushReject p ac+              Just cs -> pushAcceptanceCounts p cs ac+         in pure $ MHG $ c {acceptance = ac'}   -- 3. Accept or reject.-  accept <- mhgAccept r g-  if accept-    then do-      let !ac' = pushA p True ac-      return $ MHG $ c {link = Link y pY lY, acceptance = ac'}-    else do-      let !ac' = pushA p False ac-      return $ MHG $ c {acceptance = pushA p False ac'}+  --+  -- 3a. When rejection is inevitable, avoid calculation of the prior, the+  -- likelihood and the MHG ratio.+  case pres of+    ForceReject -> reject+    ForceAccept y -> let (pY, lY) = calcPrLh y in accept y pY lY+    (Propose y q j) ->+      if q <= 0.0 || j <= 0.0+        then reject+        else do+          -- 3b. Calculate Metropolis-Hastings-Green ratio.+          let (pY, lY) = calcPrLh y+              !r = mhgRatio (pX * lX) (pY * lY) q j+          isAccept <- mhgAccept r g+          if isAccept+            then accept y pY lY+            else reject   where-    s = prSimple p+    s = prFunction p     (Link x pX lX) = link c     pF = priorFunction c     lF = likelihoodFunction c@@ -245,12 +281,13 @@ -- At the moment this just checks whether the prior, likelihood, or posterior -- are NaN or infinite. mhgIsInValidState :: MHG a -> Bool-mhgIsInValidState a = check p || check l || check (p * l)+mhgIsInValidState a = checkSoft p || check l || check (p * l)   where     x = link $ fromMHG a     p = prior x     l = likelihood x     check v = let v' = ln v in isNaN v' || isInfinite v' || v' == 0+    checkSoft v = let v' = ln v in isNaN v' || isInfinite v'  -- Ignore the number of capabilities. I have tried a lot of stuff, but the MHG -- algorithm is just inherently sequential. Parallelization can be achieved by@@ -266,18 +303,34 @@     cc = cycle c     g = generator c -mhgAutoTune :: Int -> MHG a -> IO (MHG a)-mhgAutoTune n (MHG c) = do-  tr <- VB.map state <$> takeT n (trace c)-  return $ MHG $ c {cycle = autoTuneCycle ac tr cc}+mhgAutoTune :: TuningType -> Int -> MHG a -> IO (MHG a)+mhgAutoTune b n (MHG c) = do+  mxs <-+    if ccRequireTrace cc+      then Just . VB.map state <$> takeT n tr+      else pure Nothing+  return $ MHG $ c {cycle = autoTuneCycle b ac mxs cc}   where     ac = acceptance c     cc = cycle c+    tr = trace c  mhgResetAcceptance :: MHG a -> MHG a mhgResetAcceptance (MHG c) = MHG $ c {acceptance = resetA ac}   where     ac = acceptance c++mhgCleanAfterBurnIn :: TraceLength -> MHG a -> IO (MHG a)+mhgCleanAfterBurnIn tl (MHG c) = do+  xs <- takeT l tr+  tr' <- fromVectorT xs+  let c' = c {trace = tr'}+  pure $ MHG c'+  where+    mn = monitor c+    cc = cycle c+    tr = trace c+    l = getTraceLength Nothing tl mn cc  mhgSummarizeCycle :: IterationMode -> MHG a -> BL.ByteString mhgSummarizeCycle m (MHG c) = summarizeCycle m ac cc

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Chain.Chain -- Description :  Simple representation of a Markov chain--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -59,7 +59,7 @@ import Mcmc.Monitor import Mcmc.Prior import Mcmc.Proposal-import System.Random.MWC hiding (save)+import System.Random.Stateful import Prelude hiding (cycle)  -- | Type synonym to indicate the initial state.@@ -98,7 +98,7 @@     -- order.     acceptance :: Acceptance (Proposal a),     -- | The random number generator.-    generator :: GenIO,+    generator :: IOGenM StdGen,     --     -- Variables and functions; not saved. 

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Chain.Link -- Description :  The state combined with auxiliary variables--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -14,7 +14,7 @@ -- | -- Module      :  Mcmc.Chain.Save -- Description :  Save and load a Markov chain--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later module Mcmc.Chain.Save   ( SavedChain (..),@@ -31,7 +31,6 @@ import Data.Maybe import qualified Data.Stack.Circular as C import qualified Data.Vector as VB-import qualified Data.Vector.Unboxed as VU import Data.Word import Mcmc.Acceptance import Mcmc.Chain.Chain@@ -54,10 +53,10 @@     savedIteration :: Int,     savedTrace :: C.Stack VB.Vector (Link a),     savedAcceptance :: Acceptance Int,-    savedSeed :: VU.Vector Word32,+    savedSeed :: (Word64, Word64),     savedTuningParameters :: [Maybe (TuningParameter, AuxiliaryTuningParameters)]   }-  deriving (Eq, Read, Show)+  deriving (Eq, Show)  $(deriveJSON defaultOptions ''SavedChain) 

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Chain.Trace -- Description :  History of a Markov chain--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -12,6 +12,7 @@ module Mcmc.Chain.Trace   ( Trace,     replicateT,+    fromVectorT,     lengthT,     pushT,     headT,@@ -37,6 +38,12 @@ -- Call 'error' if the maximum size is zero or negative. replicateT :: Int -> Link a -> IO (Trace a) replicateT n l = Trace <$> C.replicate n l++-- | Create a trace from a vector. The length is determined by the vector.+--+-- Call 'error' if the vector is empty.+fromVectorT :: VB.Vector (Link a) -> IO (Trace a)+fromVectorT xs = Trace <$> C.fromVector xs  -- | Get the length of the trace. lengthT :: Trace a -> Int

@@ -4,7 +4,7 @@ -- | -- Module      :  Mcmc.Cycle -- Description :  A cycle is a list of proposals--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -15,7 +15,7 @@ module Mcmc.Cycle   ( -- * Cycles     Order (..),-    Cycle (ccProposals),+    Cycle (ccProposals, ccRequireTrace),     cycleFromList,     setOrder,     IterationMode (..),@@ -31,14 +31,13 @@ import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import Data.Default-import Data.Either import Data.List import qualified Data.Map.Strict as M import qualified Data.Vector as VB import Mcmc.Acceptance import Mcmc.Internal.Shuffle import Mcmc.Proposal-import System.Random.MWC+import System.Random.Stateful  -- | 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@@ -94,7 +93,9 @@ -- that they can be uniquely identified. data Cycle a = Cycle   { ccProposals :: [Proposal a],-    ccOrder :: Order+    ccOrder :: Order,+    -- | Does the cycle require the trace when auto tuning? See 'tRequireTrace'.+    ccRequireTrace :: Bool   }  -- | Create a 'Cycle' from a list of 'Proposal's.@@ -103,7 +104,7 @@   error "cycleFromList: Received an empty list but cannot create an empty Cycle." cycleFromList xs =   if length uniqueXs == length xs-    then Cycle xs def+    then Cycle xs def (any needsTrace xs)     else error $ "\n" ++ msg ++ "cycleFromList: Proposals are not unique."   where     uniqueXs = nub xs@@ -112,6 +113,9 @@     removedDescriptions = map (show . prDescription) removedXs     removedMsgs = zipWith (\n d -> n ++ " " ++ d) removedNames removedDescriptions     msg = unlines removedMsgs+    needsTrace p = case prTuner p of+      Nothing -> False+      Just t -> tRequireTrace t  -- | Set the order of 'Proposal's in a 'Cycle'. setOrder :: Order -> Cycle a -> Cycle a@@ -122,14 +126,17 @@   deriving (Eq)  -- | Replicate 'Proposal's according to their weights and possibly shuffle them.-prepareProposals :: IterationMode -> Cycle a -> GenIO -> IO [Proposal a]-prepareProposals m (Cycle xs o) g = case o of-  RandomO -> shuffle ps g-  SequentialO -> return ps-  RandomReversibleO -> do-    psR <- shuffle ps g-    return $ psR ++ reverse psR-  SequentialReversibleO -> return $ ps ++ reverse ps+prepareProposals :: StatefulGen g m => IterationMode -> Cycle a -> g -> m [Proposal a]+prepareProposals m (Cycle xs o _) g =+  if null ps+    then error "prepareProposals: No proposals found."+    else case o of+      RandomO -> shuffle ps g+      SequentialO -> return ps+      RandomReversibleO -> do+        psR <- shuffle ps g+        return $ psR ++ reverse psR+      SequentialReversibleO -> return $ ps ++ reverse ps   where     !ps =       concat@@ -143,7 +150,7 @@  -- The number of proposals depends on the order. getNProposalsPerCycle :: IterationMode -> Cycle a -> Int-getNProposalsPerCycle m (Cycle xs o) = case o of+getNProposalsPerCycle m (Cycle xs o _) = case o of   RandomO -> once   SequentialO -> once   RandomReversibleO -> 2 * once@@ -155,28 +162,34 @@     once = sum $ map (fromPWeight . prWeight) xs'  -- See 'tuneWithTuningParameters' and 'Tuner'.-tuneWithChainParameters :: AcceptanceRate -> VB.Vector a -> Proposal a -> Either String (Proposal a)-tuneWithChainParameters ar xs p = case prTuner p of-  Nothing -> Left "tuneWithChainParameters: Proposal is not tunable."-  Just (Tuner t fT ts fTs _) ->-    let t' = fT ar t-        ts' = fTs xs ts-     in tuneWithTuningParameters t' ts' p+tuneWithChainParameters :: TuningType -> AcceptanceRate -> Maybe (VB.Vector a) -> Proposal a -> Either String (Proposal a)+tuneWithChainParameters b ar mxs p = case prTuner p of+  Nothing -> Right p+  Just (Tuner t ts rt fT _) -> case (rt, mxs) of+    (True, Nothing) -> error "tuneWithChainParameters: trace required"+    _ ->+      let (t', ts') = fT b d ar mxs (t, ts)+       in tuneWithTuningParameters t' ts' p+      where+        d = prDimension p  -- | Calculate acceptance rates and auto tunes the 'Proposal's in the 'Cycle'. -- -- Do not change 'Proposal's that are not tuneable.-autoTuneCycle :: Acceptance (Proposal a) -> VB.Vector a -> Cycle a -> Cycle a-autoTuneCycle a xs c =-  if sort (M.keys ar) == sort ps-    then c {ccProposals = map tuneF ps}-    else error "autoTuneCycle: Proposals in map and cycle do not match."-  where-    ar = acceptanceRates a-    ps = ccProposals c-    tuneF p = case ar M.!? p of-      Just (Just x) -> fromRight p (tuneWithChainParameters x xs p)-      _ -> p+autoTuneCycle :: TuningType -> Acceptance (Proposal a) -> Maybe (VB.Vector a) -> Cycle a -> Cycle a+autoTuneCycle b a mxs c = case (ccRequireTrace c, mxs) of+  (False, Just _) -> error "autoTuneCycle: trace not required"+  (True, Nothing) -> error "autoTuneCycle: trace required"+  _ ->+    if sort (M.keys ar) == sort ps+      then c {ccProposals = map tuneF ps}+      else error "autoTuneCycle: Proposals in map and cycle do not match."+    where+      ar = acceptanceRates a+      ps = ccProposals c+      tuneF p = case ar M.!? p of+        Just (Just x) -> either error id $ tuneWithChainParameters b x mxs p+        _ -> p  -- | Horizontal line of proposal summaries. proposalHLine :: BL.ByteString

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Environment -- Description :  Runtime environment--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Internal.ByteString -- Description :  ByteString tools--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,205 +0,0 @@--- |--- Module      :  Mcmc.Internal.Gamma--- Description :  Generalized gamma function for automatic differentiation--- Copyright   :  (c) 2021 Dominik Schrempf--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  experimental--- Portability :  portable------ Creation date: Tue Jul 13 12:53:09 2021.------ The code is taken from "Numeric.SpecFunctions".-module Mcmc.Internal.Gamma-  ( logGammaG,-  )-where--import Data.Typeable-import qualified Data.Vector as VB-import Numeric.Polynomial-import Numeric.SpecFunctions-import Unsafe.Coerce--mSqrtEpsG :: RealFloat a => a-mSqrtEpsG = 1.4901161193847656e-8--mEulerMascheroniG :: RealFloat a => a-mEulerMascheroniG = 0.5772156649015328606065121---- | Generalized version of the log gamma distribution. See--- 'Numeric.SpecFunctions.logGamma'.-logGammaG :: (Typeable a, RealFloat a) => a -> a-logGammaG z-  | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z-  | otherwise = logGammaNonDouble z-{-# SPECIALIZE logGammaG :: Double -> Double #-}---- | See 'Numeric.SpecFunctions.logGamma'.-logGammaNonDouble :: RealFloat a => a -> a-logGammaNonDouble z-  | z <= 0 = 1 / 0-  | z < mSqrtEpsG = log (1 / z - mEulerMascheroniG)-  | z < 0.5 = lgamma1_15G z (z - 1) - log z-  | z < 1 = lgamma15_2G z (z - 1) - log z-  | z <= 1.5 = lgamma1_15G (z - 1) (z - 2)-  | z < 2 = lgamma15_2G (z - 1) (z - 2)-  | z < 15 = lgammaSmallG z-  | otherwise = lanczosApproxG z--lgamma1_15G :: RealFloat a => a -> a -> a-lgamma1_15G zm1 zm2 =-  r * y + r-    * ( evaluatePolynomial zm1 tableLogGamma_1_15PG-          / evaluatePolynomial zm1 tableLogGamma_1_15QG-      )-  where-    r = zm1 * zm2-    y = 0.52815341949462890625--tableLogGamma_1_15PG :: RealFloat a => VB.Vector a-tableLogGamma_1_15PG =-  VB.fromList-    [ 0.490622454069039543534e-1,-      -0.969117530159521214579e-1,-      -0.414983358359495381969e0,-      -0.406567124211938417342e0,-      -0.158413586390692192217e0,-      -0.240149820648571559892e-1,-      -0.100346687696279557415e-2-    ]-{-# NOINLINE tableLogGamma_1_15PG #-}--tableLogGamma_1_15QG :: RealFloat a => VB.Vector a-tableLogGamma_1_15QG =-  VB.fromList-    [ 1,-      0.302349829846463038743e1,-      0.348739585360723852576e1,-      0.191415588274426679201e1,-      0.507137738614363510846e0,-      0.577039722690451849648e-1,-      0.195768102601107189171e-2-    ]-{-# NOINLINE tableLogGamma_1_15QG #-}--lgamma15_2G :: RealFloat a => a -> a -> a-lgamma15_2G zm1 zm2 =-  r * y + r-    * ( evaluatePolynomial (-zm2) tableLogGamma_15_2PG-          / evaluatePolynomial (-zm2) tableLogGamma_15_2QG-      )-  where-    r = zm1 * zm2-    y = 0.452017307281494140625--tableLogGamma_15_2PG :: RealFloat a => VB.Vector a-tableLogGamma_15_2PG =-  VB.fromList-    [ -0.292329721830270012337e-1,-      0.144216267757192309184e0,-      -0.142440390738631274135e0,-      0.542809694055053558157e-1,-      -0.850535976868336437746e-2,-      0.431171342679297331241e-3-    ]-{-# NOINLINE tableLogGamma_15_2PG #-}--tableLogGamma_15_2QG :: RealFloat a => VB.Vector a-tableLogGamma_15_2QG =-  VB.fromList-    [ 1,-      -0.150169356054485044494e1,-      0.846973248876495016101e0,-      -0.220095151814995745555e0,-      0.25582797155975869989e-1,-      -0.100666795539143372762e-2,-      -0.827193521891290553639e-6-    ]-{-# NOINLINE tableLogGamma_15_2QG #-}--lgammaSmallG :: RealFloat a => a -> a-lgammaSmallG = go 0-  where-    go acc z-      | z < 3 = acc + lgamma2_3G z-      | otherwise = go (acc + log zm1) zm1-      where-        zm1 = z - 1--lgamma2_3G :: RealFloat a => a -> a-lgamma2_3G z =-  r * y + r-    * ( evaluatePolynomial zm2 tableLogGamma_2_3PG-          / evaluatePolynomial zm2 tableLogGamma_2_3QG-      )-  where-    r = zm2 * (z + 1)-    zm2 = z - 2-    y = 0.158963680267333984375e0--tableLogGamma_2_3PG :: RealFloat a => VB.Vector a-tableLogGamma_2_3PG =-  VB.fromList-    [ -0.180355685678449379109e-1,-      0.25126649619989678683e-1,-      0.494103151567532234274e-1,-      0.172491608709613993966e-1,-      -0.259453563205438108893e-3,-      -0.541009869215204396339e-3,-      -0.324588649825948492091e-4-    ]-{-# NOINLINE tableLogGamma_2_3PG #-}--tableLogGamma_2_3QG :: RealFloat a => VB.Vector a-tableLogGamma_2_3QG =-  VB.fromList-    [ 1,-      0.196202987197795200688e1,-      0.148019669424231326694e1,-      0.541391432071720958364e0,-      0.988504251128010129477e-1,-      0.82130967464889339326e-2,-      0.224936291922115757597e-3,-      -0.223352763208617092964e-6-    ]-{-# NOINLINE tableLogGamma_2_3QG #-}--lanczosApproxG :: RealFloat a => a -> a-lanczosApproxG z =-  (log (z + g - 0.5) - 1) * (z - 0.5)-    + log (evalRatioG tableLanczosG z)-  where-    g = 6.024680040776729583740234375--tableLanczosG :: RealFloat a => VB.Vector (a, a)-tableLanczosG =-  VB.fromList-    [ (56906521.91347156388090791033559122686859, 0),-      (103794043.1163445451906271053616070238554, 39916800),-      (86363131.28813859145546927288977868422342, 120543840),-      (43338889.32467613834773723740590533316085, 150917976),-      (14605578.08768506808414169982791359218571, 105258076),-      (3481712.15498064590882071018964774556468, 45995730),-      (601859.6171681098786670226533699352302507, 13339535),-      (75999.29304014542649875303443598909137092, 2637558),-      (6955.999602515376140356310115515198987526, 357423),-      (449.9445569063168119446858607650988409623, 32670),-      (19.51992788247617482847860966235652136208, 1925),-      (0.5098416655656676188125178644804694509993, 66),-      (0.006061842346248906525783753964555936883222, 1)-    ]-{-# NOINLINE tableLanczosG #-}--data LG a = LG !a !a--evalRatioG :: RealFloat a => VB.Vector (a, a) -> a -> a-evalRatioG coef x-  | x > 1 = fini $ VB.foldl' stepL (LG 0 0) coef-  | otherwise = fini $ VB.foldr' stepR (LG 0 0) coef-  where-    fini (LG num den) = num / den-    stepR (a, b) (LG num den) = LG (num * x + a) (den * x + b)-    stepL (LG num den) (a, b) = LG (num * rx + a) (den * rx + b)-    rx = recip x

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Internal.Random -- Description :  Tools for random calculations--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -12,37 +12,23 @@ -- -- Creation date: Wed Nov 25 07:14:52 2020. module Mcmc.Internal.Random-  ( splitGen,-    saveGen,+  ( saveGen,     loadGen,   ) where -import Control.Monad-import Control.Monad.Primitive-import qualified Data.Vector.Unboxed as V+import Data.IORef import Data.Word-import System.Random.MWC---- | Split a generator.------ Splitting an MWC pseudo number generator is not good practice. However, I--- have to go with this solution for now, and wait for proper support of--- spittable pseudo random number generators such as @splitmix@.-splitGen :: PrimMonad m => Int -> Gen (PrimState m) -> m [Gen (PrimState m)]-splitGen n gen-  | n <= 0 = return []-  | otherwise = do-      seeds :: [V.Vector Word32] <- replicateM n $ uniformVector gen 256-      mapM initialize seeds---- TODO: Splitmix. Remove or amend these functions as soon as split mix is used--- and is available with the statistics package.+import System.Random.Internal+import System.Random.SplitMix+import System.Random.Stateful  -- | Save a generator to a seed.-saveGen :: GenIO -> IO (V.Vector Word32)-saveGen = fmap fromSeed . save+saveGen :: IOGenM StdGen -> IO (Word64, Word64)+saveGen (IOGenM r) = do+  (StdGen g) <- readIORef r+  pure $ unseedSMGen g  -- | Load a generator from a seed.-loadGen :: V.Vector Word32 -> IO GenIO-loadGen = restore . toSeed+loadGen :: (Word64, Word64) -> IO (IOGenM StdGen)+loadGen s = newIOGenM $ StdGen $ seedSMGen' s

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Internal.Shuffle -- Description :  Shuffle a list--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -21,23 +21,26 @@ import Control.Monad.ST import qualified Data.Vector as V import qualified Data.Vector.Mutable as M-import System.Random.MWC+import System.Random.Stateful  -- Fisher-Yates shuffle. See also -- 'System.Random.MWC.Distributions.uniformPermutation' which is a little -- cleaner, in my opinion. However, I would like to move away from MWC so I -- leave the custom implementation for now.+--+-- NOTE: Moving away from MWC is not anymore true because I can use SplitMix+-- with it.  -- | Shuffle a vector.-shuffle :: [a] -> GenIO -> IO [a]+shuffle :: StatefulGen g m => [a] -> g -> m [a] shuffle xs = grabble xs (length xs)  -- | @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 -> GenIO -> IO [a]+grabble :: StatefulGen g m => [a] -> Int -> g -> m [a] grabble xs m gen = do   swaps <- forM [0 .. min (l - 1) m] $ \i -> do-    j <- uniformR (i, l) gen+    j <- uniformRM (i, l) gen     return (i, j)   return $ (V.toList . V.take m . swapElems (V.fromList xs)) swaps   where

@@ -0,0 +1,410 @@+{-# LANGUAGE ScopedTypeVariables #-}++-- |+-- Module      :  Mcmc.Internal.Gamma+-- Description :  Generalized gamma function for automatic differentiation+-- Copyright   :  2021 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Tue Jul 13 12:53:09 2021.+--+-- The code is taken from "Numeric.SpecFunctions".+module Mcmc.Internal.SpecFunctions+  ( logGammaG,+    logFactorialG,+  )+where++import Data.Typeable+import qualified Data.Vector as VB+import Numeric.Polynomial+import Numeric.SpecFunctions+import Unsafe.Coerce++mSqrtEpsG :: RealFloat a => a+mSqrtEpsG = 1.4901161193847656e-8++mEulerMascheroniG :: RealFloat a => a+mEulerMascheroniG = 0.5772156649015328606065121++-- | Generalized version of the log gamma distribution. See+-- 'Numeric.SpecFunctions.logGamma'.+logGammaG :: (Typeable a, RealFloat a) => a -> a+logGammaG z+  | typeOf z == typeOf (0 :: Double) = unsafeCoerce logGamma z+  | otherwise = logGammaNonDouble z+{-# SPECIALIZE logGammaG :: Double -> Double #-}++-- | See 'Numeric.SpecFunctions.logGamma'.+logGammaNonDouble :: RealFloat a => a -> a+logGammaNonDouble z+  | z <= 0 = 1 / 0+  | z < mSqrtEpsG = log (1 / z - mEulerMascheroniG)+  | z < 0.5 = lgamma1_15G z (z - 1) - log z+  | z < 1 = lgamma15_2G z (z - 1) - log z+  | z <= 1.5 = lgamma1_15G (z - 1) (z - 2)+  | z < 2 = lgamma15_2G (z - 1) (z - 2)+  | z < 15 = lgammaSmallG z+  | otherwise = lanczosApproxG z++lgamma1_15G :: RealFloat a => a -> a -> a+lgamma1_15G zm1 zm2 =+  r * y + r+    * ( evaluatePolynomial zm1 tableLogGamma_1_15PG+          / evaluatePolynomial zm1 tableLogGamma_1_15QG+      )+  where+    r = zm1 * zm2+    y = 0.52815341949462890625++tableLogGamma_1_15PG :: RealFloat a => VB.Vector a+tableLogGamma_1_15PG =+  VB.fromList+    [ 0.490622454069039543534e-1,+      -0.969117530159521214579e-1,+      -0.414983358359495381969e0,+      -0.406567124211938417342e0,+      -0.158413586390692192217e0,+      -0.240149820648571559892e-1,+      -0.100346687696279557415e-2+    ]+{-# NOINLINE tableLogGamma_1_15PG #-}++tableLogGamma_1_15QG :: RealFloat a => VB.Vector a+tableLogGamma_1_15QG =+  VB.fromList+    [ 1,+      0.302349829846463038743e1,+      0.348739585360723852576e1,+      0.191415588274426679201e1,+      0.507137738614363510846e0,+      0.577039722690451849648e-1,+      0.195768102601107189171e-2+    ]+{-# NOINLINE tableLogGamma_1_15QG #-}++lgamma15_2G :: RealFloat a => a -> a -> a+lgamma15_2G zm1 zm2 =+  r * y + r+    * ( evaluatePolynomial (-zm2) tableLogGamma_15_2PG+          / evaluatePolynomial (-zm2) tableLogGamma_15_2QG+      )+  where+    r = zm1 * zm2+    y = 0.452017307281494140625++tableLogGamma_15_2PG :: RealFloat a => VB.Vector a+tableLogGamma_15_2PG =+  VB.fromList+    [ -0.292329721830270012337e-1,+      0.144216267757192309184e0,+      -0.142440390738631274135e0,+      0.542809694055053558157e-1,+      -0.850535976868336437746e-2,+      0.431171342679297331241e-3+    ]+{-# NOINLINE tableLogGamma_15_2PG #-}++tableLogGamma_15_2QG :: RealFloat a => VB.Vector a+tableLogGamma_15_2QG =+  VB.fromList+    [ 1,+      -0.150169356054485044494e1,+      0.846973248876495016101e0,+      -0.220095151814995745555e0,+      0.25582797155975869989e-1,+      -0.100666795539143372762e-2,+      -0.827193521891290553639e-6+    ]+{-# NOINLINE tableLogGamma_15_2QG #-}++lgammaSmallG :: RealFloat a => a -> a+lgammaSmallG = go 0+  where+    go acc z+      | z < 3 = acc + lgamma2_3G z+      | otherwise = go (acc + log zm1) zm1+      where+        zm1 = z - 1++lgamma2_3G :: RealFloat a => a -> a+lgamma2_3G z =+  r * y + r+    * ( evaluatePolynomial zm2 tableLogGamma_2_3PG+          / evaluatePolynomial zm2 tableLogGamma_2_3QG+      )+  where+    r = zm2 * (z + 1)+    zm2 = z - 2+    y = 0.158963680267333984375e0++tableLogGamma_2_3PG :: RealFloat a => VB.Vector a+tableLogGamma_2_3PG =+  VB.fromList+    [ -0.180355685678449379109e-1,+      0.25126649619989678683e-1,+      0.494103151567532234274e-1,+      0.172491608709613993966e-1,+      -0.259453563205438108893e-3,+      -0.541009869215204396339e-3,+      -0.324588649825948492091e-4+    ]+{-# NOINLINE tableLogGamma_2_3PG #-}++tableLogGamma_2_3QG :: RealFloat a => VB.Vector a+tableLogGamma_2_3QG =+  VB.fromList+    [ 1,+      0.196202987197795200688e1,+      0.148019669424231326694e1,+      0.541391432071720958364e0,+      0.988504251128010129477e-1,+      0.82130967464889339326e-2,+      0.224936291922115757597e-3,+      -0.223352763208617092964e-6+    ]+{-# NOINLINE tableLogGamma_2_3QG #-}++lanczosApproxG :: RealFloat a => a -> a+lanczosApproxG z =+  (log (z + g - 0.5) - 1) * (z - 0.5)+    + log (evalRatioG tableLanczosG z)+  where+    g = 6.024680040776729583740234375++tableLanczosG :: RealFloat a => VB.Vector (a, a)+tableLanczosG =+  VB.fromList+    [ (56906521.91347156388090791033559122686859, 0),+      (103794043.1163445451906271053616070238554, 39916800),+      (86363131.28813859145546927288977868422342, 120543840),+      (43338889.32467613834773723740590533316085, 150917976),+      (14605578.08768506808414169982791359218571, 105258076),+      (3481712.15498064590882071018964774556468, 45995730),+      (601859.6171681098786670226533699352302507, 13339535),+      (75999.29304014542649875303443598909137092, 2637558),+      (6955.999602515376140356310115515198987526, 357423),+      (449.9445569063168119446858607650988409623, 32670),+      (19.51992788247617482847860966235652136208, 1925),+      (0.5098416655656676188125178644804694509993, 66),+      (0.006061842346248906525783753964555936883222, 1)+    ]+{-# NOINLINE tableLanczosG #-}++data LG a = LG !a !a++evalRatioG :: RealFloat a => VB.Vector (a, a) -> a -> a+evalRatioG coef x+  | x > 1 = fini $ VB.foldl' stepL (LG 0 0) coef+  | otherwise = fini $ VB.foldr' stepR (LG 0 0) coef+  where+    fini (LG num den) = num / den+    stepR (a, b) (LG num den) = LG (num * x + a) (den * x + b)+    stepL (LG num den) (a, b) = LG (num * rx + a) (den * rx + b)+    rx = recip x++-- | Generalized version of the log factorial function. See+-- 'Numeric.SpecFunctions.logFactorial'.+logFactorialG :: forall a b. (Integral a, RealFloat b, Typeable b) => a -> b+logFactorialG n+  | typeOf (undefined :: b) == typeOf (0 :: Double) = unsafeCoerce $ logFactorial n+  | otherwise = logFactorialNonDouble n+{-# SPECIALIZE logFactorialG :: Int -> Double #-}++logFactorialNonDouble :: (Integral a, RealFloat b) => a -> b+logFactorialNonDouble n+  | n < 0 = error "logFactorialNonDouble: Negative input."+  | n <= 170 = log $ VB.unsafeIndex factorialTable (fromIntegral n)+  | n < 1500 = stirling + rx * ((1 / 12) - (1 / 360) * rx * rx)+  | otherwise = stirling + (1 / 12) * rx+  where+    stirling = (x - 0.5) * log x - x + mLnSqrt2Pi+    x = fromIntegral n + 1+    rx = recip x+{-# SPECIALIZE logFactorialNonDouble :: RealFloat a => Int -> a #-}++mLnSqrt2Pi :: RealFloat a => a+mLnSqrt2Pi = 0.9189385332046727417803297364056176398613974736377834128171+{-# INLINE mLnSqrt2Pi #-}++factorialTable :: RealFloat a => VB.Vector a+{-# NOINLINE factorialTable #-}+factorialTable =+  VB.fromListN+    171+    [ 1.0,+      1.0,+      2.0,+      6.0,+      24.0,+      120.0,+      720.0,+      5040.0,+      40320.0,+      362880.0,+      3628800.0,+      3.99168e7,+      4.790016e8,+      6.2270208e9,+      8.71782912e10,+      1.307674368e12,+      2.0922789888e13,+      3.55687428096e14,+      6.402373705728e15,+      1.21645100408832e17,+      2.43290200817664e18,+      5.109094217170944e19,+      1.1240007277776077e21,+      2.5852016738884974e22,+      6.204484017332394e23,+      1.5511210043330984e25,+      4.032914611266056e26,+      1.0888869450418352e28,+      3.0488834461171384e29,+      8.841761993739702e30,+      2.6525285981219103e32,+      8.222838654177922e33,+      2.631308369336935e35,+      8.683317618811886e36,+      2.9523279903960412e38,+      1.0333147966386144e40,+      3.719933267899012e41,+      1.3763753091226343e43,+      5.23022617466601e44,+      2.0397882081197442e46,+      8.159152832478977e47,+      3.3452526613163803e49,+      1.4050061177528798e51,+      6.041526306337383e52,+      2.6582715747884485e54,+      1.1962222086548019e56,+      5.5026221598120885e57,+      2.5862324151116818e59,+      1.2413915592536073e61,+      6.082818640342675e62,+      3.0414093201713376e64,+      1.5511187532873822e66,+      8.065817517094388e67,+      4.2748832840600255e69,+      2.308436973392414e71,+      1.2696403353658275e73,+      7.109985878048634e74,+      4.0526919504877214e76,+      2.3505613312828785e78,+      1.386831185456898e80,+      8.32098711274139e81,+      5.075802138772247e83,+      3.146997326038793e85,+      1.9826083154044399e87,+      1.2688693218588415e89,+      8.24765059208247e90,+      5.44344939077443e92,+      3.647111091818868e94,+      2.4800355424368305e96,+      1.711224524281413e98,+      1.197857166996989e100,+      8.504785885678623e101,+      6.1234458376886085e103,+      4.470115461512684e105,+      3.307885441519386e107,+      2.4809140811395396e109,+      1.88549470166605e111,+      1.4518309202828586e113,+      1.1324281178206297e115,+      8.946182130782974e116,+      7.15694570462638e118,+      5.797126020747368e120,+      4.753643337012841e122,+      3.9455239697206583e124,+      3.314240134565353e126,+      2.81710411438055e128,+      2.422709538367273e130,+      2.1077572983795275e132,+      1.8548264225739844e134,+      1.650795516090846e136,+      1.4857159644817613e138,+      1.352001527678403e140,+      1.2438414054641305e142,+      1.1567725070816416e144,+      1.087366156656743e146,+      1.0329978488239058e148,+      9.916779348709496e149,+      9.619275968248211e151,+      9.426890448883246e153,+      9.332621544394413e155,+      9.332621544394415e157,+      9.425947759838358e159,+      9.614466715035125e161,+      9.902900716486179e163,+      1.0299016745145626e166,+      1.0813967582402908e168,+      1.1462805637347082e170,+      1.2265202031961378e172,+      1.3246418194518288e174,+      1.4438595832024934e176,+      1.5882455415227428e178,+      1.7629525510902446e180,+      1.974506857221074e182,+      2.2311927486598134e184,+      2.543559733472187e186,+      2.9250936934930154e188,+      3.393108684451898e190,+      3.9699371608087206e192,+      4.68452584975429e194,+      5.574585761207606e196,+      6.689502913449126e198,+      8.094298525273443e200,+      9.875044200833601e202,+      1.214630436702533e205,+      1.5061417415111406e207,+      1.8826771768889257e209,+      2.372173242880047e211,+      3.0126600184576594e213,+      3.856204823625804e215,+      4.974504222477286e217,+      6.466855489220473e219,+      8.471580690878819e221,+      1.1182486511960041e224,+      1.4872707060906857e226,+      1.9929427461615188e228,+      2.6904727073180504e230,+      3.6590428819525483e232,+      5.012888748274991e234,+      6.917786472619488e236,+      9.615723196941088e238,+      1.3462012475717523e241,+      1.898143759076171e243,+      2.6953641378881624e245,+      3.8543707171800725e247,+      5.5502938327393044e249,+      8.047926057471992e251,+      1.1749972043909107e254,+      1.7272458904546386e256,+      2.5563239178728654e258,+      3.808922637630569e260,+      5.713383956445854e262,+      8.62720977423324e264,+      1.3113358856834524e267,+      2.0063439050956823e269,+      3.0897696138473508e271,+      4.789142901463393e273,+      7.471062926282894e275,+      1.1729568794264143e278,+      1.8532718694937346e280,+      2.946702272495038e282,+      4.714723635992061e284,+      7.590705053947218e286,+      1.2296942187394494e289,+      2.0044015765453023e291,+      3.287218585534296e293,+      5.423910666131589e295,+      9.003691705778436e297,+      1.5036165148649988e300,+      2.526075744973198e302,+      4.269068009004705e304,+      7.257415615307998e306+    ]

@@ -0,0 +1,32 @@+-- |+-- Module      :  Mcmc.Jacobian+-- Description :  Types and convenience functions for Jacobians+-- Copyright   :  2022 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Tue May 31 09:37:54 2022.+module Mcmc.Jacobian+  ( Jacobian,+    JacobianG,+    JacobianFunction,+    JacobianFunctionG,+  )+where++import Numeric.Log++-- | Absolute value of the determinant of the Jacobian matrix.+type Jacobian = Log Double++-- | Generalized Jacobian.+type JacobianG a = Log a++-- | Function calculating the 'Jacobian'.+type JacobianFunction a = JacobianFunctionG a Double++-- | Function calculating the 'Jacobian'.+type JacobianFunctionG a b = a -> JacobianG b

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Likelihood -- Description :  Types and convenience functions for computing likelihoods--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -4,7 +4,7 @@ -- | -- Module      :  Mcmc.Logger -- Description :  Minimal monad logger--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -4,7 +4,7 @@ -- | -- Module      :  Mcmc.MarginalLikelihood -- Description :  Calculate the marginal likelihood--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -38,7 +38,6 @@ import Mcmc.Chain.Trace import Mcmc.Cycle import Mcmc.Environment-import Mcmc.Internal.Random import Mcmc.Likelihood import Mcmc.Logger import Mcmc.Mcmc@@ -47,7 +46,7 @@ import Mcmc.Settings import Numeric.Log hiding (sum) import System.Directory-import System.Random.MWC+import System.Random.Stateful import Text.Printf import Text.Show.Pretty import Prelude hiding (cycle)@@ -224,7 +223,7 @@   Cycle a ->   Monitor a ->   a ->-  GenIO ->+  IOGenM StdGen ->   -- For each point a vector of likelihoods stored in log domain.   ML [VU.Vector Likelihood] mlRun k xs em vb prf lhf cc mn i0 g = do@@ -271,11 +270,12 @@   Cycle a ->   Monitor a ->   a ->-  GenIO ->+  IOGenM StdGen ->   ML MarginalLikelihood-tiWrapper s prf lhf cc mn i0 g = do+tiWrapper s prf lhf cc mn i0 g0 = do   logInfoB "Path integral (thermodynamic integration)."-  [g0, g1] <- splitGen 2 g+  r1 <- splitGenM g0+  g1 <- newIOGenM r1    -- Parallel execution of both path integrals.   r <- ask@@ -354,7 +354,7 @@   Cycle a ->   Monitor a ->   a ->-  GenIO ->+  IOGenM StdGen ->   ML MarginalLikelihood sssWrapper s prf lhf cc mn i0 g = do   logInfoB "Stepping stone sampling."@@ -378,9 +378,7 @@   Cycle a ->   Monitor a ->   InitialState a ->-  -- | A source of randomness. For reproducible runs, make sure to use-  -- generators with the same seed.-  GenIO ->+  IOGenM StdGen ->   IO MarginalLikelihood marginalLikelihood s prf lhf cc mn i0 g = do   -- Initialize.@@ -401,7 +399,7 @@           ThermodynamicIntegration -> tiWrapper s prf lhf cc mn i0 g           SteppingStoneSampling -> sssWrapper s prf lhf cc mn i0 g         logInfoS $ "Marginal log likelihood: " ++ show (ln val)-        -- TODO: Simulation variance.+        -- TODO (low): Simulation variance.         logInfoS "The simulation variance is not yet available."         return val     )

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Mcmc -- Description :  Framework for running Markov chain Monte Carlo samplers--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -29,6 +29,7 @@ import Mcmc.Cycle import Mcmc.Environment import Mcmc.Logger+import Mcmc.Proposal (TuningType (LastTuningStep, NormalTuningStep)) import Mcmc.Settings import System.Exit import System.IO@@ -55,25 +56,25 @@   return $ aResetAcceptance a  mcmcExceptionHandler :: Algorithm a => Environment Settings -> a -> AsyncException -> IO b-mcmcExceptionHandler e a UserInterrupt = do+mcmcExceptionHandler e a err = do   putStrLn ""-  putStrLn "USER INTERRUPT!"+  putStrLn "INTERRUPT!"   putStrLn "Try to terminate gracefully and save chain for continuation."-  putStrLn "Press CTRL-C again to terminate now."-  putStrLn "Close output files."+  putStrLn "Press CTRL-C (again) to terminate now."+  putStrLn "Closing output files."   _ <- aCloseMonitors a   closeEnvironment e-  putStrLn "Try to save settings."+  putStrLn "Saving settings."   let s = settings e   settingsSave s-  putStrLn "Try to save compressed MCMC analysis."+  putStrLn "Saving compressed MCMC analysis."   putStrLn "For long traces, or complex objects, this may take a while."   let nm = sAnalysisName s   aSave nm a   putStrLn "Markov chain saved. Analysis can be continued."-  putStrLn "Terminate gracefully."-  exitWith $ ExitFailure 1-mcmcExceptionHandler _ _ e = throw e+  putStrLn "Graceful termination successful."+  putStrLn "Rethrowing error."+  throw err  -- XXX: Exception handling. Is it enough to mask execution of monitors and catch -- UserInterrupt during iterations?@@ -110,10 +111,13 @@   mcmcExecuteMonitors a   when (aIsInValidState a) (logWarnB "The initial state is invalid!")   a' <- mcmcBurnIn a+  logInfoS $ "Clean chain after burn in."+  let tl = sTraceLength s+  a'' <- liftIO $ aCleanAfterBurnIn tl a'   let i = fromIterations $ sIterations s   logInfoS $ "Run chain for " ++ show i ++ " iterations."-  logInfoB $ aStdMonitorHeader a'-  mcmcIterate AllProposals i a'+  logInfoB $ aStdMonitorHeader a''+  mcmcIterate AllProposals i a''  mcmcContinueRun :: Algorithm a => a -> MCMC a mcmcContinueRun a = do@@ -183,23 +187,26 @@       return a''  -- Auto tune the proposals.-mcmcAutotune :: Algorithm a => Int -> a -> MCMC a-mcmcAutotune n a = do+mcmcAutotune :: Algorithm a => TuningType -> Int -> a -> MCMC a+mcmcAutotune NormalTuningStep n a = do   logDebugB "Auto tune."-  liftIO $ aAutoTune n a+  liftIO $ aAutoTune NormalTuningStep n a+mcmcAutotune LastTuningStep n a = do+  logDebugB "Last auto tune."+  liftIO $ aAutoTune LastTuningStep n a  mcmcBurnInWithAutoTuning :: Algorithm a => IterationMode -> [Int] -> a -> MCMC a mcmcBurnInWithAutoTuning _ [] _ = error "mcmcBurnInWithAutoTuning: Empty list." mcmcBurnInWithAutoTuning m [x] a = do   -- Last round.   a' <- mcmcIterate m x a-  a'' <- mcmcAutotune x a'+  a'' <- mcmcAutotune LastTuningStep x a'   logInfoB $ aSummarizeCycle m a''   logInfoS $ "Acceptance rates calculated over the last " <> show x <> " iterations."   mcmcResetAcceptance a'' mcmcBurnInWithAutoTuning m (x : xs) a = do   a' <- mcmcIterate m x a-  a'' <- mcmcAutotune x a'+  a'' <- mcmcAutotune NormalTuningStep x a'   logDebugB $ aSummarizeCycle m a''   logDebugS $ "Acceptance rates calculated over the last " <> show x <> " iterations."   logDebugB $ aStdMonitorHeader a''

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor -- Description :  Monitor a Markov chain--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Monitor.Log -- Description :  Monitor logarithmic values--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor.Parameter -- Description :  Monitor parameters--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor.ParameterBatch -- Description :  Batch monitor parameters--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Monitor.Time -- Description :  Print time related values--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Posterior -- Description :  Types for posterior values and functions--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -3,7 +3,7 @@ -- | -- Module      :  Prior -- Description :  Types and convenience functions for computing priors--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -46,12 +46,19 @@ import Control.Monad import Data.Maybe (fromMaybe) import Data.Typeable-import Mcmc.Internal.Gamma+import Mcmc.Internal.SpecFunctions import Mcmc.Statistics.Types import Numeric.Log-import qualified Statistics.Distribution as S-import qualified Statistics.Distribution.Poisson as S +-- TODO (high): Think about using a "structure" variable.+--+-- For example,+--+-- type PriorFunctionG s a = s a -> PriorG a+--+-- Many of the prior functions would need to use the Identity Functor. This may+-- be slow.+ -- | Prior values are stored in log domain. type Prior = PriorG Double @@ -90,7 +97,7 @@  -- | Improper uniform prior; strictly less than zero. negative :: RealFloat a => PriorFunctionG a a-negative = lessThan 0+negative = lessThan 0.0 {-# SPECIALIZE negative :: PriorFunction Double #-}  -- | Exponential distributed prior.@@ -98,8 +105,8 @@ -- Call 'error' if the rate is zero or negative. exponential :: RealFloat a => Rate a -> PriorFunctionG a a exponential l x-  | l <= 0 = error "exponential: Rate is zero or negative."-  | x < 0 = error "exponential: Negative value."+  | l <= 0.0 = error "exponential: Rate is zero or negative."+  | x <= 0.0 = 0.0   | otherwise = ll * Exp (negate l * x)   where     ll = Exp $ log l@@ -110,11 +117,10 @@ -- Call 'error' if the shape or scale are zero or negative. gamma :: (Typeable a, RealFloat a) => Shape a -> Scale a -> PriorFunctionG a a gamma k t x-  | k <= 0 = error "gamma: Shape is zero or negative."-  | t <= 0 = error "gamma: Scale is zero or negative."-  | x < 0 = error "gamma: Negative value."-  | x == 0 = 0.0-  | otherwise = Exp $ log x * (k - 1) - (x / t) - logGammaG k - log t * k+  | k <= 0.0 = error "gamma: Shape is zero or negative."+  | t <= 0.0 = error "gamma: Scale is zero or negative."+  | x <= 0.0 = 0.0+  | otherwise = Exp $ log x * (k - 1.0) - (x / t) - logGammaG k - log t * k {-# SPECIALIZE gamma :: Double -> Double -> PriorFunction Double #-}  -- | See 'gamma' but parametrized using mean and variance.@@ -167,12 +173,12 @@ -- Call 'error' if the standard deviation is zero or negative. logNormal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a logNormal m s x-  | s <= 0 = error "logNormal: Standard deviation is zero or negative."-  | x <= 0 = 0+  | s <= 0.0 = error "logNormal: Standard deviation is zero or negative."+  | x <= 0.0 = 0.0   | otherwise = Exp $ t + e   where     t = negate $ mLnSqrt2Pi + log (x * s)-    a = recip $ 2 * s * s+    a = recip $ 2.0 * s * s     b = log x - m     e = negate $ a * b * b @@ -202,13 +208,11 @@ -- | Poisson distributed prior. -- -- Call 'error' if the rate is zero or negative.-poisson :: Rate Double -> PriorFunction Int+poisson :: (RealFloat a, Typeable a) => Rate a -> PriorFunctionG Int a poisson l n-  | l <= 0 = error "poisson: Rate is zero or negative."-  | n < 0 = error "poisson: Negative value."-  | otherwise = Exp $ S.logProbability d n-  where-    d = S.poisson l+  | l < 0.0 = error "poisson: Rate is zero or negative."+  | n < 0 = 0.0+  | otherwise = Exp $ log l * fromIntegral n - logFactorialG n - l  -- | Intelligent product that stops when encountering a zero. --@@ -220,5 +224,5 @@  -- The type could be generalized to any MonadPlus Integer prodM :: RealFloat a => [Log a] -> Maybe (Log a)-prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0)) 1+prodM = foldM (\ !acc x -> (acc * x) <$ guard (acc /= 0.0)) 1.0 {-# SPECIALIZE prodM :: [Log Double] -> Maybe (Log Double) #-}

@@ -5,7 +5,7 @@ -- | -- Module      :  Mcmc.Proposal -- Description :  Proposals are instruction to move around the state space--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -23,24 +23,25 @@     PSpeed (..),     Proposal (..),     KernelRatio,+    PResult (..),     Jacobian,     JacobianFunction,     (@~),     liftProposal,     liftProposalWith,-    ProposalSimple,+    PFunction,     createProposal,      -- * Tuners     Tuner (..),     Tune (..),     TuningParameter,+    TuningType (..),     TuningFunction,     AuxiliaryTuningParameters,-    AuxiliaryTuningFunction,-    defaultTuningFunctionWith,-    noTuningFunction,-    noAuxiliaryTuningFunction,+    tuningFunction,+    tuningFunctionWithAux,+    tuningFunctionOnlyAux,     tuningParameterMin,     tuningParameterMax,     tuneWithTuningParameters,@@ -52,6 +53,7 @@   ) where +import Data.Bifunctor import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import qualified Data.Double.Conversion.ByteString as BC@@ -61,8 +63,9 @@ import Lens.Micro import Mcmc.Acceptance import Mcmc.Internal.ByteString+import Mcmc.Jacobian import Numeric.Log hiding (sum)-import System.Random.MWC+import System.Random.Stateful  -- | Proposal name. newtype PName = PName {fromPName :: String}@@ -164,8 +167,8 @@     -- | The weight determines how often a 'Proposal' is executed per iteration of     -- the Markov chain.     prWeight :: PWeight,-    -- | Simple proposal without name, weight, and tuning information.-    prSimple :: ProposalSimple a,+    -- | Simple proposal function without name, weight, and tuning information.+    prFunction :: PFunction a,     -- | Tuning is disabled if set to 'Nothing'.     prTuner :: Maybe (Tuner a)   }@@ -178,24 +181,29 @@  -- | Ratio of the proposal kernels. ----- Part of the MHG acceptance ratio.------ See also 'Jacobian'.+-- For unbiased, volume preserving proposals, the values is 1.0. ----- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,--- it is more declarative to have them separate. Like so, we are constantly--- reminded: Is the Jacobian modifier different from 1.0?+-- For biased proposals, the kernel ratio is qYX / qXY, where qAB is the+-- probability density to move from A to B. type KernelRatio = Log Double --- | Absolute value of the determinant of the Jacobian matrix.------ Part of the MHG acceptance ratio.------ See also 'Jacobian'.-type Jacobian = Log Double---- | Function calculating the 'Jacobian' of a proposal.-type JacobianFunction a = a -> Jacobian+-- | Proposal result.+data PResult a+  = -- | Accept the new value regardless of the prior, likelihood or Jacobian.+    ForceAccept !a+  | -- | Reject the proposal regardless of the prior, likelihood or Jacobian.+    ForceReject+  | -- | Propose a new value.+    --+    -- In order to calculate the Metropolis-Hastings-Green ratio, we need to know+    -- the ratio of the backward to forward kernels (the 'KernelRatio' or the+    -- probability masses or probability densities) and the 'Jacobian'.+    --+    -- NOTE: Actually the 'Jacobian' should be part of the 'KernelRatio'. However,+    -- it is more declarative to have them separate. Like so, we are constantly+    -- reminded: Is the Jacobian modifier different from 1.0?+    Propose !a !KernelRatio !Jacobian+  deriving (Show, Eq)  -- | Lift a proposal from one data type to another. --@@ -219,55 +227,52 @@  -- | Lift a proposal from one data type to another. ----- A function to calculate the Jacobian has to be provided (but see '(@~)').+-- A function to calculate the Jacobian has to be provided. If the Jabobian is+-- 1.0, use 'liftProposal'. -- -- For further reference, please see the [example -- @Pair@](https://github.com/dschrempf/mcmc/blob/master/mcmc-examples/Pair/Pair.hs). liftProposalWith :: JacobianFunction b -> Lens' b a -> Proposal a -> Proposal b liftProposalWith jf l (Proposal n r d p w s t) =-  Proposal n r d p w (liftProposalSimpleWith jf l s) (liftTunerWith jf l <$> t)+  Proposal n r d p w (liftPFunctionWith jf l s) (liftTunerWith jf l <$> t) --- | Simple proposal without tuning information.+-- TODO @Dominik (high, feature): Proposals should be aware of: Is this Burn in, or not?++-- | Simple proposal function without tuning information. -- -- Instruction about randomly moving from the current state to a new state, -- given some source of randomness. ----- In order to calculate the Metropolis-Hastings-Green ratio, we need to know--- the ratio of the backward to forward kernels (the 'KernelRatio' or the--- probability masses or probability densities) and the 'Jacobian'.------ 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, KernelRatio, Jacobian)+-- Maybe report acceptance counts internal to the proposal (e.g., used by+-- proposals based on Hamiltonian dynamics).+type PFunction a = a -> IOGenM StdGen -> IO (PResult a, Maybe AcceptanceCounts) --- Lift a simple proposal from one data type to another.-liftProposalSimpleWith :: JacobianFunction b -> Lens' b a -> ProposalSimple a -> ProposalSimple b-liftProposalSimpleWith jf l s = s'+-- Lift a proposal function from one data type to another.+liftPFunctionWith :: JacobianFunction b -> Lens' b a -> PFunction a -> PFunction b+liftPFunctionWith jf l s = s'   where     s' y g = do-      (x', r, j) <- s (y ^. l) g-      let y' = set l x' y-          jxy = jf y-          jyx = jf y'-          j' = j * jyx / jxy-      return (y', r, j')+      (pr, ac) <- s (y ^. l) g+      let pr' = case pr of+            ForceAccept x' -> ForceAccept $ set l x' y+            ForceReject -> ForceReject+            Propose x' r j ->+              let y' = set l x' y+                  jxy = jf y+                  jyx = jf y'+                  j' = j * jyx / jxy+               in Propose y' r j'+      pure (pr', ac)  -- | Create a proposal with a single tuning parameter. ----- Proposals with arbitrary tuning parameters have to be created manually. See+-- Proposals with auxiliary tuning parameters have to be created manually. See -- 'Tuner' for more information, and 'Mcmc.Proposal.Hamiltonian' for an example. createProposal ::   -- | Description of the proposal type and parameters.   PDescription ->-  -- | Function creating a simple proposal for a given tuning parameter.-  (TuningParameter -> ProposalSimple a) ->+  -- | Function creating a simple proposal function for a given tuning parameter.+  (TuningParameter -> PFunction a) ->   -- | Speed.   PSpeed ->   -- | Dimension.@@ -282,36 +287,36 @@ createProposal r f s d n w Tune =   Proposal n r s d w (f 1.0) (Just tuner)   where-    fT = defaultTuningFunctionWith d-    fTs = noAuxiliaryTuningFunction+    fT = tuningFunction     g t _ = Right $ f t-    tuner = Tuner 1.0 fT VU.empty fTs g+    tuner = Tuner 1.0 VU.empty False fT g createProposal r f s d n w NoTune =   Proposal n r s d w (f 1.0) Nothing  -- | Required information to tune 'Proposal's. data Tuner a = Tuner   { tTuningParameter :: TuningParameter,-    tTuningFunction :: TuningFunction,     tAuxiliaryTuningParameters :: AuxiliaryTuningParameters,-    tAuxiliaryTuningFunction :: AuxiliaryTuningFunction a,+    -- | Does the tuner require the trace over the last tuning period?+    tRequireTrace :: Bool,+    tTuningFunction :: TuningFunction a,     -- | Given the tuning parameter, and the auxiliary tuning parameters, get-    -- the tuned simple proposal.+    -- the tuned propose function.     --     -- Should return 'Left' if the vector of auxiliary tuning parameters is     -- invalid.-    tGetSimpleProposal ::+    tGetPFunction ::       TuningParameter ->       AuxiliaryTuningParameters ->-      Either String (ProposalSimple a)+      Either String (PFunction a)   }  -- Lift tuner from one data type to another. liftTunerWith :: JacobianFunction b -> Lens' b a -> Tuner a -> Tuner b-liftTunerWith jf l (Tuner p fP ps fPs g) = Tuner p fP ps fPs' g'+liftTunerWith jf l (Tuner p ps nt fP g) = Tuner p ps nt fP' g'   where-    fPs' = fPs . VB.map (^. l)-    g' x xs = liftProposalSimpleWith jf l <$> g x xs+    fP' b d r = fP b d r . fmap (VB.map (^. l))+    g' x xs = liftPFunctionWith jf l <$> g x xs  -- | Tune proposal? data Tune = Tune | NoTune@@ -323,38 +328,53 @@ -- expected acceptance rate; and vice versa. type TuningParameter = Double --- | Compute new tuning parameter from a given acceptance rate and the old--- tuning parameter.-type TuningFunction = AcceptanceRate -> TuningParameter -> TuningParameter+-- | The last tuning step may be special.+data TuningType = NormalTuningStep | LastTuningStep --- | Auxiliary tuning parameters; vector may be empty.+-- | Compute new tuning parameters.+type TuningFunction a =+  TuningType ->+  PDimension ->+  -- | Acceptance rate of last tuning period.+  AcceptanceRate ->+  -- | Trace of last tuning period. Only available when requested by proposal.+  Maybe (VB.Vector a) ->+  (TuningParameter, AuxiliaryTuningParameters) ->+  (TuningParameter, AuxiliaryTuningParameters)++-- | Auxiliary tuning parameters. -- -- Auxiliary tuning parameters are not shown in proposal summaries.+--+-- Vector may be empty. type AuxiliaryTuningParameters = VU.Vector TuningParameter --- | Compute new auxiliary tuning parameters from a given trace and the old--- auxiliary tuning parameters.-type AuxiliaryTuningFunction a = VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters+tuningFunctionSimple :: PDimension -> AcceptanceRate -> TuningParameter -> TuningParameter+tuningFunctionSimple d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t  -- | Default tuning function. ----- Subject to change.-defaultTuningFunctionWith ::-  -- Optimal acceptance rate.-  PDimension ->-  TuningFunction-defaultTuningFunctionWith d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t+-- The default tuning function only uses the acceptance rate. In particular, it+-- does not handle auxiliary tuning parameters and ignores the actual samples+-- attained during the last tuning period.+tuningFunction :: TuningFunction a+tuningFunction _ d r _ = bimap (tuningFunctionSimple d r) id --- | Do not tune.------ Useful if auxiliary tuning parameters are tuned, but not the main tuning--- parameter.-noTuningFunction :: TuningFunction-noTuningFunction _ = id+-- | Also tune auxiliary tuning parameters.+tuningFunctionWithAux ::+  -- | Auxiliary tuning function.+  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->+  TuningFunction a+tuningFunctionWithAux _ _ _ _ Nothing = error "tuningFunctionWithAux: empty trace"+tuningFunctionWithAux f b d r (Just xs) = bimap (tuningFunctionSimple d r) (f b xs) --- | Do not tune auxiliary parameters.-noAuxiliaryTuningFunction :: AuxiliaryTuningFunction a-noAuxiliaryTuningFunction _ = id+-- | Only tune auxiliary tuning parameters.+tuningFunctionOnlyAux ::+  -- | Auxiliary tuning function.+  (TuningType -> VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters) ->+  TuningFunction a+tuningFunctionOnlyAux _ _ _ _ Nothing = error "tuningFunctionOnlyAux: empty trace"+tuningFunctionOnlyAux f b _ _ (Just xs) = bimap id (f b xs)  -- IDEA: Per proposal type tuning parameter boundaries. For example, a sliding -- proposal with a large tuning parameter is not a problem. But then, if the@@ -391,14 +411,14 @@   Either String (Proposal a) tuneWithTuningParameters t ts p = case prTuner p of   Nothing -> Left "tuneWithTuningParameters: Proposal is not tunable."-  Just (Tuner _ fT _ fTs g) ->+  Just (Tuner _ _ nt fT g) ->     -- Ensure that the tuning parameter is strictly positive and well bounded.     let t' = max tuningParameterMin t         t'' = min tuningParameterMax t'         psE = g t'' ts      in case psE of           Left err -> Left $ "tune: " <> err-          Right ps -> Right $ p {prSimple = ps, prTuner = Just $ Tuner t'' fT ts fTs g}+          Right ps -> Right $ p {prFunction = ps, prTuner = Just $ Tuner t'' ts nt fT g}  -- | See 'PDimension'. getOptimalRate :: PDimension -> Double

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Proposal.Bactrian -- Description :  Bactrian proposals--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -23,8 +23,8 @@ import Numeric.Log import Statistics.Distribution import Statistics.Distribution.Normal-import System.Random.MWC import System.Random.MWC.Distributions+import System.Random.Stateful  -- | Type synonym indicating the spike parameter. type SpikeParameter = Double@@ -32,7 +32,7 @@ genBactrian ::   SpikeParameter ->   StandardDeviation Double ->-  GenIO ->+  IOGenM StdGen ->   IO Double genBactrian m s g = do   let mn = m * s@@ -55,21 +55,20 @@ bactrianAdditive ::   SpikeParameter ->   StandardDeviation Double ->-  ProposalSimple Double+  PFunction Double bactrianAdditive m s x g = do   dx <- genBactrian m s g-  return (x + dx, 1.0, 1.0)+  pure (Propose (x + dx) 1.0 1.0, Nothing) --- bactrianSimple lens spike stdDev tune forwardOp backwardOp-bactrianAdditiveSimple ::+bactrianAdditivePFunction ::   SpikeParameter ->   StandardDeviation Double ->   TuningParameter ->-  ProposalSimple Double-bactrianAdditiveSimple m s t-  | m < 0 = error "bactrianAdditiveSimple: Spike parameter negative."-  | m >= 1 = error "bactrianAdditiveSimple: Spike parameter 1.0 or larger."-  | s <= 0 = error "bactrianAdditiveSimple: Standard deviation 0.0 or smaller."+  PFunction Double+bactrianAdditivePFunction m s t+  | m < 0 = error "bactrianAdditivePFunction: Spike parameter negative."+  | m >= 1 = error "bactrianAdditivePFunction: Spike parameter 1.0 or larger."+  | s <= 0 = error "bactrianAdditivePFunction: Standard deviation 0.0 or smaller."   | otherwise = bactrianAdditive m (t * s)  -- | Additive symmetric proposal with kernel similar to the silhouette of a@@ -90,7 +89,7 @@   PWeight ->   Tune ->   Proposal Double-slideBactrian m s = createProposal description (bactrianAdditiveSimple m s) PFast (PDimension 1)+slideBactrian m s = createProposal description (bactrianAdditivePFunction m s) PFast (PDimension 1)   where     description = PDescription $ "Slide Bactrian; spike: " ++ show m ++ ", sd: " ++ show s @@ -106,24 +105,24 @@ bactrianMult ::   SpikeParameter ->   StandardDeviation Double ->-  ProposalSimple Double+  PFunction Double bactrianMult m s x g = do   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)+  pure (Propose (x * u) (qYX / qXY) jac, Nothing) -bactrianMultSimple ::+bactrianMultPFunction ::   SpikeParameter ->   StandardDeviation Double ->   TuningParameter ->-  ProposalSimple Double-bactrianMultSimple m s t-  | m < 0 = error "bactrianMultSimple: Spike parameter negative."-  | m >= 1 = error "bactrianMultSimple: Spike parameter 1.0 or larger."-  | s <= 0 = error "bactrianMultSimple: Standard deviation 0.0 or smaller."+  PFunction Double+bactrianMultPFunction m s t+  | m < 0 = error "bactrianMultPFunction: Spike parameter negative."+  | m >= 1 = error "bactrianMultPFunction: Spike parameter 1.0 or larger."+  | s <= 0 = error "bactrianMultPFunction: Standard deviation 0.0 or smaller."   | otherwise = bactrianMult m (t * s)  -- | Multiplicative proposal with kernel similar to the silhouette of a Bactrian@@ -137,6 +136,6 @@   PWeight ->   Tune ->   Proposal Double-scaleBactrian m s = createProposal description (bactrianMultSimple m s) PFast (PDimension 1)+scaleBactrian m s = createProposal description (bactrianMultPFunction m s) PFast (PDimension 1)   where     description = PDescription $ "Scale Bactrian; spike: " ++ show m <> ", sd: " <> show s

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Proposal.Generic -- Description :  Generic interface for creating proposals--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -28,7 +28,7 @@ -- 1. Let \(D\) be a continuous probability distribution on \(\mathbb{D}\); --    sample an auxiliary variable \(epsilon \sim D\). ----- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). Propose a+-- 2. Suppose \(\odot : \mathbb{X} \times \mathbb{D} \to \mathbb{X}\). PFunction a --    new state \(x' = x \odot \epsilon\). -- -- 3. If the proposal is unbiased, the Metropolis-Hastings-Green ratio can@@ -71,7 +71,7 @@   -- Conversion to log domain is necessary, because some determinants of   -- Jacobians are very small (or large).   Maybe (a -> Double -> Jacobian) ->-  ProposalSimple a+  PFunction a genericContinuous d f mInv mJac x g = do   u <- genContVar d g   let r = case mInv of@@ -83,7 +83,7 @@       j = case mJac of         Nothing -> 1.0         Just fJac -> fJac x u-  return (x `f` u, r, j)+  pure (Propose (x `f` u) r j, Nothing) {-# INLINEABLE genericContinuous #-}  -- | Generic function to create proposals for discrete parameters (e.g., 'Int').@@ -103,7 +103,7 @@   --   -- Only required for biased proposals.   Maybe (Int -> Int) ->-  ProposalSimple a+  PFunction a genericDiscrete d f mfInv x g = do   u <- genDiscreteVar d g   let r = case mfInv of@@ -112,5 +112,5 @@           let qXY = Exp $ logProbability d u               qYX = Exp $ logProbability d (fInv u)            in qYX / qXY-  return (x `f` u, r, 1.0)+  pure (Propose (x `f` u) r 1.0, Nothing) {-# INLINEABLE genericDiscrete #-}

@@ -1,601 +0,0 @@-{-# LANGUAGE TupleSections #-}-{-# LANGUAGE UndecidableInstances #-}---- |--- Module      :  Mcmc.Proposal.Hamiltonian--- Description :  Hamiltonian Monte Carlo proposal--- Copyright   :  (c) 2021 Dominik Schrempf--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  experimental--- Portability :  portable------ Creation date: Mon Jul  5 12:59:42 2021.------ The Hamiltonian Monte Carlo (HMC) proposal.------ For references, see:------ - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using---   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),---   Handbook of Markov Chain Monte Carlo (2011), CRC press.------ - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data---   analysis (2014), CRC Press.------ - [3] Review by Betancourt and notes: Betancourt, M., A conceptual---   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).------ NOTE on implementation:------ - The implementation assumes the existence of the 'Gradient'. Like so, the---   user can use automatic or manual differentiation, depending on the problem---   at hand.------ - The Hamiltonian proposal acts on a vector of storable 'Values'. Functions---   converting the state to and from this vector have to be provided. See---   'HSettings'.------ - The desired acceptance rate is 0.65, although the dimension of the proposal---   is high.------ - The speed of this proposal can change drastically when tuned because the---   leapfrog trajectory length is changed.------ - The Hamiltonian proposal is agnostic of the actual prior and likelihood---   functions, and so, points with zero posterior probability cannot be---   detected. This affects models with constrained parameters. See Gelman p.---   303. This problem can be ameliorated by providing a 'Validate' function so---   that the proposal can gracefully fail as soon as the state becomes invalid.-module Mcmc.Proposal.Hamiltonian-  ( Values,-    Gradient,-    Validate,-    Masses,-    LeapfrogTrajectoryLength,-    LeapfrogScalingFactor,-    HTuneLeapfrog (..),-    HTuneMasses (..),-    HTune (..),-    HSettings (..),-    hamiltonian,-  )-where--import qualified Data.Vector as VB-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Unboxed as VU-import Mcmc.Proposal-import qualified Numeric.LinearAlgebra as L-import Numeric.Log-import Numeric.MathFunctions.Constants-import qualified Statistics.Covariance as S-import qualified Statistics.Function as S-import qualified Statistics.Sample as S-import System.Random.MWC---- TODO: No-U-turn sampler (NUTS). Ameliorates necessity to determine the--- leapfrog trajectory length L. (I think this is a necessary extension.)---- TODO: Riemannian adaptation: State-dependent mass matrix. (Seems a little bit--- of an overkill.)---- | The Hamiltonian proposal acts on a vector of floating point values.-type Values = L.Vector Double---- | Gradient of the log posterior function.------ The gradient has to be provided for the complete state. The reason is that--- the gradient may change if parameters untouched by the Hamiltonian proposal--- are altered by other proposals.-type Gradient a = a -> a---- | Function validating the state.------ Useful if parameters are constrained.------ Also the validity of the state may depend on parameters untouched by the--- Hamiltonian proposal.-type Validate a = a -> Bool---- | Parameter mass matrix.------ The masses roughly describe how reluctant the particles move through the--- state space. If a parameter has higher mass, the momentum in this direction--- will be changed less by the provided gradient, than when the same parameter--- has lower mass. Off-diagonal entries describe the covariance structure. If--- two parameters are negatively correlated, their generated initial momenta are--- likely to have opposite signs.------ The proposal is more efficient if masses are assigned according to the--- inverse (co)-variance structure of the posterior function. That is,--- parameters changing on larger scales should have lower masses than parameters--- changing on lower scales. In particular, the optimal entries of the diagonal--- of the mass matrix are the inverted variances of the parameters distributed--- according to the posterior function.------ Of course, the scales of the parameters of the posterior function are usually--- unknown. Often, it is sufficient to------ - set the diagonal entries of the mass matrix to identical values roughly---   scaled with the inverted estimated average variance of the posterior---   function; or even to------ - set all diagonal entries of the mass matrix to 1.0, and all other entries---   to 0.0, and trust the tuning algorithm (see 'HTune') to find the correct---   values.-type Masses = L.Herm Double---- | Mean leapfrog trajectory length \(L\).------ Number of leapfrog steps per proposal.------ To avoid problems with ergodicity, the actual number of leapfrog steps is--- sampled per proposal from a discrete uniform distribution over the interval--- \([\text{floor}(0.8L),\text{ceiling}(1.2L)]\).------ For a discussion of ergodicity and reasons why randomization is important,--- see [1] p. 15; also mentioned in [2] p. 304.------ Usually set to 10, but larger values may be desirable.------ NOTE: To avoid errors, the left bound has an additional hard minimum of 1,--- and the right bound is required to be larger equal than the left bound.------ NOTE: Call 'error' if value is less than 1.-type LeapfrogTrajectoryLength = Int---- | Mean of leapfrog scaling factor \(\epsilon\).------ Determines the size of each leapfrog step.------ To avoid problems with ergodicity, the actual leapfrog scaling factor is--- sampled per proposal from a continuous uniform distribution over the interval--- \((0.8\epsilon,1.2\epsilon]\).------ For a discussion of ergodicity and reasons why randomization is important,--- see [1] p. 15; also mentioned in [2] p. 304.------ Usually set such that \( L \epsilon = 1.0 \), but smaller values may be--- required if acceptance rates are low.------ NOTE: Call 'error' if value is zero or negative.-type LeapfrogScalingFactor = Double---- Internal. Values; target state containing parameters.-type Positions = Values---- Internal. Momenta of the parameters.-type Momenta = L.Vector Double---- | Tune leapfrog parameters?-data HTuneLeapfrog-  = HNoTuneLeapfrog-  | -- | We expect that the larger the leapfrog scaling factor the lower the-    -- acceptance ratio. Consequently, if the acceptance rate is too low, the-    -- leapfrog scaling factor is decreased and vice versa. Further, the leapfrog-    -- trajectory length is scaled such that the product of the leapfrog scaling-    -- factor and leapfrog trajectory length stays roughly constant.-    HTuneLeapfrog-  deriving (Eq, Show)---- | Tune masses?------ The masses are tuned according to the (co)variances of the parameters--- obtained from the posterior distribution over the last auto tuning interval.-data HTuneMasses-  = HNoTuneMasses-  | -- | Diagonal only: The variances of the parameters are calculated and the-    -- masses are amended using the old masses and the inverted variances. If, for-    -- a specific coordinate, the sample size is 60 or lower, or if the calculated-    -- variance is out of predefined bounds [1e-6, 1e6], the mass of the affected-    -- position is not changed.-    HTuneDiagonalMassesOnly-  | -- | All masses: The covariance matrix of the parameters is estimated and the-    -- inverted matrix (sometimes called precision matrix) is used as mass matrix.-    -- This procedure is error prone, but models with high correlations between-    -- parameters it is necessary to tune off-diagonal entries. The full mass-    -- matrix is only tuned if more than 200 samples are available. For these-    -- reasons, when tuning all masses it is recommended to use tuning settings-    -- such as-    ---    -- @-    -- BurnInWithCustomAutoTuning ([10, 20 .. 200] ++ replicate 5 500)-    -- @-    HTuneAllMasses-  deriving (Eq, Show)---- | Tuning settings.-data HTune = HTune HTuneLeapfrog HTuneMasses-  deriving (Eq, Show)---- | Specifications of the Hamilton Monte Carlo proposal.-data HSettings a = HSettings-  { -- | The sample state is used for error checks and to calculate the dimension-    -- of the proposal.-    hSample :: a,-    -- | Extract values to be manipulated by the Hamiltonian proposal from the-    -- state.-    hToVector :: a -> Values,-    -- | Put those values back into the state.-    hFromVectorWith :: a -> Values -> a,-    hGradient :: Gradient a,-    hMaybeValidate :: Maybe (Validate a),-    hMasses :: Masses,-    hLeapfrogTrajectoryLength :: LeapfrogTrajectoryLength,-    hLeapfrogScalingFactor :: LeapfrogScalingFactor,-    hTune :: HTune-  }--checkHSettings :: Eq a => HSettings a -> Maybe String-checkHSettings (HSettings x toVec fromVec _ _ masses l eps _)-  | any (<= 0) diagonalMasses = Just "checkHSettings: Some diagonal entries of the mass matrix are zero or negative."-  | nrows /= ncols = Just "checkHSettings: Mass matrix is not square."-  | fromVec x xVec /= x = Just "checkHSettings: 'fromVectorWith x (toVector x) /= x' for sample state."-  | L.size xVec /= nrows = Just "checkHSettings: Mass matrix has different size than 'toVector x', where x is sample state."-  | l < 1 = Just "checkHSettings: Leapfrog trajectory length is zero or negative."-  | eps <= 0 = Just "checkHSettings: Leapfrog scaling factor is zero or negative."-  | otherwise = Nothing-  where-    ms = L.unSym masses-    diagonalMasses = L.toList $ L.takeDiag ms-    nrows = L.rows ms-    ncols = L.cols ms-    xVec = toVec x---- Internal. Mean vector containing zeroes.-type HMu = L.Vector Double---- Internal. Symmetric, inverted mass matrix.-type HMassesInv = L.Herm Double---- Internal. Symmetric, inverted mass matrix scaled with the leapfrog step size--- epsilon.-type HMassesInvEps = L.Herm Double---- Internal. Logarithm of the determinant of the mass matrix.-type HLogDetMasses = Double---- Internal data type containing memoized values.-data HData = HData-  { _hMu :: HMu,-    _hMassesInv :: HMassesInv,-    _hMassesInvEps :: HMassesInvEps,-    _hLogDetMasses :: HLogDetMasses-  }---- Call 'error' if the determinant of the covariance matrix is negative.-getHData :: HSettings a -> HData-getHData s =-  -- The multivariate normal distribution requires a positive definite matrix-  -- with positive determinant.-  if sign == 1.0-    then HData mu massesInvH massesInvEpsH logDetMasses-    else-      let msg =-            "hamiltonianSimple: Determinant of covariance matrix is negative."-              <> " The logarithm of the absolute value of the determinant is: "-              <> show logDetMasses-              <> "."-       in error msg-  where-    ms = hMasses s-    nrows = L.rows $ L.unSym ms-    mu = L.fromList $ replicate nrows 0.0-    (massesInv, (logDetMasses, sign)) = L.invlndet $ L.unSym ms-    -- In theory we can trust that the matrix is symmetric here, because the-    -- inverse of a symmetric matrix is symmetric. However, one may want to-    -- implement a check anyways.-    massesInvH = L.trustSym massesInv-    eps = hLeapfrogScalingFactor s-    massesInvEpsH = L.scale eps massesInvH--generateMomenta ::-  -- Provided so that it does not have to be recreated.-  HMu ->-  Masses ->-  GenIO ->-  IO Momenta-generateMomenta mu masses gen = do-  seed <- uniformM gen :: IO Int-  let momenta = L.gaussianSample seed 1 mu masses-  return $ L.flatten momenta---- Prior distribution of momenta.------ Log of density of multivariate normal distribution with given parameters.--- https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function.-logDensityMultivariateNormal ::-  -- Mean vector.-  L.Vector Double ->-  -- Inverted covariance matrix.-  L.Herm Double ->-  -- Logarithm of the determinant of the covariance matrix.-  Double ->-  -- Value vector.-  L.Vector Double ->-  Log Double-logDensityMultivariateNormal mu sigmaInvH logDetSigma xs =-  Exp $ c + (-0.5) * (logDetSigma + ((dxs L.<# sigmaInv) L.<.> dxs))-  where-    dxs = xs - mu-    k = fromIntegral $ L.size mu-    c = negate $ m_ln_sqrt_2_pi * k-    sigmaInv = L.unSym sigmaInvH--leapfrog ::-  Gradient Positions ->-  Maybe (Validate Positions) ->-  HMassesInvEps ->-  LeapfrogTrajectoryLength ->-  LeapfrogScalingFactor ->-  Positions ->-  Momenta ->-  -- Maybe (Positions', Momenta'); fail if state is not valid.-  Maybe (Positions, Momenta)-leapfrog grad mVal hMassesInvEps l eps theta phi = do-  let -- The first half step of the momenta.-      phiHalf = leapfrogStepMomenta (0.5 * eps) grad theta phi-  -- L-1 full steps. This gives the positions theta_{L-1}, and the momenta-  -- phi_{L-1/2}.-  (thetaLM1, phiLM1Half) <- go (l - 1) (Just (theta, phiHalf))-  -- The last full step of the positions.-  thetaL <- valF $ leapfrogStepPositions hMassesInvEps thetaLM1 phiLM1Half-  let -- The last half step of the momenta.-      phiL = leapfrogStepMomenta (0.5 * eps) grad thetaL phiLM1Half-  return (thetaL, phiL)-  where-    valF x = case mVal of-      Nothing -> Just x-      Just f -> if f x then Just x else Nothing-    go _ Nothing = Nothing-    go 0 (Just (t, p)) = Just (t, p)-    go n (Just (t, p)) =-      let t' = leapfrogStepPositions hMassesInvEps t p-          p' = leapfrogStepMomenta eps grad t' p-          r = (,p') <$> valF t'-       in go (n - 1) r--leapfrogStepMomenta ::-  -- Step size (usually a multiple of the leapfrog scaling factor).-  Double ->-  Gradient Positions ->-  -- Current positions.-  Positions ->-  -- Current momenta.-  Momenta ->-  -- New momenta.-  Momenta-leapfrogStepMomenta eps grad theta phi = phi + L.scale eps (grad theta)--leapfrogStepPositions ::-  HMassesInvEps ->-  -- Current positions.-  Positions ->-  -- Current momenta.-  Momenta ->-  -- New positions.-  Positions-leapfrogStepPositions hMassesInvEps theta phi = theta + (L.unSym hMassesInvEps L.#> phi)--massesToTuningParameters :: Masses -> AuxiliaryTuningParameters-massesToTuningParameters = VB.convert . L.flatten . L.unSym--tuningParametersToMasses ::-  -- Dimension of the mass matrix.-  Int ->-  AuxiliaryTuningParameters ->-  Masses-tuningParametersToMasses d = L.trustSym . L.reshape d . VB.convert--hTuningParametersToSettings ::-  HSettings a ->-  TuningParameter ->-  AuxiliaryTuningParameters ->-  Either String (HSettings a)-hTuningParametersToSettings s t ts-  | nTsNotOK =-      Left "hTuningParametersToSettings: Auxiliary variables do not have correct dimension."-  | otherwise =-      Right $-        s-          { hMasses = msTuned,-            hLeapfrogTrajectoryLength = lTuned,-            hLeapfrogScalingFactor = eTuned-          }-  where-    ms = hMasses s-    d = L.rows $ L.unSym ms-    l = hLeapfrogTrajectoryLength s-    e = hLeapfrogScalingFactor s-    (HTune tlf tms) = hTune s-    nTsNotOK =-      let nTs = VU.length ts-       in case tms of-            HNoTuneMasses -> nTs /= 0-            _ -> nTs /= d * d-    msTuned = case tms of-      HNoTuneMasses -> ms-      _ -> tuningParametersToMasses d ts-    -- The larger epsilon, the larger the proposal step size and the lower the-    -- expected acceptance ratio.-    ---    -- Further, we roughly keep \( L * \epsilon = 1.0 \). The equation is not-    -- correct, because we pull L closer to the original value to keep the-    -- runtime somewhat acceptable.-    (lTuned, eTuned) = case tlf of-      HNoTuneLeapfrog -> (l, e)-      HTuneLeapfrog -> (ceiling $ fromIntegral l / (t ** 0.9) :: Int, t * e)--hamiltonianSimpleWithTuningParameters ::-  HSettings a ->-  TuningParameter ->-  AuxiliaryTuningParameters ->-  Either String (ProposalSimple a)-hamiltonianSimpleWithTuningParameters s t ts =-  hamiltonianSimple <$> hTuningParametersToSettings s t ts---- The inverted covariance matrix and the log determinant of the covariance--- matrix are calculated by 'hamiltonianSimple'.-hamiltonianSimpleWithMemoizedCovariance ::-  HSettings a ->-  HData ->-  ProposalSimple a-hamiltonianSimpleWithMemoizedCovariance st dt x g = do-  phi <- generateMomenta mu masses g-  lRan <- uniformR (lL, lR) g-  eRan <- uniformR (eL, eR) g-  case leapfrog gradientVec mValVec massesInvEps lRan eRan theta phi of-    Nothing -> return (x, 0.0, 1.0)-    Just (theta', phi') ->-      let -- Prior of momenta.-          prPhi = logDensityMultivariateNormal mu massesInv logDetMasses phi-          -- NOTE: Neal page 12: In order for the proposal to be in detailed-          -- balance, the momenta have to be negated before proposing the new-          -- value. This is not required here since the prior involves a-          -- multivariate normal distribution with means 0.-          prPhi' = logDensityMultivariateNormal mu massesInv logDetMasses phi'-          kernelR = prPhi' / prPhi-       in return (fromVec x theta', kernelR, 1.0)-  where-    (HSettings _ toVec fromVec gradient mVal masses l e _) = st-    theta = toVec x-    lL = maximum [1 :: Int, floor $ (0.8 :: Double) * fromIntegral l]-    lR = maximum [lL, ceiling $ (1.2 :: Double) * fromIntegral l]-    eL = 0.8 * e-    eR = 1.2 * e-    (HData mu massesInv massesInvEps logDetMasses) = dt-    -- Vectorize the gradient and validation functions.-    gradientVec = toVec . gradient . fromVec x-    mValVec = mVal >>= (\f -> return $ f . fromVec x)--hamiltonianSimple ::-  HSettings a ->-  ProposalSimple a-hamiltonianSimple s = hamiltonianSimpleWithMemoizedCovariance s hd-  where-    hd = getHData s---- If changed, also change help text of 'HTuneMasses'.-massMin :: Double-massMin = 1e-6---- If changed, also change help text of 'HTuneMasses'.-massMax :: Double-massMax = 1e6---- Minimal number of unique samples required for tuning the diagonal entries of--- the mass matrix.------ If changed, also change help text of 'HTuneMasses'.-samplesMinDiagonal :: Int-samplesMinDiagonal = 61---- Minimal number of samples required for tuning all entries of the mass matrix.------ If changed, also change help text of 'HTuneMasses'.-samplesMinAll :: Int-samplesMinAll = 201--getSampleSize :: VS.Vector Double -> Int-getSampleSize = VS.length . VS.uniq . S.gsort---- Diagonal elements are variances which are strictly positive.-getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double-getNewMassDiagonalWithRescue sampleSize massOld massEstimate-  | sampleSize < samplesMinDiagonal = massOld-  -- NaN and negative masses could be errors.-  | isNaN massEstimate = massOld-  | massEstimate <= 0 = massOld-  | massMin > massNew = massMin-  | massNew > massMax = massMax-  | otherwise = massNew-  where-    massNewSqrt = recip 3 * (sqrt massOld + 2 * sqrt massEstimate)-    massNew = massNewSqrt ** 2---- NOTE: Here, we lose time because we convert the states to vectors again,--- something that has already been done. But then, auto tuning is not a runtime--- determining factor.-tuneDiagonalMassesOnly ::-  Int ->-  (a -> Positions) ->-  AuxiliaryTuningFunction a-tuneDiagonalMassesOnly dim toVec xs ts-  -- If not enough data is available, do not tune.-  | VB.length xs < samplesMinDiagonal = ts-  | otherwise =-      -- Replace the diagonal.-      massesToTuningParameters $-        L.trustSym $ massesOld - L.diag massesDiagonalOld + L.diag massesDiagonalNew-  where-    -- xs: Each vector entry contains all parameter values of one iteration.-    -- xs': Each row contains all parameter values of one iteration.-    xs' = L.fromRows $ VB.toList $ VB.map toVec xs-    sampleSizes = VS.fromList $ map getSampleSize $ L.toColumns xs'-    massesOld = L.unSym $ tuningParametersToMasses dim ts-    massesDiagonalOld = L.takeDiag massesOld-    massesDiagonalEstimate = VS.fromList $ map (recip . S.variance) $ L.toColumns xs'-    massesDiagonalNew =-      VS.zipWith3-        getNewMassDiagonalWithRescue-        sampleSizes-        massesDiagonalOld-        massesDiagonalEstimate---- NOTE: Here, we lose time because we convert the states to vectors again,--- something that has already been done. But then, auto tuning is not a runtime--- determining factor.-tuneAllMasses ::-  Int ->-  (a -> Positions) ->-  AuxiliaryTuningFunction a-tuneAllMasses dim toVec xs ts-  -- If not enough data is available, do not tune.-  | VB.length xs < samplesMinDiagonal = ts-  -- If not enough data is available, only the diagonal masses are tuned.-  | VB.length xs < samplesMinAll = fallbackDiagonal-  | L.rank xs' /= dim = fallbackDiagonal-  | otherwise = massesToTuningParameters $ L.trustSym massesNew-  where-    fallbackDiagonal = tuneDiagonalMassesOnly dim toVec xs ts-    -- xs: Each vector entry contains all parameter values of one iteration.-    -- xs': Each row contains all parameter values of one iteration.-    xs' = L.fromRows $ VB.toList $ VB.map toVec xs-    (_, ss, xsNormalized) = S.scale xs'-    -- sigmaNormalized = L.unSym $ either error id $ S.oracleApproximatingShrinkage xsNormalized-    sigmaNormalized = L.unSym $ either error fst $ S.graphicalLasso 0.5 xsNormalized-    sigma = S.rescaleSWith ss sigmaNormalized-    massesNew = L.inv sigma---- | Hamiltonian Monte Carlo proposal.-hamiltonian ::-  Eq a =>-  HSettings a ->-  PName ->-  PWeight ->-  Proposal a-hamiltonian s n w = case checkHSettings s of-  Just err -> error err-  Nothing ->-    let desc = PDescription "Hamiltonian Monte Carlo (HMC)"-        toVec = hToVector s-        dim = (L.size $ toVec $ hSample s)-        pDim = PSpecial dim 0.65-        ts = massesToTuningParameters (hMasses s)-        ps = hamiltonianSimple s-        hamiltonianWith = Proposal n desc PSlow pDim w ps-        tSet@(HTune tlf tms) = hTune s-        tFun = case tlf of-          HNoTuneLeapfrog -> noTuningFunction-          HTuneLeapfrog -> defaultTuningFunctionWith pDim-        tFunAux = case tms of-          HNoTuneMasses -> noAuxiliaryTuningFunction-          HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly dim toVec-          HTuneAllMasses -> tuneAllMasses dim toVec-     in case tSet of-          (HTune HNoTuneLeapfrog HNoTuneMasses) -> hamiltonianWith Nothing-          _ ->-            let tuner = Tuner 1.0 tFun ts tFunAux (hamiltonianSimpleWithTuningParameters s)-             in hamiltonianWith $ Just tuner

@@ -0,0 +1,215 @@+{-# LANGUAGE RankNTypes #-}++-- |+-- Module      :  Mcmc.Proposal.Hamiltonian.Common+-- Description :  Code shared by proposals based on Hamiltonian dynamics+-- Copyright   :  2022 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Fri Jun  3 10:27:12 2022.+--+-- See "Mcmc.Proposal.Hamiltonian.Hamiltonian".+--+-- References:+--+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.+--+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data+--   analysis (2014), CRC Press.+--+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).+--+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of+--   Machine Learning Research.+module Mcmc.Proposal.Hamiltonian.Common+  ( -- * Basic types+    Positions,+    Momenta,+    Masses,+    LeapfrogTrajectoryLength,+    LeapfrogScalingFactor,+    LeapfrogSimulationLength,++    -- * Tuning+    HTuneLeapfrog (..),+    HTuneMasses (..),+    HTuningConf (..),++    -- * Structure of state+    HStructure (..),++    -- * Target function+    HTarget (..),+  )+where++import Data.Typeable+import qualified Data.Vector.Storable as VS+import Mcmc.Jacobian+import Mcmc.Likelihood+import Mcmc.Prior+import qualified Numeric.LinearAlgebra as L++-- NOTE: Implementing the Riemannian adaptation (state-dependent mass matrix).+-- seems a little bit of an overkill.++-- | The Hamiltonian proposal acts on a vector of floating point values referred+-- to as positions.+--+-- The positions can represent the complete state or a subset of the state of+-- the Markov chain; see 'HStructure'.+--+-- The length of the position vector determines the size of the squared mass+-- matrix; see 'Masses'.+type Positions = VS.Vector Double++-- | Momenta of the 'Positions'.+type Momenta = VS.Vector Double++-- | Mass matrix.+--+-- The masses roughly describe how reluctant the particles move through the+-- state space. If a parameter has higher mass the velocity in this direction+-- will be changed less by the provided gradient than when the same parameter+-- has lower mass. Off-diagonal entries describe the covariance structure. If+-- two parameters are negatively correlated, their generated initial momenta are+-- likely to have opposite signs.+--+-- The matrix is square with the number of rows and columns being equal to the+-- length of 'Positions'.+--+-- The proposal is more efficient if masses are assigned according to the+-- inverse (co)-variance structure of the posterior function. That is,+-- parameters changing on larger scales should have lower masses than parameters+-- changing on lower scales. In particular, the optimal entries of the diagonal+-- of the mass matrix are the inverted variances of the parameters distributed+-- according to the posterior function.+--+-- Of course, the scales of the parameters of the posterior function are usually+-- unknown. Often, it is sufficient to+--+-- - set the diagonal entries of the mass matrix to identical values roughly+--   scaled with the inverted estimated average variance of the posterior+--   function; or even to+--+-- - set all diagonal entries of the mass matrix to 1.0, and all other entries+--   to 0.0, and trust the tuning algorithm (see 'HTuningConf') to find the+--   correct values.+type Masses = L.Herm Double++-- | Mean leapfrog trajectory length \(L\).+--+-- Number of leapfrog steps per proposal.+--+-- To avoid problems with ergodicity, the actual number of leapfrog steps is+-- sampled per proposal from a discrete uniform distribution over the interval+-- \([\text{floor}(0.9L),\text{ceiling}(1.1L)]\).+--+-- For a discussion of ergodicity and reasons why randomization is important,+-- see [1] p. 15; also mentioned in [2] p. 304.+--+-- To avoid errors, the left bound of the interval has an additional hard+-- minimum of 1, and the right bound is required to be larger equal than the+-- left bound.+type LeapfrogTrajectoryLength = Int++-- | Mean of leapfrog scaling factor \(\epsilon\).+--+-- Determines the size of each leapfrog step.+--+-- To avoid problems with ergodicity, the actual leapfrog scaling factor is+-- sampled per proposal from a continuous uniform distribution over the interval+-- \((0.9\epsilon,1.1\epsilon]\).+--+-- For a discussion of ergodicity and reasons why randomization is important,+-- see [1] p. 15; also mentioned in [2] p. 304.+--+-- Call 'error' if value is zero or negative.+type LeapfrogScalingFactor = Double++-- | Product of 'LeapfrogTrajectoryLength' and 'LeapfrogScalingFactor'.+--+-- A good value is hard to find and varies between applications. For example,+-- see Figure 6 in [4].+--+-- Call 'error' if value is zero or negative.+type LeapfrogSimulationLength = Double++-- | Tune leapfrog parameters?+data HTuneLeapfrog+  = HNoTuneLeapfrog+  | -- | We expect that the larger the leapfrog scaling factor the lower the+    -- acceptance ratio. Consequently, if the acceptance rate is too low, the+    -- leapfrog scaling factor is decreased and vice versa. Further, the leapfrog+    -- trajectory length is scaled such that the product of the leapfrog scaling+    -- factor and leapfrog trajectory length stays constant.+    HTuneLeapfrog+  deriving (Eq, Show)++-- | Tune masses?+--+-- The masses are tuned according to the (co)variances of the parameters+-- obtained from the posterior distribution over the last auto tuning interval.+data HTuneMasses+  = HNoTuneMasses+  | -- | Diagonal only: The variances of the parameters are calculated and the+    -- masses are amended using the old masses and the inverted variances. If,+    -- for a specific coordinate, the sample size is 60 or lower, or if the+    -- calculated variance is out of predefined bounds [1e-8, 1e8], the mass of+    -- the affected position is not changed.+    HTuneDiagonalMassesOnly+  | -- | All masses: The covariance matrix of the parameters is estimated and+    -- the inverted matrix (sometimes called precision matrix) is used as mass+    -- matrix. This procedure is error prone, but models with high correlations+    -- between parameters strongly profit from tuning off-diagonal entries. The+    -- full mass matrix is only tuned if equal or more than "(n_masses ++    -- max(n_masses, 61)" samples are available. For these reasons, when tuning+    -- all masses it is recommended to use tuning settings such as+    --+    -- @+    -- BurnInWithCustomAutoTuning [10, 20 .. 200] [200, 220 .. 500]+    -- @+    HTuneAllMasses+  deriving (Eq, Show)++-- | Tuning configuration of the Hamilton Monte Carlo proposal.+data HTuningConf = HTuningConf HTuneLeapfrog HTuneMasses+  deriving (Eq, Show)++-- | The Hamilton Monte Carlo proposal requires information about the structure+-- of the state, which is denoted as @s@.+--+-- Please also refer to the top level module documentation of+-- "Mcmc.Proposal.Hamiltonian.Hamiltonian".+data HStructure s = HStructure+  { -- | The sample state is used for error checks and to calculate the dimension+    -- of the proposal.+    hSample :: s Double,+    -- | Extract a a subset of values to be manipulated by the Hamiltonian+    -- proposal from the complete state.+    hToVector :: s Double -> Positions,+    -- | Put those values back into the complete state.+    hFromVectorWith :: s Double -> Positions -> s Double+  }++-- | The target is composed of the prior, likelihood, and jacobian functions.+--+-- The structure of the state is denoted as @s@.+--+-- Please also refer to the top level module documentation.+data HTarget s = HTarget+  { -- | Function computing the log prior.+    hPrior :: forall a. (RealFloat a, Typeable a) => Maybe (PriorFunctionG (s a) a),+    -- | Function computing the log likelihood.+    hLikelihood :: forall a. (RealFloat a, Typeable a) => LikelihoodFunctionG (s a) a,+    -- | Function computing the log of the Jacobian.+    hJacobian :: forall a. (RealFloat a, Typeable a) => Maybe (JacobianFunctionG (s a) a)+  }

@@ -0,0 +1,210 @@+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}++-- |+-- Module      :  Mcmc.Proposal.Hamiltonian.Hamiltonian+-- Description :  Hamiltonian Monte Carlo proposal+-- Copyright   :  2021 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Mon Jul  5 12:59:42 2021.+--+-- The Hamiltonian Monte Carlo (HMC) proposal.+--+-- The HMC proposal acts on 'Positions', a vector of floating point values. The+-- manipulated values can represent the complete state, or a subset of the+-- complete state. Functions converting the state to and from this vector have+-- to be provided; see 'HStructure'.+--+-- Even though the proposal may only act on a subset of the complete state, the+-- prior, likelihood, and Jacobian functions of the complete state have to be+-- provided; see 'HTarget'. This is because parameters not manipulated by the+-- HMC proposal still influence the prior, likelihood and Jacobian functions.+--+-- The points given above have implications on how the HMC proposal is handled:+-- Do not use 'liftProposalWith', 'liftProposal', or '(@~)' with the HMC+-- proposal; instead provide proper conversion functions with 'HStructure'.+--+-- The gradient of the log target function is calculated using automatic+-- differentiation; see the excellent+-- [ad](https://hackage.haskell.org/package/ad) package.+--+-- The desired acceptance rate is 0.65, although the dimension of the proposal+-- is high.+--+-- The speed of this proposal changes drastically with the leapfrog trajectory+-- length and the leapfrog scaling factor. Hence, the speed will change during+-- burn in.+--+-- References:+--+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.+--+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data+--   analysis (2014), CRC Press.+--+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).+--+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of+--   Machine Learning Research.+module Mcmc.Proposal.Hamiltonian.Hamiltonian+  ( -- * Hamiltonian Monte Carlo proposal+    HParams (..),+    defaultHParams,+    hamiltonian,+  )+where++import Data.Bifunctor+import Mcmc.Acceptance+import Mcmc.Algorithm.MHG+import Mcmc.Proposal+import Mcmc.Proposal.Hamiltonian.Common+import Mcmc.Proposal.Hamiltonian.Internal+import Mcmc.Proposal.Hamiltonian.Masses+import Numeric.AD.Double+import qualified Numeric.LinearAlgebra as L+import Numeric.Log+import System.Random.Stateful++-- | Parameters of the Hamilton Monte Carlo proposal.+--+-- If a parameter is 'Nothing', a default value is used (see 'defaultHParams').+data HParams = HParams+  { hLeapfrogScalingFactor :: Maybe LeapfrogScalingFactor,+    hLeapfrogSimulationLength :: Maybe LeapfrogSimulationLength,+    hMasses :: Maybe Masses+  }+  deriving (Show)++-- | Default parameters.+--+-- - Estimate a reasonable leapfrog scaling factor using Algorithm 4 [4]. If all+--   fails, use 0.1.+--+-- - Leapfrog simulation length is set to 0.5.+--+-- - The mass matrix is set to the identity matrix.+defaultHParams :: HParams+defaultHParams = HParams Nothing Nothing Nothing++hamiltonianPFunctionWithTuningParameters ::+  Traversable s =>+  Dimension ->+  HStructure s ->+  (s Double -> Target) ->+  TuningParameter ->+  AuxiliaryTuningParameters ->+  Either String (PFunction (s Double))+hamiltonianPFunctionWithTuningParameters d hstruct targetWith _ ts = do+  hParamsI <- fromAuxiliaryTuningParameters d ts+  pure $ hamiltonianPFunction hParamsI hstruct targetWith++-- TODO @Dominik (high, issue): Acceptance counts. How to combine with values+-- reported here and from the NUTS sampler.++-- TODO @Dominik (high, feature): The expected acceptance counts should not be+-- calculated after burn in. Rather, the actual acceptance counts should be+-- reported. For this to work, the proposal needs to know if it is in "burn in+-- phase" or not.++-- The inverted covariance matrix and the log determinant of the covariance+-- matrix are calculated by 'hamiltonianPFunction'.+hamiltonianPFunction ::+  HParamsI ->+  HStructure s ->+  (s Double -> Target) ->+  PFunction (s Double)+hamiltonianPFunction hparamsi hstruct targetWith x g = do+  p <- generateMomenta mus ms g+  eRan <- uniformRM (eL, eR) g+  -- NOTE: The NUTS paper does not sample l since l varies naturally because+  -- of epsilon. I still think it should vary because otherwise, there may be+  -- dragons due to periodicity.+  let lM = la / eRan+      lL = maximum [1 :: Int, floor $ 0.9 * lM]+      lR = maximum [lL, ceiling $ 1.1 * lM]+  lRan <- uniformRM (lL, lR) g+  case leapfrog (targetWith x) msI lRan eRan q p of+    Nothing -> pure (ForceReject, Just $ AcceptanceCounts 0 100)+    -- Check if next state is accepted here, because the Jacobian is included in+    -- the target function. If not: pure (x, 0.0, 1.0).+    Just (q', p', prQ, prQ') -> do+      let -- Prior of momenta.+          prP = exponentialKineticEnergy msI p+          prP' = exponentialKineticEnergy msI p'+          r = prQ' * prP' / (prQ * prP)+      accept <- mhgAccept r g+      -- NOTE: For example, Neal page 12: In order for the Hamiltonian proposal+      -- to be in detailed balance, the momenta have to be negated before+      -- proposing the new value. That is, the negated momenta would guide the+      -- chain back to the previous state. However, we are only interested in+      -- the positions, and are not even storing the momenta.+      let pr = if accept then ForceAccept (fromVec x q') else ForceReject+          ar = exp $ ln r+          getCounts s = max 0 $ min 100 $ round $ s * 100+          ac =+            if ar >= 0+              then let cs = getCounts ar in AcceptanceCounts cs (100 - cs)+              else error $ "hamiltonianPFunction: Acceptance rate negative."+      pure (pr, Just ac)+  where+    (HParamsI e la ms _ _ msI mus) = hparamsi+    (HStructure _ toVec fromVec) = hstruct+    q = toVec x+    eL = 0.9 * e+    eR = 1.1 * e++-- | Hamiltonian Monte Carlo proposal.+--+-- The structure of the state is denoted as @s@.+--+-- May call 'error' during initialization.+hamiltonian ::+  Traversable s =>+  HParams ->+  HTuningConf ->+  HStructure s ->+  HTarget s ->+  PName ->+  PWeight ->+  Proposal (s Double)+hamiltonian hparams htconf hstruct htarget n w =+  let -- Misc.+      desc = PDescription "Hamiltonian Monte Carlo (HMC)"+      (HStructure sample toVec fromVec) = hstruct+      dim = L.size $ toVec sample+      -- See bottom of page 1615 in [4].+      pDim = PSpecial dim 0.65+      -- Vectorize and derive the target function.+      (HTarget mPrF lhF mJcF) = htarget+      tF y = case (mPrF, mJcF) of+        (Nothing, Nothing) -> lhF y+        (Just prF, Nothing) -> prF y * lhF y+        (Nothing, Just jcF) -> lhF y * jcF y+        (Just prF, Just jcF) -> prF y * lhF y * jcF y+      tFnG = grad' (ln . tF)+      targetWith x = bimap Exp toVec . tFnG . fromVec x+      (HParams mEps mLa mMs) = hparams+      hParamsI =+        either error id $+          hParamsIWith (targetWith sample) (toVec sample) mEps mLa mMs+      ps = hamiltonianPFunction hParamsI hstruct targetWith+      hamiltonianWith = Proposal n desc PSlow pDim w ps+      -- Tuning.+      ts = toAuxiliaryTuningParameters hParamsI+      tuner = do+        tfun <- hTuningFunctionWith dim toVec htconf+        let pfun = hamiltonianPFunctionWithTuningParameters dim hstruct targetWith+        pure $ Tuner 1.0 ts True tfun pfun+   in case checkHStructureWith (hpsMasses hParamsI) hstruct of+        Just err -> error err+        Nothing -> hamiltonianWith tuner

@@ -0,0 +1,413 @@+-- Creation date: Thu Jun  9 15:12:39 2022.+--+-- See "Mcmc.Proposal.Hamiltonian.Hamiltonian".+--+-- References:+--+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.+--+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data+--   analysis (2014), CRC Press.+--+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).+--+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of+--   Machine Learning Research.++-- |+-- Module      :  Mcmc.Proposal.Hamiltonian.Internal+-- Description :  Internal definitions related to Hamiltonian dynamics+-- Copyright   :  2022 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+module Mcmc.Proposal.Hamiltonian.Internal+  ( -- * Parameters+    HParamsI (..),+    hParamsIWith,++    -- * Tuning+    toAuxiliaryTuningParameters,+    fromAuxiliaryTuningParameters,+    findReasonableEpsilon,+    hTuningFunctionWith,++    -- * Structure of state+    checkHStructureWith,++    -- * Hamiltonian dynamics+    generateMomenta,+    exponentialKineticEnergy,++    -- * Leapfrog integrator+    Target,+    leapfrog,+  )+where++import Control.Monad+import Control.Monad.ST+import Data.Foldable+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Unboxed as VU+import Mcmc.Proposal+import Mcmc.Proposal.Hamiltonian.Common+import Mcmc.Proposal.Hamiltonian.Masses+import qualified Numeric.LinearAlgebra as L+import Numeric.Log+import System.Random.Stateful++-- Variable tuning parameters.+--+-- See Algorithm 5 or 6 in [4].+data TParamsVar = TParamsVar+  { -- \bar{eps} of Algorithm 5 or 6.+    tpvLeapfrogScalingFactorMean :: LeapfrogScalingFactor,+    -- H_i of Algorithm 5 or 6.+    tpvHStatistics :: Double,+    -- m of Algorithm 5 or 6.+    tpvCurrentTuningStep :: Double+  }+  deriving (Show)++tParamsVar :: TParamsVar+tParamsVar = TParamsVar 1.0 0.0 1.0++-- Fixed tuning parameters.+--+-- See Algorithm 5 and 6 in [4].+data TParamsFixed = TParamsFixed+  { tpfEps0 :: Double,+    tpfMu :: Double,+    tpfGa :: Double,+    tpfT0 :: Double,+    tpfKa :: Double+  }+  deriving (Show)++-- The default tuning parameters in [4] are:+--+--   mu = log $ 10 * eps+--   ga = 0.05+--   t0 = 10+--   ka = 0.75+--+-- However, these default tuning parameters will be off, because the authors+-- suggesting these values tune the proposal after every single iteration.+--+-- The following values are tweaked for our case, where tuning does not happen+-- after each iteration. Of course, we could tune the leapfrog parameters after+-- each generation. Even the mass parameters could be tuned each iteration when+-- the masses are estimated from more past iterations spanning many tuning+-- intervals.+--+-- NOTE: In theory, these we could expose these internal tuning parameters to+-- the user.+tParamsFixedWith :: LeapfrogScalingFactor -> TParamsFixed+tParamsFixedWith eps = TParamsFixed eps mu ga t0 ka+  where+    mu = log $ 10 * eps+    ga = 0.1+    t0 = 3+    ka = 0.5++-- All internal parameters.+data HParamsI = HParamsI+  { hpsLeapfrogScalingFactor :: LeapfrogScalingFactor,+    hpsLeapfrogSimulationLength :: LeapfrogSimulationLength,+    hpsMasses :: Masses,+    hpsTParamsVar :: TParamsVar,+    hpsTParamsFixed :: TParamsFixed,+    hpsMassesI :: MassesI,+    hpsMu :: Mu+  }+  deriving (Show)++-- NOTE: If changed, amend help text of 'defaultHParams', and 'defaultNParams'.+defaultLeapfrogScalingFactor :: LeapfrogScalingFactor+defaultLeapfrogScalingFactor = 0.1++-- NOTE: If changed, amend help text of 'defaultHParams'.+defaultLeapfrogSimulationLength :: LeapfrogSimulationLength+defaultLeapfrogSimulationLength = 0.5++-- NOTE: If changed, amend help text of 'defaultHParams'.+defaultMassesWith :: Int -> Masses+defaultMassesWith d = L.trustSym $ L.ident d++-- Instantiate all internal parameters.+hParamsIWith ::+  Target ->+  Positions ->+  Maybe LeapfrogScalingFactor ->+  Maybe LeapfrogSimulationLength ->+  Maybe Masses ->+  Either String HParamsI+hParamsIWith htarget p mEps mLa mMs = do+  d <- case VS.length p of+    0 -> eWith "Empty position vector."+    d -> Right d+  ms <- case mMs of+    Nothing -> Right $ defaultMassesWith d+    Just ms -> do+      let ms' = cleanMatrix $ L.unSym ms+          diagonalMs = L.toList $ L.takeDiag ms'+      when (any (<= 0) diagonalMs) $ eWith "Some diagonal masses are zero or negative."+      let nrows = L.rows ms'+          ncols = L.cols ms'+      when (nrows /= ncols) $ eWith "Mass matrix is not square."+      Right ms+  let msI = getMassesI ms+      mus = getMus ms+  la <- case mLa of+    Nothing -> Right defaultLeapfrogSimulationLength+    Just l+      | l <= 0 -> eWith "Leapfrog simulation length is zero or negative."+      | otherwise -> Right l+  eps <- case mEps of+    Nothing -> Right $ runST $ do+      -- NOTE: This is not random. However, I do not want to provide a generator+      -- when creating the proposal.+      g <- newSTGenM $ mkStdGen 42+      findReasonableEpsilon htarget ms p g+    Just e+      | e <= 0 -> eWith "Leapfrog scaling factor is zero or negative."+      | otherwise -> Right e+  let tParamsFixed = tParamsFixedWith eps+  pure $ HParamsI eps la ms tParamsVar tParamsFixed msI mus+  where+    eWith m = Left $ "hParamsIWith: " <> m++-- Save internal parameters.+toAuxiliaryTuningParameters :: HParamsI -> AuxiliaryTuningParameters+toAuxiliaryTuningParameters (HParamsI eps la ms tpv tpf _ _) =+  -- Put masses to the end. Like so, conversion is easier.+  VU.fromList $ eps : la : epsMean : h : m : eps0 : mu : ga : t0 : ka : msL+  where+    (TParamsVar epsMean h m) = tpv+    (TParamsFixed eps0 mu ga t0 ka) = tpf+    msL = VU.toList $ massesToVector ms++-- Load internal parameters.+fromAuxiliaryTuningParameters :: Dimension -> AuxiliaryTuningParameters -> Either String HParamsI+fromAuxiliaryTuningParameters d xs+  | (d * d) + 10 /= len = Left "fromAuxiliaryTuningParameters: Dimension mismatch."+  | fromIntegral (d * d) /= lenMs = Left "fromAuxiliaryTuningParameters: Masses dimension mismatch."+  | otherwise = case VU.toList $ VU.take 10 xs of+      [eps, la, epsMean, h, m, eps0, mu, ga, t0, ka] ->+        let tpv = TParamsVar epsMean h m+            tpf = TParamsFixed eps0 mu ga t0 ka+         in Right $ HParamsI eps la ms tpv tpf msI mus+      -- To please the exhaustive pattern match checker.+      _ -> Left "fromAuxiliaryTuningParameters: Impossible dimension mismatch."+  where+    len = VU.length xs+    msV = VU.drop 10 xs+    lenMs = VU.length msV+    ms = vectorToMasses d msV+    msI = getMassesI ms+    mus = getMus ms++-- See Algorithm 4 in [4].+findReasonableEpsilon ::+  StatefulGen g m =>+  Target ->+  Masses ->+  Positions ->+  g ->+  m LeapfrogScalingFactor+findReasonableEpsilon t ms q g = do+  p <- generateMomenta mu ms g+  case leapfrog t msI 1 eI q p of+    Nothing -> pure defaultLeapfrogScalingFactor+    Just (_, p', prQ, prQ') -> do+      let expEKin = exponentialKineticEnergy msI p+          expEKin' = exponentialKineticEnergy msI p'+          rI :: Double+          rI = exp $ ln $ prQ' * expEKin' / (prQ * expEKin)+          a :: Double+          a = if rI > 0.5 then 1 else (-1)+          go e r =+            if r ** a > 2 ** (negate a)+              then case leapfrog t msI 1 e q p of+                Nothing -> e+                Just (_, p'', _, prQ'') ->+                  let expEKin'' = exponentialKineticEnergy msI p''+                      r' :: Double+                      r' = exp $ ln $ prQ'' * expEKin'' / (prQ * expEKin)+                      e' = (2 ** a) * e+                   in go e' r'+              else e+      pure $ go eI rI+  where+    eI = 1.0+    msI = getMassesI ms+    mu = getMus ms++hTuningFunctionWith ::+  Dimension ->+  -- Conversion from value to vector.+  (a -> Positions) ->+  HTuningConf ->+  Maybe (TuningFunction a)+hTuningFunctionWith n toVec (HTuningConf lc mc) = case (lc, mc) of+  (HNoTuneLeapfrog, HNoTuneMasses) -> Nothing+  (_, _) -> Just $+    \tt pdim ar mxs (_, ts) ->+      case mxs of+        Nothing -> error "hTuningFunctionWith: empty trace"+        Just xs ->+          let (HParamsI eps la ms tpv tpf msI mus) =+                -- NOTE: Use error here, because a dimension mismatch is a serious bug.+                either error id $ fromAuxiliaryTuningParameters n ts+              (TParamsVar epsMean h m) = tpv+              (TParamsFixed eps0 mu ga t0 ka) = tpf+              (ms', msI') = case mc of+                HNoTuneMasses -> (ms, msI)+                HTuneDiagonalMassesOnly -> tuneDiagonalMassesOnly toVec xs (ms, msI)+                HTuneAllMasses -> tuneAllMasses toVec xs (ms, msI)+              (eps'', epsMean'', h'') = case lc of+                HNoTuneLeapfrog -> (eps, epsMean, h)+                HTuneLeapfrog ->+                  let delta = getOptimalRate pdim+                      c = recip $ m + t0+                      h' = (1.0 - c) * h + c * (delta - ar)+                      logEps' = mu - (sqrt m / ga) * h'+                      eps' = exp logEps'+                      mMKa = m ** (negate ka)+                      epsMean' = exp $ mMKa * logEps' + (1 - mMKa) * log epsMean+                   in (eps', epsMean', h')+              eps''' = case tt of+                NormalTuningStep -> eps''+                LastTuningStep -> epsMean''+              tpv' = TParamsVar epsMean'' h'' (m + 1.0)+           in (eps''' / eps0, toAuxiliaryTuningParameters $ HParamsI eps''' la ms' tpv' tpf msI' mus)++checkHStructureWith :: Foldable s => Masses -> HStructure s -> Maybe String+checkHStructureWith ms (HStructure x toVec fromVec)+  | toList (fromVec x xVec) /= toList x = eWith "'fromVectorWith x (toVector x) /= x' for sample state."+  | L.size xVec /= nrows = eWith "Mass matrix and 'toVector x' have different sizes for sample state."+  | otherwise = Nothing+  where+    eWith m = Just $ "checkHStructureWith: " <> m+    nrows = L.rows $ L.unSym ms+    xVec = toVec x++-- Generate momenta for a new iteration.+generateMomenta ::+  StatefulGen g m =>+  Mu ->+  Masses ->+  g ->+  m Momenta+generateMomenta mu masses gen = do+  seed <- uniformM gen+  let momenta = L.gaussianSample seed 1 mu masses+  return $ L.flatten momenta++-- Compute exponent of kinetic energy.+--+-- Use a general matrix which has special representations for diagonal and+-- sparse matrices, both of which are really useful here.+exponentialKineticEnergy ::+  MassesI ->+  Momenta ->+  Log Double+exponentialKineticEnergy msI xs =+  -- NOTE: Because of numerical errors, the following formulas exhibit different+  -- traces (although the posterior appears to be the same):+  -- - This one we cannot use with general matrices:+  --   Exp $ (-0.5) * ((xs L.#> msI) L.<.> xs)+  Exp $ (-0.5) * (xs L.<.> (msI L.!#> xs))++-- Function calculating target value and gradient.+--+-- The function acts on the subset of the state manipulated by the proposal but+-- the value and gradient have to be calculated for the complete state. The+-- reason is that parameters untouched by the Hamiltonian proposal may affect+-- the result or the gradient.+--+-- Make sure that the value is calculated lazily because many times, only the+-- gradient is required.+type Target = Positions -> (Log Double, Positions)++-- Leapfrog integrator.+leapfrog ::+  Target ->+  MassesI ->+  --+  LeapfrogTrajectoryLength ->+  LeapfrogScalingFactor ->+  --+  Positions ->+  Momenta ->+  -- | (New positions, new momenta, old target, new target).+  --+  -- Fail if state is not valid.+  Maybe (Positions, Momenta, Log Double, Log Double)+leapfrog tF msI l eps q p = do+  -- The first half step of the momenta.+  (x, pHalf) <-+    let (x, pHalf) = leapfrogStepMomenta (0.5 * eps) tF q p+     in if x > 0.0+          then Just (x, pHalf)+          else Nothing+  -- L-1 full steps for positions and momenta. This gives the positions q_{L-1},+  -- and the momenta p_{L-1/2}.+  (qLM1, pLM1Half) <- go (l - 1) $ Just $ (q, pHalf)+  -- The last full step of the positions.+  let qL = leapfrogStepPositions msI eps qLM1 pLM1Half+  -- The last half step of the momenta.+  (x', pL) <-+    let (x', pL) = leapfrogStepMomenta (0.5 * eps) tF qL pLM1Half+     in if x' > 0.0+          then Just (x', pL)+          else Nothing+  return (qL, pL, x, x')+  where+    go _ Nothing = Nothing+    go n (Just (qs, ps))+      | n <= 0 = Just (qs, ps)+      | otherwise =+          let qs' = leapfrogStepPositions msI eps qs ps+              (x, ps') = leapfrogStepMomenta eps tF qs' p+           in if x > 0.0+                then go (n - 1) $ Just $ (qs', ps')+                else Nothing++leapfrogStepMomenta ::+  LeapfrogScalingFactor ->+  Target ->+  -- Current positions.+  Positions ->+  -- Current momenta.+  Momenta ->+  -- New momenta; also return value target function to be collected at the end+  -- of the leapfrog integration.+  (Log Double, Momenta)+leapfrogStepMomenta eps tf q p = (x, p + L.scale eps g)+  where+    (x, g) = tf q++leapfrogStepPositions ::+  MassesI ->+  LeapfrogScalingFactor ->+  -- Current positions.+  Positions ->+  -- Current momenta.+  Momenta ->+  -- New positions.+  Positions+-- NOTE: Because of numerical errors, the following formulas exhibit different+-- traces (although the posterior appears to be the same):+-- 1. This one we cannot use with general matrices:+--    leapfrogStepPositions msI eps q p = q + (L.scale eps msI L.!#> p)+-- 2. This one seems to be more numerically unstable:+--    leapfrogStepPositions msI eps q p = q + L.scale eps (msI L.!#> p)+leapfrogStepPositions msI eps q p = q + (msI L.!#> L.scale eps p)

@@ -0,0 +1,329 @@+-- |+-- Module      :  Mcmc.Proposal.Hamiltonian.Masses+-- Description :  Mass matrices+-- Copyright   :  2022 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Tue Jun 14 10:09:24 2022.+module Mcmc.Proposal.Hamiltonian.Masses+  ( Mu,+    MassesI,+    toGMatrix,+    cleanMatrix,+    getMassesI,+    getMus,+    Dimension,+    vectorToMasses,+    massesToVector,++    -- * Tuning+    tuneDiagonalMassesOnly,+    tuneAllMasses,+  )+where++import Data.Maybe+import qualified Data.Vector as VB+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Unboxed as VU+import Mcmc.Proposal.Hamiltonian.Common+import qualified Numeric.LinearAlgebra as L+import qualified Statistics.Covariance as S+import qualified Statistics.Function as S+import qualified Statistics.Sample as S++-- Mean vector containing zeroes. We save this vector because it is required+-- when sampling from the multivariate normal distribution.+type Mu = L.Vector Double++-- General, symmetric, inverted mass matrix.+type MassesI = L.GMatrix++-- Purge masses and inverted masses (excluding the diagonal) strictly smaller+-- than the precision.+--+-- If changed, also change help text of 'HTuneMasses', which is indirectly+-- affected via 'massMin', and 'massMax'.+precision :: Double+precision = 1e-8++isDiag :: L.Matrix Double -> Bool+isDiag xs = abs (sumDiag - sumFull) < precision+  where+    xsAbs = L.cmap abs xs+    sumDiag = L.sumElements (L.takeDiag xsAbs)+    sumFull = L.sumElements xsAbs++-- Consider a matrix sparse if less than (5 * number of rows) elements are+-- non-zero.+isSparse :: L.Matrix Double -> Bool+isSparse xs = nNonZero < fromIntegral nMax+  where+    n = L.rows xs+    m = min 5 n+    nMax = n * m+    f x = if abs x >= precision then 1 else 0 :: Double+    xsInd = L.cmap f xs+    nNonZero = L.sumElements xsInd++toAssocMatrix :: L.Matrix Double -> L.AssocMatrix+toAssocMatrix xs+  | n /= m = error "toAssocMatrix: Matrix not square."+  | otherwise =+      [ ((i, j), e)+        | i <- [0 .. (n - 1)],+          j <- [0 .. (n - 1)],+          let e = xs `L.atIndex` (i, j),+          abs e >= precision+      ]+  where+    n = L.rows xs+    m = L.cols xs++toGMatrix :: L.Matrix Double -> L.GMatrix+toGMatrix xs+  | n == 0 || m == 0 = error "toGMatrix: Matrix empty."+  | n /= m = error "toGMatrix: Matrix not square."+  | isDiag xs = L.mkDiagR n m $ L.takeDiag xs+  | isSparse xs = L.mkSparse $ toAssocMatrix xs+  | otherwise = L.mkDense xs+  where+    n = L.rows xs+    m = L.cols xs++-- Diagonal:+-- - NaN values are set to 1.+-- - Negative values are set to 1.+-- - Small elements are set to 'precision'.+--+-- Off-diagonal:+-- - NaN values are set to 0.+-- - Elements with absolute values strictly smaller than 'precision' are purged.+--+-- We are permissive with negative and NaN values because adequate masses are+-- crucial. The Hamiltonian algorithms also work when the masses are off.+cleanMatrix :: L.Matrix Double -> L.Matrix Double+cleanMatrix xs =+  (L.diag $ L.cmap cleanDiag xsDiag) + L.cmap cleanOffDiag xsOffDiag+  where+    xsDiag = L.takeDiag xs+    cleanDiag x+      | isNaN x = 1+      | x < 0 = 1+      -- The strict comparison is important.+      | x < precision = precision+      | otherwise = x+    xsOffDiag = xs - L.diag xsDiag+    cleanOffDiag x+      | isNaN x = 0+      -- The strict comparison is important.+      | abs x < precision = 0+      | otherwise = x++getMassesI :: L.Herm Double -> L.GMatrix+getMassesI xs+  | n == 0 || m == 0 = error "getMassesI: Matrix empty."+  | n /= m = error "getMassesI: Matrix not square."+  | sign /= 1.0 = error "getMassesI: Determinant of matrix is negative."+  | otherwise = toGMatrix $ cleanMatrix $ xsI+  where+    xs' = L.unSym $ xs+    n = L.rows xs'+    m = L.cols xs'+    (xsI, (_, sign)) = L.invlndet xs'++getMus :: Masses -> L.Vector Double+getMus xs+  | n == 0 || m == 0 = error "getMu: Matrix empty."+  | n /= m = error "getMu: Matrix not square."+  | otherwise = L.fromList $ replicate n 0.0+  where+    xs' = L.unSym xs+    n = L.rows xs'+    m = L.cols xs'++-- Dimension of the proposal.+type Dimension = Int++massesToVector :: Masses -> VU.Vector Double+massesToVector = VU.convert . L.flatten . L.unSym++vectorToMasses :: Dimension -> VU.Vector Double -> Masses+vectorToMasses d = L.trustSym . L.reshape d . VU.convert++-- If changed, also change help text of 'HTuneMasses'.+massMin :: Double+massMin = precision++-- If changed, also change help text of 'HTuneMasses'.+massMax :: Double+massMax = recip precision++-- Minimal number of unique samples required for tuning the diagonal entries of+-- the mass matrix.+--+-- NOTE: If changed, also change help text of 'HTuneMasses'.+samplesDiagonalMin :: Int+samplesDiagonalMin = 61++-- Minimal number of samples required for tuning all entries of the mass matrix.+--+-- NOTE: If changed, also change help text of 'HTuneMasses'.+samplesAllMinWith :: Dimension -> Int+samplesAllMinWith d = samplesDiagonalMin + max samplesDiagonalMin d++getSampleSize :: VS.Vector Double -> Int+getSampleSize = VS.length . VS.uniq . S.gsort++getNewMassDiagonalWithRescue :: Int -> Double -> Double -> Double+getNewMassDiagonalWithRescue sampleSize massOld massEstimate+  | sampleSize < samplesDiagonalMin = massOld+  -- Be permissive with NaN and negative diagonal masses. Diagonal masses are+  -- variances which are strictly positive.+  | isNaN massEstimate = massOld+  | massEstimate <= 0 = massOld+  | massMin > massNew = massMin+  | massNew > massMax = massMax+  | otherwise = massNew+  where+    massNewSqrt = recip 3 * (sqrt massOld + 2 * sqrt massEstimate)+    massNew = massNewSqrt ** 2++-- The Cholesky decomposition, which is performed when sampling new momenta with+-- 'generateMomenta', requires a positive definite covariance matrix. The+-- Graphical Lasso algorithm finds positive definite covariance matrices, but+-- sometimes positive definiteness is violated because of numerical errors.+-- Further, when non-diagonal masses are already non-zero, the tuning of+-- diagonal masses only may violate positive definiteness.+--+-- Find the closest positive definite matrix of a given matrix.+--+-- See https://gist.github.com/fasiha/fdb5cec2054e6f1c6ae35476045a0bbd.+findClosestPositiveDefiniteMatrix :: L.Matrix Double -> L.Matrix Double+findClosestPositiveDefiniteMatrix a+  | n == 0 || m == 0 = error "findClosestPositiveDefiniteMatrix: Matrix empty."+  | n /= m = error "findClosestPositiveDefiniteMatrix: Matrix not square."+  | isPositiveDefinite a = a+  | otherwise = go a3 1+  where+    n = L.rows a+    m = L.cols a+    b = L.unSym $ L.sym a+    (_, s, v) = L.svd b+    h = (L.tr v) L.<> (L.diag s L.<> v)+    a2 = L.scale 0.5 (b + h)+    a3 = L.unSym $ L.sym a2+    isPositiveDefinite = isJust . L.mbChol . L.trustSym+    --+    i = L.ident n+    -- See https://hackage.haskell.org/package/ieee754-0.8.0/docs/src/Numeric-IEEE.html#line-177.+    eps = 2.2204460492503131e-16+    go x k+      | isPositiveDefinite x = x+      | otherwise =+          let minEig = L.minElement $ L.cmap L.realPart $ L.eigenvalues x+              nu = negate minEig * (k ** 2) + eps+              x' = x + L.scale nu i+           in go x' (k + 1)++tuneDiagonalMassesOnly ::+  -- Conversion from value to vector.+  (a -> Positions) ->+  -- Value vector.+  VB.Vector a ->+  -- Old mass matrix, and inverted mass matrix.+  (Masses, MassesI) ->+  -- new mass matrix, and inverted mass matrix.+  (Masses, MassesI)+-- NOTE: Here, we lose time because we convert the states to vectors again,+-- something that has already been done. But then, auto tuning is not a runtime+-- determining factor.+tuneDiagonalMassesOnly toVec xs (ms, msI)+  -- If not enough data is available, do not tune.+  | VB.length xs < samplesDiagonalMin = (ms, msI)+  | dimState /= dimMs = error "tuneDiagonalMassesOnly: Dimension mismatch."+  -- Replace the diagonal.+  | otherwise =+      let msDirty = msOld - L.diag msDiagonalOld + L.diag msDiagonalNew+          -- Positive definite matrices are symmetric.+          ms' = L.trustSym $ findClosestPositiveDefiniteMatrix $ cleanMatrix msDirty+          msI' = getMassesI ms'+       in (ms', msI')+  where+    -- xs: Each element contains all parameters of one iteration.+    -- xs': Each element is a vector containing all parameters changed by the+    -- proposal of one iteration.+    xs' = VB.map toVec xs+    -- xs'': Matrix with each row containing all parameter values changed by the+    -- proposal of one iteration.+    xs'' = L.fromRows $ VB.toList xs'+    -- We can safely use 'VB.head' here since the length of 'xs' must be larger+    -- than 'samplesDiagonalMin'.+    dimState = VS.length $ VB.head xs'+    sampleSizes = VS.fromList $ map getSampleSize $ L.toColumns xs''+    msOld = L.unSym ms+    dimMs = L.rows msOld+    msDiagonalOld = L.takeDiag msOld+    msDiagonalEstimate = VS.fromList $ map (recip . S.variance) $ L.toColumns xs''+    msDiagonalNew =+      VS.zipWith3+        getNewMassDiagonalWithRescue+        sampleSizes+        msDiagonalOld+        msDiagonalEstimate++-- This value was carefully tuned using the example "hamiltonian".+defaultGraphicalLassoPenalty :: Double+defaultGraphicalLassoPenalty = 0.4++tuneAllMasses ::+  -- Conversion from value to vector.+  (a -> Positions) ->+  -- Value vector.+  VB.Vector a ->+  -- Old mass matrix, and inverted mass matrix.+  (Masses, MassesI) ->+  -- New mass matrix, and inverted mass matrix.+  (Masses, MassesI)+-- NOTE: Here, we lose time because we convert the states to vectors again,+-- something that has already been done. But then, auto tuning is not a runtime+-- determining factor.+tuneAllMasses toVec xs (ms, msI)+  -- If not enough data is available, do not tune.+  | VB.length xs < samplesDiagonalMin = (ms, msI)+  -- If not enough data is available, only the diagonal masses are tuned.+  | VB.length xs < samplesAllMinWith dimMs = fallbackDiagonal+  | L.rank xs'' /= dimState = fallbackDiagonal+  | dimState /= dimMs = error "tuneAllMasses: Dimension mismatch."+  | otherwise = (msPD', msPDI')+  where+    fallbackDiagonal = tuneDiagonalMassesOnly toVec xs (ms, msI)+    -- xs: Each element contains all parameters of one iteration.+    -- xs': Each element is a vector containing all parameters changed by the+    -- proposal of one iteration.+    xs' = VB.map toVec xs+    -- xs'': Matrix with each row containing all parameter values changed by the+    -- proposal of one iteration.+    xs'' = L.fromRows $ VB.toList xs'+    -- We can safely use 'VB.head' here since the length of 'xs' must be larger+    -- than 'samplesDiagonalMin'.+    dimState = VS.length $ VB.head xs'+    dimMs = L.rows $ L.unSym ms+    (_, ss, xsNormalized) = S.scale xs''+    sigmaNormalized =+      L.unSym $+        either error fst $+          S.graphicalLasso defaultGraphicalLassoPenalty xsNormalized+    -- Sigma is the inverted mass matrix.+    sigma = S.rescaleSWith ss sigmaNormalized+    msI' = toGMatrix sigma+    -- The masses should be positive definite, but sometimes they happen to be+    -- not because of numerical errors.+    ms' = L.sym $ cleanMatrix $ L.inv sigma+    -- Positive definite matrices are symmetric.+    msPD' = L.trustSym $ findClosestPositiveDefiniteMatrix $ L.unSym ms'+    msPDI' = if L.unSym ms' == L.unSym msPD' then msI' else getMassesI msPD'

@@ -0,0 +1,305 @@+-- |+-- Module      :  Mcmc.Proposal.Hamiltonian.Nuts+-- Description :  No-U-Turn sampler (NUTS)+-- Copyright   :  2022 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Fri May 27 09:58:23 2022.+--+-- For a general introduction to Hamiltonian proposals, see+-- "Mcmc.Proposal.Hamiltonian.Hamiltonian".+--+-- This module implements the No-U-Turn Sampler (NUTS), as described in [4].+--+-- Work in progress.+--+-- References:+--+-- - [1] Chapter 5 of Handbook of Monte Carlo: Neal, R. M., MCMC Using+--   Hamiltonian Dynamics, In S. Brooks, A. Gelman, G. Jones, & X. Meng (Eds.),+--   Handbook of Markov Chain Monte Carlo (2011), CRC press.+--+-- - [2] Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B., Bayesian data+--   analysis (2014), CRC Press.+--+-- - [3] Review by Betancourt and notes: Betancourt, M., A conceptual+--   introduction to Hamiltonian Monte Carlo, arXiv, 1701–02434 (2017).+--+-- - [4] Matthew D. Hoffman, Andrew Gelman (2014) The No-U-Turn Sampler:+--   Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, Journal of+--   Machine Learning Research.+module Mcmc.Proposal.Hamiltonian.Nuts+  ( NParams (..),+    defaultNParams,+    nuts,+  )+where++import Data.Bifunctor+import Mcmc.Acceptance+import Mcmc.Proposal+import Mcmc.Proposal.Hamiltonian.Common+import Mcmc.Proposal.Hamiltonian.Internal+import Mcmc.Proposal.Hamiltonian.Masses+import Numeric.AD.Double+import qualified Numeric.LinearAlgebra as L+import Numeric.Log+import System.Random.Stateful++-- Internal; Slice variable 'u'.+type SliceVariable = Log Double++-- Internal; Forward is True.+type Direction = Bool++-- Internal; Doubling step number 'j'.+type DoublingStep = Int++-- Internal; Number of leapfrog steps within the slice 'n'.+type NStepsOk = Int++-- Internal; Estimated acceptance rate \(\alpha\)'.+type Alpha = Log Double++-- Internal; Number of accepted steps.+type NAlpha = Int++-- Internal; Well, that's fun, isn't it? Have a look at Algorithm 3 in [4].+type BuildTreeReturnType = (Positions, Momenta, Positions, Momenta, Positions, NStepsOk, Alpha, NAlpha)++-- Constant determining largest allowed leapfrog integration error. See+-- discussion around Equation (3) in [4].+deltaMax :: Log Double+deltaMax = Exp 1000++-- Second function in Algorithm 3 and Algorithm 6, respectively in [4].+buildTreeWith ::+  -- The exponent of the total energy of the starting state is used to+  -- calcaulate the expected acceptance rate 'Alpha'.+  Log Double ->+  MassesI ->+  Target ->+  IOGenM StdGen ->+  --+  Positions ->+  Momenta ->+  SliceVariable ->+  Direction ->+  DoublingStep ->+  LeapfrogScalingFactor ->+  IO (Maybe BuildTreeReturnType)+buildTreeWith expETot0 msI tfun g q p u v j e+  | j <= 0 =+      -- Move backwards or forwards?+      let e' = if v then e else negate e+       in case leapfrog tfun msI 1 e' q p of+            Nothing -> pure Nothing+            Just (q', p', _, expEPot') ->+              if errorIsSmall+                then pure $ Just (q', p', q', p', q', n', alpha, 1)+                else pure Nothing+              where+                expEKin' = exponentialKineticEnergy msI p'+                expETot' = expEPot' * expEKin'+                n' = if u <= expEPot' * expEKin' then 1 else 0+                errorIsSmall = u < deltaMax * expETot'+                alpha' = expETot' / expETot0+                alpha = min 1.0 alpha'++  -- Recursive case. This is complicated because the algorithm is written for an+  -- imperative language, and because we have two stacked monads.+  | otherwise = do+      mr <- buildTree q p u v (j - 1) e+      case mr of+        Nothing -> pure Nothing+        -- Here, the suffixes 'm' and 'p' stand for minus and plus, respectively.+        Just (qm, pm, qp, pp, q', n', a', na') -> do+          mr' <-+            if v+              then -- Forwards.+              do+                mr'' <- buildTree qp pp u v (j - 1) e+                case mr'' of+                  Nothing -> pure Nothing+                  Just (_, _, qp', pp', q'', n'', a'', na'') ->+                    pure $ Just (qm, pm, qp', pp', q'', n'', a'', na'')+              else -- Backwards.+              do+                mr'' <- buildTree qm pm u v (j - 1) e+                case mr'' of+                  Nothing -> pure Nothing+                  Just (qm', pm', _, _, q'', n'', a'', na'') ->+                    pure $ Just (qm', pm', qp, pp, q'', n'', a'', na'')+          case mr' of+            Nothing -> pure Nothing+            Just (qm'', pm'', qp'', pp'', q''', n''', a''', na''') -> do+              b <- uniformRM (0, 1) g :: IO Double+              let q'''' = if b < fromIntegral n''' / (fromIntegral $ n' + n''') then q''' else q'+                  a'''' = a' + a'''+                  na'''' = na' + na'''+                  n'''' = n' + n'''+                  -- Important: Check for U-turn. This formula differs from the+                  -- formula using indicator functions in Algorithm 3. However,+                  -- check Equation (4).+                  isUTurn = let dq = (qp'' - qm'') in (dq * pm'' < 0) || (dq * pp'' < 0)+              if isUTurn+                then pure Nothing+                else pure $ Just (qm'', pm'', qp'', pp'', q'''', n'''', a'''', na'''')+  where+    buildTree = buildTreeWith expETot0 msI tfun g++-- | Paramters of the NUTS proposal.+--+-- Includes tuning parameters and tuning configuration.+data NParams = NParams+  { nLeapfrogScalingFactor :: Maybe LeapfrogScalingFactor,+    nMasses :: Maybe Masses+  }+  deriving (Show)++-- | Default parameters.+--+-- - Estimate a reasonable leapfrog scaling factor using Algorithm 4 [4]. If all+--   fails, use 0.1.+--+-- - The mass matrix is set to the identity matrix.+defaultNParams :: NParams+defaultNParams = NParams Nothing Nothing++nutsPFunctionWithTuningParameters ::+  Traversable s =>+  Dimension ->+  HStructure s ->+  (s Double -> Target) ->+  TuningParameter ->+  AuxiliaryTuningParameters ->+  Either String (PFunction (s Double))+nutsPFunctionWithTuningParameters d hstruct targetWith _ ts = do+  hParamsI <- fromAuxiliaryTuningParameters d ts+  pure $ nutsPFunction hParamsI hstruct targetWith++data IsNew+  = Old+  | OldWith {_acceptanceCountsOld :: AcceptanceCounts}+  | NewWith {_acceptanceCountsNew :: AcceptanceCounts}++-- First function in Algorithm 3.+nutsPFunction ::+  HParamsI ->+  HStructure s ->+  (s Double -> Target) ->+  PFunction (s Double)+nutsPFunction hparamsi hstruct targetWith x g = do+  p <- generateMomenta mus ms g+  uZeroOne <- uniformRM (0, 1) g :: IO Double+  -- NOTE (runtime): Here we need the target function value from the previous+  -- step. For now, I just recalculate the value, but this is, of course, slow!+  -- However, if other proposals have changed the state inbetween, we do need to+  -- recalculate this value.+  let q = toVec x+      expEPot = fst $ target q+      expEKin = exponentialKineticEnergy msI p+      expETot = expEPot * expEKin+      uZeroOneL = Exp $ log uZeroOne+      u = expETot * uZeroOneL+  let -- Recursive case. This is complicated because the algorithm is written for an+      -- imperative language, and because we have two stacked monads.+      --+      -- Here, the suffixes 'm' and 'p' stand for minus and plus, respectively.+      go qm pm qp pp j y n isNew = do+        v <- uniformM g :: IO Direction+        mr' <-+          if v+            then -- Forwards.+            do+              mr <- buildTreeWith expETot msI target g qp pp u v j e+              case mr of+                Nothing -> pure Nothing+                Just (_, _, qp', pp', y', n', a, na) -> pure $ Just (qm, pm, qp', pp', y', n', a, na)+            else -- Backwards.+            do+              mr <- buildTreeWith expETot msI target g qm pm u v j e+              case mr of+                Nothing -> pure Nothing+                Just (qm', pm', _, _, y', n', a, na) -> pure $ Just (qm', pm', qp, pp, y', n', a, na)+        case mr' of+          Nothing -> pure (y, isNew)+          Just (qm'', pm'', qp'', pp'', y'', n'', a, na) -> do+            let r = fromIntegral n'' / fromIntegral n :: Double+                ar = (exp $ ln a) / fromIntegral na :: Double+                getCounts s = max 0 $ min 100 $ round $ s * 100+                ac =+                  if ar >= 0+                    then let cs = getCounts ar in AcceptanceCounts cs (100 - cs)+                    else error $ "nutsPFunction: Acceptance rate negative."+            isAccept <-+              if r > 1.0+                then pure True+                else do+                  b <- uniformRM (0, 1) g+                  pure $ b < r+            let (y''', isNew') = if isAccept then (y'', NewWith ac) else (y, OldWith ac)+                isUTurn = let dq = (qp'' - qm'') in (dq * pm'' < 0) || (dq * pp'' < 0)+            if isUTurn+              then pure (y''', isNew')+              else go qm'' pm'' qp'' pp'' (j + 1) y''' (n + n'') isNew'+  (x', isNew) <- go q p q p 0 q 1 Old+  pure $ case isNew of+    Old -> (ForceReject, Just $ AcceptanceCounts 0 100)+    OldWith ac -> (ForceReject, Just $ ac)+    NewWith ac -> (ForceAccept $ fromVec x', Just $ ac)+  where+    (HParamsI e _ ms _ _ msI mus) = hparamsi+    (HStructure _ toVec fromVecWith) = hstruct+    fromVec = fromVecWith x+    target = targetWith x++-- | No U-turn Hamiltonian Monte Carlo sampler (NUTS).+--+-- The structure of the state is denoted as @s@.+--+-- May call 'error' during initialization.+nuts ::+  Traversable s =>+  NParams ->+  HTuningConf ->+  HStructure s ->+  HTarget s ->+  PName ->+  PWeight ->+  Proposal (s Double)+nuts nparams htconf hstruct htarget n w =+  let -- Misc.+      desc = PDescription "No U-turn sampler (NUTS)"+      (HStructure sample toVec fromVec) = hstruct+      dim = L.size $ toVec sample+      -- See bottom of page 1616 in [4].+      pDim = PSpecial dim 0.6+      -- Vectorize and derive the target function.+      (HTarget mPrF lhF mJcF) = htarget+      tF y = case (mPrF, mJcF) of+        (Nothing, Nothing) -> lhF y+        (Just prF, Nothing) -> prF y * lhF y+        (Nothing, Just jcF) -> lhF y * jcF y+        (Just prF, Just jcF) -> prF y * lhF y * jcF y+      tFnG = grad' (ln . tF)+      targetWith x = bimap Exp toVec . tFnG . fromVec x+      (NParams mEps mMs) = nparams+      hParamsI =+        either error id $+          hParamsIWith (targetWith sample) (toVec sample) mEps Nothing mMs+      ps = nutsPFunction hParamsI hstruct targetWith+      nutsWith = Proposal n desc PSlow pDim w ps+      -- Tuning.+      ts = toAuxiliaryTuningParameters hParamsI+      tuner = do+        tfun <- hTuningFunctionWith dim toVec htconf+        let pfun = nutsPFunctionWithTuningParameters dim hstruct targetWith+        pure $ Tuner 1.0 ts True tfun pfun+   in case checkHStructureWith (hpsMasses hParamsI) hstruct of+        Just err -> error err+        Nothing -> nutsWith tuner

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Proposal.Scale -- Description :  Multiplicative proposals--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -26,8 +26,8 @@  -- The actual proposal with tuning parameter. The tuning parameter does not -- change the mean.-scaleSimple :: Shape Double -> Scale Double -> TuningParameter -> ProposalSimple Double-scaleSimple k th t =+scalePFunction :: Shape Double -> Scale Double -> TuningParameter -> PFunction Double+scalePFunction k th t =   genericContinuous     (gammaDistr (k / t) (th * t))     (*)@@ -52,7 +52,7 @@   PWeight ->   Tune ->   Proposal Double-scale k th = createProposal description (scaleSimple k th) PFast (PDimension 1)+scale k th = createProposal description (scalePFunction k th) PFast (PDimension 1)   where     description = PDescription $ "Scale; shape: " ++ show k ++ ", scale: " ++ show th @@ -66,16 +66,16 @@   PWeight ->   Tune ->   Proposal Double-scaleUnbiased k = createProposal description (scaleSimple k (1 / k)) PFast (PDimension 1)+scaleUnbiased k = createProposal description (scalePFunction k (1 / k)) PFast (PDimension 1)   where     description = PDescription $ "Scale unbiased; shape: " ++ show k -scaleContrarilySimple ::+scaleContrarilyPFunction ::   Shape Double ->   Scale Double ->   TuningParameter ->-  ProposalSimple (Double, Double)-scaleContrarilySimple k th t =+  PFunction (Double, Double)+scaleContrarilyPFunction k th t =   genericContinuous     (gammaDistr (k / t) (th * t))     contra@@ -96,6 +96,6 @@   PWeight ->   Tune ->   Proposal (Double, Double)-scaleContrarily k th = createProposal description (scaleContrarilySimple k th) PFast (PDimension 2)+scaleContrarily k th = createProposal description (scaleContrarilyPFunction k th) PFast (PDimension 2)   where     description = PDescription $ "Scale contrariliy; shape: " ++ show k ++ ", scale: " ++ show th

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Proposal.Simplex -- Description :  Proposals on simplices--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -33,8 +33,6 @@ 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.@@ -100,16 +98,14 @@ -- 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 :: TuningParameter -> ProposalSimple Simplex-dirichletSimple t (SimplexUnsafe xs) g = do+dirichletPFunction :: TuningParameter -> PFunction Simplex+dirichletPFunction 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@@ -118,9 +114,7 @@         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)+  pure (Propose (SimplexUnsafe ys) r 1.0, Nothing)   where     tf = getTuningFunction t @@ -140,7 +134,7 @@ -- For high dimensional simplices, this proposal may have low acceptance rates. -- In this case, please see the coordinate wise 'beta' proposal. dirichlet :: PDimension -> PName -> PWeight -> Tune -> Proposal Simplex-dirichlet = createProposal (PDescription "Dirichlet") dirichletSimple PFast+dirichlet = createProposal (PDescription "Dirichlet") dirichletPFunction PFast  -- The tuning parameter is the inverted mean of the shape values. --@@ -149,8 +143,8 @@ -- analyses. -- -- See also the 'dirichlet' proposal.-betaSimple :: Dimension -> TuningParameter -> ProposalSimple Simplex-betaSimple i t (SimplexUnsafe xs) g = do+betaPFunction :: Dimension -> TuningParameter -> PFunction Simplex+betaPFunction 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.@@ -179,7 +173,8 @@       -- 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)+      y = either error id $ simplexFromVector ys+  pure (Propose y r jac, Nothing)   where     xI = xs V.! i     xsSum = V.sum xs@@ -206,6 +201,6 @@ -- This proposal has been assigned a dimension of 2. See the discussion at -- 'PDimension'. beta :: Dimension -> PName -> PWeight -> Tune -> Proposal Simplex-beta i = createProposal description (betaSimple i) PFast (PDimension 2)+beta i = createProposal description (betaPFunction i) PFast (PDimension 2)   where     description = PDescription $ "Beta; coordinate: " ++ show i

@@ -3,7 +3,7 @@ -- | -- Module      :  Mcmc.Proposal.Slide -- Description :  Additive proposals--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -26,8 +26,8 @@ import Statistics.Distribution.Uniform  -- The actual proposal with tuning parameter.-slideSimple :: Mean Double -> StandardDeviation Double -> TuningParameter -> ProposalSimple Double-slideSimple m s t =+slidePFunction :: Mean Double -> StandardDeviation Double -> TuningParameter -> PFunction Double+slidePFunction m s t =   genericContinuous (normalDistr m (s * t)) (+) (Just negate) Nothing  -- | Additive proposal.@@ -40,13 +40,13 @@   PWeight ->   Tune ->   Proposal Double-slide m s = createProposal description (slideSimple m s) PFast (PDimension 1)+slide m s = createProposal description (slidePFunction m s) PFast (PDimension 1)   where     description = PDescription $ "Slide; mean: " ++ show m ++ ", sd: " ++ show s  -- The actual proposal with tuning parameter.-slideSymmetricSimple :: StandardDeviation Double -> TuningParameter -> ProposalSimple Double-slideSymmetricSimple s t =+slideSymmetricPFunction :: StandardDeviation Double -> TuningParameter -> PFunction Double+slideSymmetricPFunction s t =   genericContinuous (normalDistr 0.0 (s * t)) (+) Nothing Nothing  -- | See 'slide'.@@ -60,13 +60,13 @@   PWeight ->   Tune ->   Proposal Double-slideSymmetric s = createProposal description (slideSymmetricSimple s) PFast (PDimension 1)+slideSymmetric s = createProposal description (slideSymmetricPFunction s) PFast (PDimension 1)   where     description = PDescription $ "Slide symmetric; sd: " ++ show s  -- The actual proposal with tuning parameter.-slideUniformSimple :: Size -> TuningParameter -> ProposalSimple Double-slideUniformSimple d t =+slideUniformPFunction :: Size -> TuningParameter -> PFunction Double+slideUniformPFunction d t =   genericContinuous (uniformDistr (-t * d) (t * d)) (+) Nothing Nothing  -- | See 'slide'.@@ -80,19 +80,19 @@   PWeight ->   Tune ->   Proposal Double-slideUniformSymmetric d = createProposal description (slideUniformSimple d) PFast (PDimension 1)+slideUniformSymmetric d = createProposal description (slideUniformPFunction d) PFast (PDimension 1)   where     description = PDescription $ "Slide uniform symmetric; delta: " ++ show d  contra :: (Double, Double) -> Double -> (Double, Double) contra (x, y) u = (x + u, y - u) -slideContrarilySimple ::+slideContrarilyPFunction ::   Mean Double ->   StandardDeviation Double ->   TuningParameter ->-  ProposalSimple (Double, Double)-slideContrarilySimple m s t =+  PFunction (Double, Double)+slideContrarilyPFunction m s t =   genericContinuous (normalDistr m (s * t)) contra (Just negate) Nothing  -- | See 'slide'.@@ -108,6 +108,6 @@   PWeight ->   Tune ->   Proposal (Double, Double)-slideContrarily m s = createProposal description (slideContrarilySimple m s) PFast (PDimension 2)+slideContrarily m s = createProposal description (slideContrarilyPFunction m s) PFast (PDimension 2)   where     description = PDescription $ "Slide contrarily; mean: " ++ show m ++ ", sd: " ++ show s

@@ -5,7 +5,7 @@ -- | -- Module      :  Mcmc.Settings -- Description :  Settings of Markov chain Monte Carlo samplers--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | -- Module      :  Mcmc.Statistics.Types -- Description :  Types indicating properties of distributions--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com

@@ -1,7 +1,7 @@ -- | --   Module      :  Mcmc.ProposalSpec --   Description :  Unit tests for Mcmc.Proposal---   Copyright   :  (c) Dominik Schrempf, 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -17,7 +17,7 @@ import Mcmc.Cycle import Mcmc.Proposal import Mcmc.Proposal.Slide-import System.Random.MWC+import System.Random.Stateful import Test.Hspec  p1 :: Proposal Double@@ -34,7 +34,7 @@   describe "prepareProposals" $     it "returns the correct number of proposals in a cycle" $       do-        g <- create+        g <- newIOGenM $ mkStdGen 0         l1 <- length <$> prepareProposals AllProposals c g         l1 `shouldBe` 4         l2 <- length <$> prepareProposals AllProposals (setOrder RandomReversibleO c) g

@@ -1,7 +1,7 @@ -- | --   Module      :  Mcmc.SaveSpec --   Description :  Unit tests for Mcmc.Save---   Copyright   :  (c) Dominik Schrempf, 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -14,13 +14,14 @@   ) where +import Data.IORef (readIORef) import Mcmc import Mcmc.Chain.Chain import Mcmc.Chain.Save import Mcmc.Chain.Trace import Statistics.Distribution import Statistics.Distribution.Normal-import qualified System.Random.MWC as R+import System.Random.Stateful import Test.Hspec  trueMean :: Double@@ -65,7 +66,7 @@   describe "save and load" $     it "doesn't change the MCMC chain" $       do-        gen <- R.create+        gen <- newIOGenM $ mkStdGen 0         a <- mhg settings noPrior lh proposals mon 0 gen         c <- fromMHG <$> mcmc settings a         savedChain <- toSavedChain c@@ -87,17 +88,17 @@         frozenT2 <- freezeT (trace c)         frozenT2' <- freezeT (trace c')         frozenT2 `shouldBe` frozenT2'-        g2 <- R.save $ generator r-        g2' <- R.save $ generator r'+        g2 <- readIORef $ unIOGenM $ generator r+        g2' <- readIORef $ unIOGenM $ generator r'         g2 `shouldBe` g2'    describe "mhContinue" $     it "mcmc 50 + mcmcContinue 50 == mcmc 100" $       do-        gen1 <- R.create+        gen1 <- newIOGenM $ mkStdGen 0         a1 <- mhg settings noPrior lh proposals mon 0 gen1         r1 <- fromMHG <$> mcmc settings a1-        gen2 <- R.create+        gen2 <- newIOGenM $ mkStdGen 0         let settings' = settings {sIterations = Iterations 50}         a2 <- mhg settings' noPrior lh proposals mon 0 gen2         r2' <- mcmc settings' a2@@ -107,6 +108,6 @@         frozenT1 <- freezeT (trace r1)         frozenT2 <- freezeT (trace r2)         frozenT1 `shouldBe` frozenT2-        g <- R.save $ generator r1-        g' <- R.save $ generator r2+        g <- readIORef $ unIOGenM $ generator r1+        g' <- readIORef $ unIOGenM $ generator r2         g `shouldBe` g'