diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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. diff --git a/mcmc.cabal b/mcmc.cabal --- a/mcmc.cabal +++ b/mcmc.cabal @@ -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 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 diff --git a/src/Mcmc/Chain/Save.hs b/src/Mcmc/Chain/Save.hs --- a/src/Mcmc/Chain/Save.hs +++ b/src/Mcmc/Chain/Save.hs @@ -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 ] diff --git a/src/Mcmc/Cycle.hs b/src/Mcmc/Cycle.hs --- a/src/Mcmc/Cycle.hs +++ b/src/Mcmc/Cycle.hs @@ -171,7 +171,7 @@ (prName p) (prDescription p) (prWeight p) - (tGetTuningParameter <$> prTuner p) + (tTuningParameter <$> prTuner p) (prDimension p) (ar p) | p <- ps diff --git a/src/Mcmc/Proposal.hs b/src/Mcmc/Proposal.hs --- a/src/Mcmc/Proposal.hs +++ b/src/Mcmc/Proposal.hs @@ -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 diff --git a/src/Mcmc/Proposal/Hamiltonian.hs b/src/Mcmc/Proposal/Hamiltonian.hs --- a/src/Mcmc/Proposal/Hamiltonian.hs +++ b/src/Mcmc/Proposal/Hamiltonian.hs @@ -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