mcmc 0.6.0.0 → 0.6.1.0
raw patch · 6 files changed
+434/−257 lines, 6 filesdep +covariancedep +hmatrixdep −matrices
Dependencies added: covariance, hmatrix
Dependencies removed: matrices
Files
- ChangeLog.md +7/−0
- mcmc.cabal +3/−2
- src/Mcmc/Chain/Save.hs +1/−1
- src/Mcmc/Cycle.hs +1/−1
- src/Mcmc/Proposal.hs +21/−10
- src/Mcmc/Proposal/Hamiltonian.hs +401/−243
ChangeLog.md view
@@ -5,6 +5,13 @@ ## Unreleased changes +## 0.6.1.0++- Revamp Hamiltonian proposal (storable vectors).+- Use mass matrices; allow tuning of all masses (covariance estimation using+ specialized estimators).++ ## 0.6.0.0 - Improve documentation.
mcmc.cabal view
@@ -1,6 +1,6 @@ cabal-version: 2.2 name: mcmc-version: 0.6.0.0+version: 0.6.1.0 synopsis: Sample from a posterior using Markov chain Monte Carlo description: Please see the README on GitHub at <https://github.com/dschrempf/mcmc#readme> category: Math, Statistics@@ -84,14 +84,15 @@ , bytestring , circular , containers+ , covariance , data-default , deepseq , directory , dirichlet , double-conversion+ , hmatrix , log-domain , math-functions- , matrices , microlens , mwc-random , monad-parallel
src/Mcmc/Chain/Save.hs view
@@ -72,7 +72,7 @@ ps = ccProposals cc ac' = transformKeysA ps [0 ..] ac ts =- [ (\t -> (tGetTuningParameter t, tGetAuxiliaryTuningParameters t)) <$> mt+ [ (\t -> (tTuningParameter t, tAuxiliaryTuningParameters t)) <$> mt | mt <- map prTuner ps ]
src/Mcmc/Cycle.hs view
@@ -171,7 +171,7 @@ (prName p) (prDescription p) (prWeight p)- (tGetTuningParameter <$> prTuner p)+ (tTuningParameter <$> prTuner p) (prDimension p) (ar p) | p <- ps
src/Mcmc/Proposal.hs view
@@ -30,10 +30,12 @@ ProposalSimple, Tuner (..), Tune (..),- defaultTuningFunction,- createProposal, TuningParameter, AuxiliaryTuningParameters,+ defaultTuningFunctionWith,+ noTuningFunction,+ noAuxiliaryTuningFunction,+ createProposal, tuningParameterMin, tuningParameterMax, tuneWithTuningParameters,@@ -51,6 +53,7 @@ import qualified Data.Double.Conversion.ByteString as BC import Data.Function import qualified Data.Vector as VB+import qualified Data.Vector.Unboxed as VU import Lens.Micro import Lens.Micro.Extras import Mcmc.Acceptance@@ -241,11 +244,11 @@ -- | Required information to tune 'Proposal's. data Tuner a = Tuner- { tGetTuningParameter :: TuningParameter,+ { tTuningParameter :: TuningParameter, -- | Instruction about how to compute new tuning parameter from a given -- acceptance rate and the old tuning parameter. tComputeTuningParameter :: AcceptanceRate -> TuningParameter -> TuningParameter,- tGetAuxiliaryTuningParameters :: AuxiliaryTuningParameters,+ tAuxiliaryTuningParameters :: AuxiliaryTuningParameters, -- | Instruction about how to compute new auxiliary tuning parameters from a -- given trace and the old auxiliary tuning parameters. tComputeAuxiliaryTuningParameters ::@@ -281,21 +284,29 @@ type TuningParameter = Double -- | Auxiliary tuning parameters; vector may be empty.-type AuxiliaryTuningParameters = VB.Vector TuningParameter+type AuxiliaryTuningParameters = VU.Vector TuningParameter -- | Default tuning function. -- -- Subject to change.-defaultTuningFunction ::+defaultTuningFunctionWith :: -- Optimal acceptance rate. PDimension -> AcceptanceRate -> TuningParameter -> TuningParameter-defaultTuningFunction d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t+defaultTuningFunctionWith d r t = let rO = getOptimalRate d in exp (2 * (r - rO)) * t +-- | Do no tune.+--+-- Useful if auxiliary tuning parameters are tuned, but not the main tuning+-- parameter.+noTuningFunction :: AcceptanceRate -> TuningParameter -> TuningParameter+noTuningFunction _ = id++-- | Do not tune auxiliary parameters. noAuxiliaryTuningFunction :: VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters-noAuxiliaryTuningFunction _ ts = ts+noAuxiliaryTuningFunction _ = id -- | Create a proposal with a single tuning parameter. --@@ -318,10 +329,10 @@ createProposal r f d n w Tune = Proposal n r d w (f 1.0) (Just tuner) where- fT = defaultTuningFunction d+ fT = defaultTuningFunctionWith d fTs = noAuxiliaryTuningFunction g t _ = Right $ f t- tuner = Tuner 1.0 fT VB.empty fTs g+ tuner = Tuner 1.0 fT VU.empty fTs g createProposal r f d n w NoTune = Proposal n r d w (f 1.0) Nothing
src/Mcmc/Proposal/Hamiltonian.hs view
@@ -29,105 +29,121 @@ -- -- 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 implementation assumes the existence of the 'Gradient'. Like so, the+-- user can use automatic or manual differentiation, depending on the problem+-- at hand. ----- - The state needs to be list like or 'Traversable' so that the structure of--- the state space is available. A 'Traversable' constraint on the data type--- is nice because it is more general than, for example, a list, and--- user-defined data structures can be used.+-- - The 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 state needs to have a zip-like 'Applicative' instance so that--- - matrix/vector operations can be performed.-+-- - 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- ( Gradient,+ ( Values,+ Gradient,+ Validate, Masses, LeapfrogTrajectoryLength, LeapfrogScalingFactor,+ HTuneLeapfrog (..),+ HTuneMasses (..), HTune (..), HSettings (..), hamiltonian, ) where -import Data.Foldable-import qualified Data.Matrix as M-import Data.Maybe-import Data.Traversable import qualified Data.Vector as VB-import Mcmc.Prior+import 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 Statistics.Distribution-import Statistics.Distribution.Normal+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: At the moment, the HMC proposal is agnostic of the prior and--- likelihood, that is, the posterior function. This means, that it cannot know--- when it reaches a point with zero posterior probability. This also affects--- restricted or constrained parameters. See Gelman p. 303.+-- TODO: No-U-turn sampler (NUTS). Ameliorates necessity to determine the+-- leapfrog trajectory length L. (I think this is a necessary extension.) --- TODO: No-U-turn sampler.+-- TODO: Riemannian adaptation: State-dependent mass matrix. (Seems a little bit+-- of an overkill.) --- TODO: Riemannian adaptation.+-- | The Hamiltonian proposal acts on a vector of floating point values.+type Values = L.Vector Double -- | Gradient of the log posterior function.-type Gradient f = f Double -> f Double+--+-- 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.-type Validate f = f Double -> Bool---- | Masses of parameters. ----- NOTE: Full specification of a mass matrix including off-diagonal elements is--- not supported.------ NOTE: Parameters without masses ('Nothing') are not changed by the+-- 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 particle moves through the+-- 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.+-- 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, and for a diagonal mass matrix, the--- optimal masses are the inverted variances of the parameters distributed+-- 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 masses to identical values roughly scaled with the inverted--- estimated average variance of the posterior function; or even 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 masses to 1.0, and trust the tuning algorithm (see--- 'HTuneMassesAndLeapfrog') to find the correct values.-type Masses f = f (Maybe Double)+-- - 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 proposal from a discrete uniform distribution over the interval+-- 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. ----- Usually set to 10, but larger values may be desirable.+-- NOTE: Call 'error' if value is less than 1. type LeapfrogTrajectoryLength = Int -- | Mean of leapfrog scaling factor \(\epsilon\).@@ -143,101 +159,186 @@ -- -- 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 --- Target state containing parameters.-type Positions f = f Double+-- Internal. Values; target state containing parameters.+type Positions = Values --- Momenta of the parameters.-type Momenta f = f (Maybe Double)+-- Internal. Momenta of the parameters.+type Momenta = L.Vector Double --- | Tuning settings.--------- Tuning of leapfrog parameters:------ We expect that the larger the leapfrog step size the larger the proposal step--- size and the lower the acceptance ratio. Consequently, if the acceptance rate--- is too low, the leapfrog step size is decreased and vice versa. Further, the--- leapfrog trajectory length is scaled such that the product of the leapfrog--- step size and trajectory length stays constant.------ Tuning of masses:+-- | 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 variances of all parameters of the posterior distribution obtained over--- the last auto tuning interval is calculated and the masses are amended using--- the old masses and the inverted variances. If, for a specific coordinate, the--- sample size is too low, or if the calculated variance is out of predefined--- bounds, the mass of the affected position is not changed.-data HTune- = -- | Tune masses and leapfrog parameters.- HTuneMassesAndLeapfrog- | -- | Tune leapfrog parameters only.- HTuneLeapfrogOnly- | -- | Do not tune at all.- HNoTune+-- 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) --- | Specifications for Hamilton Monte Carlo proposal.-data HSettings f = HSettings- { hGradient :: Gradient f,- hMaybeValidate :: Maybe (Validate f),- hMasses :: Masses f,+-- | Tuning settings.+data HTune = HTune HTuneLeapfrog HTuneMasses+ deriving (Eq, Show)++-- | Specifications of the Hamilton Monte Carlo proposal.+data HSettings a = HSettings+ { -- | 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 :: Foldable f => HSettings f -> Maybe String-checkHSettings (HSettings _ _ masses l eps _)- | any f masses = Just "checkHSettings: One or more masses are zero or negative."+checkHSettings :: Eq a => a -> HSettings a -> Maybe String+checkHSettings x (HSettings 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- f (Just m) = m <= 0- f Nothing = False+ 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 ::- Traversable f =>- Masses f ->+ -- Provided so that it does not have to be recreated.+ HMu ->+ Masses -> GenIO ->- IO (Momenta f)-generateMomenta masses gen = traverse (generateWith gen) masses- where- generateWith g (Just m) = let d = normalDistr 0 (sqrt m) in Just <$> genContVar d g- generateWith _ Nothing = pure Nothing+ IO Momenta+generateMomenta mu masses gen = do+ seed <- uniformM gen :: IO Int+ let momenta = L.gaussianSample seed 1 mu masses+ return $ L.flatten momenta -priorMomenta ::- (Applicative f, Foldable f) =>- Masses f ->- Momenta f ->- Prior-priorMomenta masses phi = foldl' (*) 1.0 $ f <$> masses <*> phi+-- 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- f (Just m) (Just p) = let d = normalDistr 0 (sqrt m) in Exp $ logDensity d p- f Nothing Nothing = 1.0- f _ _ = error "priorMomenta: Got just a mass and no momentum or the other way around."+ dxs = xs - mu+ k = fromIntegral $ L.size mu+ c = negate $ m_ln_sqrt_2_pi * k+ sigmaInv = L.unSym sigmaInvH leapfrog ::- Applicative f =>- Gradient f ->- Maybe (Validate f) ->- Masses f ->+ Gradient Positions ->+ Maybe (Validate Positions) ->+ HMassesInvEps -> LeapfrogTrajectoryLength -> LeapfrogScalingFactor ->- Positions f ->- Momenta f ->- -- Maybe (Positions', Momenta').- Maybe (Positions f, Momenta f)-leapfrog grad mVal masses l eps theta phi = do+ 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 eps masses thetaLM1 phiLM1Half+ 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)@@ -248,201 +349,258 @@ go _ Nothing = Nothing go 0 (Just (t, p)) = Just (t, p) go n (Just (t, p)) =- let t' = leapfrogStepPositions eps masses t p+ let t' = leapfrogStepPositions hMassesInvEps t p p' = leapfrogStepMomenta 1.0 eps grad t' p r = (,p') <$> valF t' in go (n - 1) r leapfrogStepMomenta ::- Applicative f => -- Size of step (half or full step). Double -> LeapfrogScalingFactor ->- Gradient f ->+ Gradient Positions -> -- Current positions.- Positions f ->+ Positions -> -- Current momenta.- Momenta f ->+ Momenta -> -- New momenta.- Momenta f-leapfrogStepMomenta xi eps grad theta phi = phi <+. ((xi * eps) .* grad theta)- where- (<+.) :: Applicative f => f (Maybe Double) -> f Double -> f (Maybe Double)- (<+.) xs ys = f <$> xs <*> ys- f Nothing _ = Nothing- f (Just x) y = Just $ x + y+ Momenta+leapfrogStepMomenta xi eps grad theta phi = phi + L.scale (xi * eps) (grad theta) leapfrogStepPositions ::- Applicative f =>- LeapfrogScalingFactor ->- Masses f ->+ HMassesInvEps -> -- Current positions.- Positions f ->+ Positions -> -- Current momenta.- Momenta f ->- Positions f--- The arguments are flipped to encounter the maybe momentum.-leapfrogStepPositions eps masses theta phi = theta <+. (mScaledReversed .*> phi)- where- (<+.) :: Applicative f => f Double -> f (Maybe Double) -> f Double- (<+.) xs ys = f <$> xs <*> ys- f x Nothing = x- f x (Just y) = x + y- mScaledReversed = (fmap . fmap) ((* eps) . (** (-1))) masses- (.*>) :: Applicative f => f (Maybe Double) -> f (Maybe Double) -> f (Maybe Double)- (.*>) xs ys = g <$> xs <*> ys- g (Just x) (Just y) = Just $ x * y- g Nothing Nothing = Nothing- g _ _ = error "leapfrogStepPositions: Got just a mass and no momentum or the other way around."---- Scalar-vector multiplication.-(.*) :: Applicative f => Double -> f Double -> f Double-(.*) x ys = (* x) <$> ys+ Momenta ->+ -- New positions.+ Positions+leapfrogStepPositions hMassesInvEps theta phi = theta + (L.unSym hMassesInvEps L.#> phi) --- NOTE: Fixed parameters without mass have a tuning parameter of NaN.-massesToTuningParameters :: Foldable f => Masses f -> AuxiliaryTuningParameters-massesToTuningParameters = VB.fromList . map (fromMaybe nan) . toList- where- nan = 0 / 0+massesToTuningParameters :: Masses -> AuxiliaryTuningParameters+massesToTuningParameters = VB.convert . L.flatten . L.unSym --- We need the structure in order to fill it with the given parameters. tuningParametersToMasses ::- Traversable f =>+ -- Dimension of the mass matrix.+ Int -> AuxiliaryTuningParameters ->- Masses f ->- Either String (Masses f)-tuningParametersToMasses xs ms =- if null xs'- then sequenceA msE- else Left "tuningParametersToMasses: Too many values."- where- (xs', msE) = mapAccumL setValue (VB.toList xs) ms- setValue [] _ = ([], Left "tuningParametersToMasses: Too few values.")- -- NOTE: Recover fixed parameters and unset their mass.- setValue (y : ys) _ = let y' = if isNaN y then Nothing else Just y in (ys, Right y')+ Masses+tuningParametersToMasses d = L.trustSym . L.reshape d . VB.convert hTuningParametersToSettings ::- Traversable f =>+ HSettings a -> TuningParameter -> AuxiliaryTuningParameters ->- HSettings f ->- Either String (HSettings f)-hTuningParametersToSettings t ts (HSettings g v m l e tn) =- if tn == HTuneMassesAndLeapfrog- then case tuningParametersToMasses ts m of- Left err -> Left err- Right m' -> Right $ HSettings g v m' lTuned eTuned tn- else Right $ HSettings g v m lTuned eTuned tn+ 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 = ceiling $ fromIntegral l / (t ** 0.9) :: Int- eTuned = t * e+ (lTuned, eTuned) = case tlf of+ HNoTuneLeapfrog -> (l, e)+ HTuneLeapfrog -> (ceiling $ fromIntegral l / (t ** 0.9) :: Int, t * e) hamiltonianSimpleWithTuningParameters ::- (Applicative f, Traversable f) =>- HSettings f ->+ HSettings a -> TuningParameter -> AuxiliaryTuningParameters ->- Either String (ProposalSimple (Positions f))-hamiltonianSimpleWithTuningParameters s t ts = case hTuningParametersToSettings t ts s of- Left err -> Left err- Right s' -> Right $ hamiltonianSimple s'+ Either String (ProposalSimple a)+hamiltonianSimpleWithTuningParameters s t ts =+ hamiltonianSimple <$> hTuningParametersToSettings s t ts -hamiltonianSimple ::- (Applicative f, Traversable f) =>- HSettings f ->- ProposalSimple (Positions f)-hamiltonianSimple (HSettings gradient mVal masses l e _) theta g = do- phi <- generateMomenta masses g+-- 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 gradient mVal masses lRan eRan theta phi of- Nothing -> return (theta, 0.0, 1.0)+ case leapfrog gradientVec mValVec massesInvEps lRan eRan theta phi of+ Nothing -> return (x, 0.0, 1.0) Just (theta', phi') ->- let prPhi = priorMomenta masses 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 normal distributions- -- centered around 0. However, if the multivariate normal distribution is- -- used, it makes a difference.- prPhi' = priorMomenta masses phi'+ -- 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 (theta', kernelR, 1.0)+ 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) -minVariance :: Double-minVariance = 1e-6+hamiltonianSimple ::+ HSettings a ->+ ProposalSimple a+hamiltonianSimple s = hamiltonianSimpleWithMemoizedCovariance s hd+ where+ hd = getHData s -maxVariance :: Double-maxVariance = 1e6+-- If changed, also change help text of 'HTuneMasses'.+massMin :: Double+massMin = 1e-6 -minSamples :: Int-minSamples = 60+-- If changed, also change help text of 'HTuneMasses'.+massMax :: Double+massMax = 1e6 -computeAuxiliaryTuningParameters ::- Foldable f =>- VB.Vector (Positions f) ->+-- 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) ->+ VB.Vector a -> AuxiliaryTuningParameters -> AuxiliaryTuningParameters-computeAuxiliaryTuningParameters xss ts =- VB.zipWith (\t -> rescueWith t . calcSamplesAndVariance) ts xssT+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- -- TODO: Improve matrix transposition.- xssT = VB.fromList $ M.toColumns $ M.fromLists $ VB.toList $ VB.map toList xss- calcSamplesAndVariance xs = (VB.length $ VB.uniq $ S.gsort xs, S.variance xs)- rescueWith t (sampleSize, var) =- if var < minVariance || maxVariance < var || sampleSize < minSamples- then -- then traceShow ("Rescue with " <> show t) t- t- else- let t' = sqrt (t * recip var)- in -- in traceShow ("Old mass " <> show t <> " new mass " <> show t') t'- t'+ -- 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) ->+ VB.Vector a ->+ AuxiliaryTuningParameters ->+ AuxiliaryTuningParameters+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.rescaleWith ss sigmaNormalized+ massesNew = L.inv sigma+ -- | Hamiltonian Monte Carlo proposal.------ The 'Applicative' and 'Traversable' instances are used for element-wise--- operations.------ Assume a zip-like 'Applicative' instance so that cardinality remains--- constant.------ NOTE: The desired acceptance rate is 0.65, although the dimension of the--- proposal is high.------ NOTE: The speed of this proposal can change drastically when tuned because--- the leapfrog trajectory length is changed. hamiltonian ::- (Applicative f, Traversable f) =>- -- | The sample state is used to calculate the dimension of the proposal.- f Double ->- HSettings f ->+ Eq a =>+ -- | The sample state is used for error checks and to calculate the dimension+ -- of the proposal.+ a ->+ HSettings a -> PName -> PWeight ->- Proposal (f Double)-hamiltonian x s n w = case checkHSettings s of+ Proposal a+hamiltonian x s n w = case checkHSettings x s of Just err -> error err Nothing -> let desc = PDescription "Hamiltonian Monte Carlo (HMC)"- dim = PSpecial (length x) 0.65+ toVec = hToVector s+ dim = (L.size $ toVec x)+ pDim = PSpecial dim 0.65 ts = massesToTuningParameters (hMasses s) ps = hamiltonianSimple s- p' = Proposal n desc dim w ps- fT = defaultTuningFunction dim- tS = hTune s- fTs =- if tS == HTuneMassesAndLeapfrog- then computeAuxiliaryTuningParameters- else \_ xs -> xs- in case tS of- HNoTune -> p' Nothing- _ -> p' $ Just $ Tuner 1.0 fT ts fTs (hamiltonianSimpleWithTuningParameters s)+ hamiltonianWith = Proposal n desc 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