commit 78372fbccb60397a57bc0b33dd66f03f5fc64f9f
parent ed8667057890443dae758a6e58ca00799c9b2186
Author: Jared Tobin <jared@jtobin.ca>
Date: Thu, 8 Oct 2015 08:43:11 +1300
Pick-up.
Diffstat:
6 files changed, 260 insertions(+), 153 deletions(-)
diff --git a/Examples/Example.hs b/Examples/Example.hs
@@ -17,7 +17,7 @@ rosenbrock :: Target
rosenbrock = Target lRosenbrock glRosenbrock
tunables :: Tunables
-tunables = Tunables 0.05 20
+tunables = Tunables 0.001 50
options :: Options
options = Options 5000 (fromList [1, 1])
diff --git a/Examples/Rosenbrock_HMC.hs b/Examples/Rosenbrock_HMC.hs
@@ -0,0 +1,47 @@
+import Numeric.MCMC.Hamiltonian
+import System.Random.MWC
+import System.Environment
+import System.Exit
+import System.IO
+import Numeric.AD
+import Control.Monad
+
+target :: RealFloat a => [a] -> a
+target [x0, x1] = (-1)*(5*(x1 - x0^2)^2 + 0.05*(1 - x0)^2)
+
+gTarget :: [Double] -> [Double]
+gTarget = grad target
+
+main = do
+ args <- getArgs
+ when (args == []) $ do
+ putStrLn "(hasty-hamiltonian) Rosenbrock density "
+ putStrLn "Usage: ./Rosenbrock_HMC <numSteps> <thinEvery> <nDisc> <stepSize> <inits>"
+ putStrLn " "
+ putStrLn "numSteps : Number of Markov chain iterations to run. "
+ putStrLn "thinEvery : Print every nth iteration. "
+ putStrLn "nDisc : Number of discretizing steps to take. "
+ putStrLn "stepSize : Perturbation scaling parameter. "
+ putStrLn "inits : Filepath containing points at which to "
+ putStrLn " initialize the chain. "
+ exitSuccess
+
+ inits <- fmap (map read . words) (readFile (args !! 4)) :: IO [Double]
+
+ let nepochs = read (head args) :: Int
+ thinEvery = read (args !! 1) :: Int
+ nDisc = read (args !! 2) :: Double
+ eps = read (args !! 3) :: Double
+ params = Options target gTarget nDisc eps
+ config = MarkovChain inits 0
+
+ g <- create
+ results <- runChain params nepochs thinEvery config g
+
+ hPutStrLn stderr $
+ let nAcc = accepts results
+ total = nepochs
+ in show nAcc ++ " / " ++ show total ++ " (" ++
+ show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
+ ") proposals accepted"
+
diff --git a/Examples/SPDE_HMC.hs b/Examples/SPDE_HMC.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE BangPatterns #-}
+
+import Numeric.MCMC.Hamiltonian
+import System.Random.MWC
+import System.Environment
+import System.Exit
+import System.IO
+import Control.Monad
+import Numeric.AD
+
+target :: RealFloat a => [a] -> a
+target xs = go 0 0 xs
+ where go t0 t1 [] = (- t0 / (2*h)) - (0.5 * h * t1)
+ go t0 t1 (u:us:uss) = go (t0 + (us - u)^2) (t1 + v (us + u)) uss
+ h = 1 / fromIntegral (length xs)
+ v x = (1 - x^2)^2
+{-# INLINE target #-}
+
+gTarget :: [Double] -> [Double]
+gTarget = grad target
+{-# INLINE gTarget #-}
+
+main = do
+ args <- getArgs
+ when (args == []) $ do
+ putStrLn "(hasty-hamiltonian) Stochastic partial differential equation "
+ putStrLn "Usage: ./SPDE_HMC <numSteps> <inits> <thinEvery> "
+ putStrLn " "
+ putStrLn "numSteps : Number of Markov chain iterations to run. "
+ putStrLn "thinEvery : Print every nth iteration. "
+ putStrLn "nDisc : Number of discretizing steps to take. "
+ putStrLn "stepSize : Perturbation scaling parameter. "
+ putStrLn "inits : Filepath containing points at which to "
+ putStrLn " initialize the chain. "
+ exitSuccess
+
+ inits <- fmap (map read . words) (readFile (args !! 4)) :: IO [Double]
+
+ let nepochs = read (head args) :: Int
+ thinEvery = read (args !! 1) :: Int
+ nDisc = read (args !! 2) :: Double
+ eps = read (args !! 3) :: Double
+ params = Options target gTarget nDisc eps
+ config = MarkovChain inits 0
+
+ g <- create
+ results <- runChain params nepochs thinEvery config g
+
+ hPutStrLn stderr $
+ let nAcc = accepts results
+ total = nepochs
+ in show nAcc ++ " / " ++ show total ++ " (" ++
+ show ((fromIntegral nAcc / fromIntegral total) :: Float) ++
+ ") proposals accepted"
+
diff --git a/Numeric/MCMC/Hamiltonian.hs b/Numeric/MCMC/Hamiltonian.hs
@@ -1,139 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
--- | Hamiltonian Monte Carlo. See, ex: Neal (2012)
--- http://arxiv.org/pdf/1206.1901.pdf.
-
-module Numeric.MCMC.Hamiltonian (
- Target(..)
- , Tunables(..)
- , Options(..)
- , hamiltonian
- , hmc
- ) where
-
-import Control.Applicative
-import Control.Monad
-import Control.Monad.Primitive
-import Control.Monad.Trans.Class
-import Control.Monad.Trans.State.Strict
-import Data.Vector (Vector)
-import qualified Data.Vector as V
-import System.Random.MWC
-import System.Random.MWC.Distributions
-
-
--- | The target we want to sample from. Consists of a log target and its
--- gradient.
-data Target = Target {
- lTarget :: Parameters -> Double -- ^ log target
- , glTarget :: Parameters -> Parameters -- ^ gradient of log target
- }
-
-type Parameters = Vector Double
-
-data Tunables = Tunables {
- stepSize :: !Double -- ^ step size for a given proposal
- , leapfrogs :: !Int -- ^ number of leapfrog steps to take
- } deriving Show
-
-data Options = Options {
- epochs :: !Int -- ^ number of epochs to iterate the chain
- , initial :: !Parameters -- ^ start location
- } deriving Show
-
-type Particle = (Parameters, Parameters)
-
-type Transition m = StateT Parameters m Parameters
-
--- | The Hamiltonian transition operator.
-hamiltonian
- :: PrimMonad m
- => Target
- -> Tunables
- -> Gen (PrimState m)
- -> Transition m
-hamiltonian target tunables g = do
- q0 <- get
- r0 <- V.replicateM (V.length q0) (lift $ standard g)
- zc <- lift $ uniform g
- let (q, r) = leapfrogIntegrator target tunables (q0, r0)
- next = nextState target (q0, q) (r0, r) zc
- put next
- return next
-
--- | The leapfrog or Stormer-Verlet integrator.
-leapfrogIntegrator :: Target -> Tunables -> Particle -> Particle
-leapfrogIntegrator target tunables (q0, r0) = go q0 r0 l where
- l = leapfrogs tunables
- go q r 0 = (q, r)
- go q r n =
- let (q1, r1) = leapfrog target tunables (q, r)
- in go q1 r1 (pred n)
-
--- | A single iteration of the leapfrog integrator.
-leapfrog :: Target -> Tunables -> Particle -> Particle
-leapfrog target tunables (q, r) = (qf, rf) where
- rm = adjustMomentum target tunables (q, r)
- qf = adjustPosition tunables (rm, q)
- rf = adjustMomentum target tunables (qf, rm)
-
--- | Adjust momentum according to a half-leapfrog step.
-adjustMomentum :: Target -> Tunables -> Particle -> Parameters
-adjustMomentum target tunables (q, r) = r .+ ((0.5 * e) .* g q) where
- e = stepSize tunables
- g = glTarget target
-
--- | Adjust position according to a half-leapfrog step.
-adjustPosition :: Tunables -> Particle -> Parameters
-adjustPosition tunables (r, q) = q .+ (e .* r) where
- e = stepSize tunables
-
--- | Scalar-vector multiplication.
-(.*) :: Double -> Parameters -> Parameters
-z .* xs = (* z) <$> xs
-
--- | Scalar-vector addition.
-(.+) :: Parameters -> Parameters -> Parameters
-xs .+ ys = V.zipWith (+) xs ys
-
--- | The next state of the Markov chain.
-nextState
- :: Target
- -> Particle
- -> Particle
- -> Double
- -> Parameters
-nextState target position momentum z
- | z < pAccept = snd position
- | otherwise = fst position
- where
- pAccept = acceptProb target position momentum
-
--- | A target augmented by momentum auxilliary variables.
-auxilliaryTarget :: Target -> Particle -> Double
-auxilliaryTarget target (t, r) = f t - 0.5 * innerProduct r r where
- f = lTarget target
-
--- | The acceptance probability of a move.
-acceptProb :: Target -> Particle -> Particle -> Double
-acceptProb target (q0, q1) (r0, r1) = exp . min 0 $
- auxilliaryTarget target (q1, r1) - auxilliaryTarget target (q0, r0)
-
--- | Simple inner product.
-innerProduct :: Parameters -> Parameters -> Double
-innerProduct xs ys = V.sum $ V.zipWith (*) xs ys
-
--- | Run a chain of HMC.
-hmc
- :: PrimMonad m
- => Target
- -> Tunables
- -> Options
- -> Gen (PrimState m)
- -> m [Parameters]
-hmc target tunables options g = do
- let h = hamiltonian target tunables g
- n = epochs options
- q = initial options
- evalStateT (replicateM n h) q
-
diff --git a/Numeric/Sampling/Hamiltonian.hs b/Numeric/Sampling/Hamiltonian.hs
@@ -0,0 +1,143 @@
+{-# OPTIONS_GHC -Wall #-}
+
+-- | Hamiltonian Monte Carlo. See, ex: Neal (2012)
+-- http://arxiv.org/pdf/1206.1901.pdf.
+
+module Numeric.Sampling.Hamiltonian (
+ Target(..)
+ , Tunables(..)
+ , Options(..)
+ , hamiltonian
+ , hmc
+ ) where
+
+import Control.Applicative
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.Trans.Class
+import Control.Monad.Trans.State.Strict
+import qualified Data.Foldable as Foldable
+import Data.Sampling.Types
+import Data.Vector (Vector)
+import qualified Data.Vector as V
+import System.Random.MWC
+import System.Random.MWC.Distributions
+
+data Tunables = Tunables {
+ stepSize :: !Double -- ^ step size for a given proposal
+ , leapfrogs :: !Int -- ^ number of leapfrog steps to take
+ } deriving Show
+
+data Options = Options {
+ epochs :: !Int -- ^ number of epochs to iterate the chain
+ , initial :: !Parameters -- ^ start location
+ } deriving Show
+
+type Particle = (Parameters, Parameters)
+
+type Transition m = StateT Parameters m Parameters
+
+-- | The Hamiltonian transition operator.
+hamiltonian
+ :: PrimMonad m
+ => Target
+ -> Tunables
+ -> Gen (PrimState m)
+ -> Transition m
+hamiltonian target tunables g = do
+ q0 <- get
+ r0 <- V.replicateM (V.length q0) (lift $ standard g)
+ zc <- lift $ uniform g
+ let (q, r) = leapfrogIntegrator target tunables (q0, r0)
+ next = nextState target (q0, q) (r0, r) zc
+ put next
+ return next
+{-# INLINE hamiltonian #-}
+
+-- | The leapfrog or Stormer-Verlet integrator.
+leapfrogIntegrator :: Target -> Tunables -> Particle -> Particle
+leapfrogIntegrator target tunables (q0, r0) = go q0 r0 l where
+ l = leapfrogs tunables
+ go q r 0 = (q, r)
+ go q r n =
+ let (q1, r1) = leapfrog target tunables (q, r)
+ in go q1 r1 (pred n)
+{-# INLINE leapfrogIntegrator #-}
+
+-- | A single iteration of the leapfrog integrator.
+leapfrog :: Target -> Tunables -> Particle -> Particle
+leapfrog target tunables (q, r) = (qf, rf) where
+ rm = adjustMomentum target tunables (q, r)
+ qf = adjustPosition tunables (rm, q)
+ rf = adjustMomentum target tunables (qf, rm)
+{-# INLINE leapfrog #-}
+
+-- | Adjust momentum according to a half-leapfrog step.
+adjustMomentum :: Target -> Tunables -> Particle -> Parameters
+adjustMomentum target tunables (q, r) = r .+ ((0.5 * e) .* g q) where
+ e = stepSize tunables
+ g = glTarget target
+{-# INLINE adjustMomentum #-}
+
+-- | Adjust position according to a half-leapfrog step.
+adjustPosition :: Tunables -> Particle -> Parameters
+adjustPosition tunables (r, q) = q .+ (e .* r) where
+ e = stepSize tunables
+{-# INLINE adjustPosition #-}
+
+-- | Scalar-vector multiplication.
+(.*) :: Double -> Parameters -> Parameters
+z .* xs = (* z) <$> xs
+{-# INLINE (.*) #-}
+
+-- | Scalar-vector addition.
+(.+) :: Parameters -> Parameters -> Parameters
+xs .+ ys = V.zipWith (+) xs ys
+{-# INLINE (.+) #-}
+
+-- | The next state of the Markov chain.
+nextState
+ :: Target
+ -> Particle
+ -> Particle
+ -> Double
+ -> Parameters
+nextState target position momentum z
+ | z < pAccept = snd position
+ | otherwise = fst position
+ where
+ pAccept = acceptProb target position momentum
+{-# INLINE nextState #-}
+
+-- | A target augmented by momentum auxilliary variables.
+auxilliaryTarget :: Target -> Particle -> Double
+auxilliaryTarget target (t, r) = f t - 0.5 * innerProduct r r where
+ f = lTarget target
+{-# INLINE auxilliaryTarget #-}
+
+-- | The acceptance probability of a move.
+acceptProb :: Target -> Particle -> Particle -> Double
+acceptProb target (q0, q1) (r0, r1) = exp . min 0 $
+ auxilliaryTarget target (q1, r1) - auxilliaryTarget target (q0, r0)
+{-# INLINE acceptProb #-}
+
+-- | Simple inner product.
+innerProduct :: Parameters -> Parameters -> Double
+innerProduct xs ys = V.sum $ V.zipWith (*) xs ys
+{-# INLINE innerProduct #-}
+
+-- | Run a chain of HMC.
+hmc
+ :: PrimMonad m
+ => Target
+ -> Tunables
+ -> Options
+ -> Gen (PrimState m)
+ -> m [Parameters]
+hmc target tunables options g = do
+ let h = hamiltonian target tunables g
+ n = epochs options
+ q = initial options
+ evalStateT (replicateM n h) q
+{-# INLINE hmc #-}
+
diff --git a/hasty-hamiltonian.cabal b/hasty-hamiltonian.cabal
@@ -1,5 +1,5 @@
name: hasty-hamiltonian
-version: 0.2.0.0
+version: 0.3.0.0
synopsis: Speedy traversal through parameter space.
homepage: http://jtobin.github.com/hasty-hamiltonian
license: BSD3
@@ -8,7 +8,7 @@ author: Jared Tobin
maintainer: jared@jtobin.ca
category: Numeric
build-type: Simple
-cabal-version: >=1.10
+cabal-version: >=1.18
Description:
Speedy gradient-based traversal through parameter space.
@@ -20,14 +20,15 @@ Source-repository head
library
default-language: Haskell2010
exposed-modules:
- Numeric.MCMC.Hamiltonian
+ Numeric.Sampling.Hamiltonian
build-depends:
- base >= 4.7.0.1 && < 5
- , mtl >= 2.2.1
- , mwc-random >= 0.13.2.0
- , primitive >= 0.5.3.0
- , transformers >= 0.4.1.0
- , vector >= 0.10.11.0
+ base
+ , mtl
+ , mcmc-types
+ , mwc-random
+ , primitive
+ , transformers
+ , vector
ghc-options:
-Wall
@@ -39,9 +40,9 @@ Test-Suite hasty-examples
ghc-options:
-threaded -rtsopts -Wall -O2
build-depends:
- ad >= 4.2.1
- , base >= 4.7.0.1 && < 5
- , hasty-hamiltonian >= 0.2.0.0
+ ad
+ , base
+ , hasty-hamiltonian
, mwc-random >= 0.13.2.0
- , vector >= 0.10.11.0
+ , vector
