commit c6d45d0031e8044eb2d20fdb6dbd95d04a9de4ec
Author: Jared Tobin <jared@jtobin.ca>
Date: Sun, 4 Nov 2012 22:32:10 +1300
Initial commit.
Diffstat:
7 files changed, 190 insertions(+), 0 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,10 @@
+working
+*swp
+.git
+.DS_Store
+dist
+Setup.hs
+*trace*
+demos
+output
+
diff --git a/Examples/code/Rosenbrock_HMC.hs b/Examples/code/Rosenbrock_HMC.hs
@@ -0,0 +1,45 @@
+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> <nDisc> <stepSize> <inits>"
+ putStrLn " "
+ putStrLn "numSteps : Number of Markov chain iterations to run. "
+ 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 !! 3)) :: IO [Double]
+
+ let nepochs = read (head args) :: Int
+ nDisc = read (args !! 1) :: Double
+ eps = read (args !! 2) :: Double
+ params = Options target gTarget nDisc eps
+ config = MarkovChain inits 0
+
+ g <- create
+ results <- runChain params nepochs 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/data/input/inits-2d.dat b/Examples/data/input/inits-2d.dat
@@ -0,0 +1 @@
+5.0 40.0
diff --git a/LICENSE b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2012, Jared Tobin
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+
+ * Neither the name of Jared Tobin nor the names of other
+ contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Numeric/MCMC/Hamiltonian.hs b/Numeric/MCMC/Hamiltonian.hs
@@ -0,0 +1,81 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE BangPatterns #-}
+
+module Numeric.MCMC.Hamiltonian where
+
+import Control.Monad.Reader
+import Control.Monad.Primitive
+import System.Random.MWC
+import System.Random.MWC.Distributions
+
+-- | State of the Markov chain. Current parameter values are held in 'theta',
+-- while accepts counts the number of proposals accepted.
+data MarkovChain = MarkovChain { theta :: [Double]
+ , accepts :: {-# UNPACK #-} !Int }
+
+-- | Options for the chain. The target (expected to be a log density), its
+-- gradient, and a step size tuning parameter.
+data Options = Options { _target :: [Double] -> Double
+ , _gTarget :: [Double] -> [Double]
+ , _nLeaps :: {-# UNPACK #-} !Double
+ , _eps :: {-# UNPACK #-} !Double }
+
+-- | A result with this type has a view of the chain options.
+type ViewsOptions = ReaderT Options
+
+-- | Display the current state.
+instance Show MarkovChain where
+ show config = filter (`notElem` "[]") $ show (theta config)
+
+-- | The 'leapfrog' or Stormer-Verlet discretizer.
+leapfrog :: Monad m
+ => [Double]
+ -> [Double]
+ -> ViewsOptions m ([Double], [Double])
+leapfrog t0 r0 = do
+ Options _ gTarget ndisc e <- ask
+ let go !t !r n
+ | n == 0 = (t, r)
+ | otherwise = let rm = zipWith (+) r (map (* (0.5*e)) (gTarget t))
+ tt = zipWith (+) t (map (*e) rm)
+ rt = zipWith (+) rm (map (* (0.5*e)) (gTarget tt))
+ in go tt rt (n - 1)
+ return $! go t0 r0 ndisc
+
+-- | Perform a Metropolis accept/reject step.
+metropolisStep :: PrimMonad m
+ => MarkovChain
+ -> Gen (PrimState m)
+ -> ViewsOptions m MarkovChain
+metropolisStep state g = do
+ Options target _ _ _ <- ask
+ let t0 = theta state
+ nacc = accepts state
+ r0 <- replicateM (length t0) (lift $ standard g)
+ zc <- lift $ uniformR (0 :: Double, 1 :: Double) g
+ (t, r) <- leapfrog t0 r0
+
+ let mc = if zc < min 1 (exp arRatio)
+ then (t, 1)
+ else (t0, 0)
+
+ xs <.> ys = sum $ zipWith (*) xs ys
+ arRatio = target t + 0.5*(r0 <.> r0)
+ - target t0 - 0.5*(r <.> r)
+
+ return $! MarkovChain (fst mc) (nacc + snd mc)
+
+-- | Diffuse through states.
+runChain :: Options -- Options of the Markov chain.
+ -> Int -- Number of epochs to iterate the chain.
+ -> MarkovChain -- Initial state of the Markov chain.
+ -> Gen RealWorld -- MWC PRNG
+ -> IO MarkovChain -- End state of the Markov chain, wrapped in IO.
+runChain params nepochs initConfig g
+ | nepochs == 0 = return initConfig
+ | otherwise = do
+ result <- runReaderT (metropolisStep initConfig g) params
+ print result
+ runChain params (nepochs - 1) result g
+
+
diff --git a/README.md b/README.md
@@ -0,0 +1,2 @@
+Speedy, gradient-based traversal through parameter space.
+
diff --git a/hasty-hamiltonian.cabal b/hasty-hamiltonian.cabal
@@ -0,0 +1,21 @@
+-- Initial hasty-hamiltonian.cabal generated by cabal init. For further
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name: hasty-hamiltonian
+version: 0.1.0.0
+synopsis: Speedy traversal through parameter space.
+-- description:
+homepage: http://jtobin.github.com/hasty-hamiltonian
+license: BSD3
+license-file: LICENSE
+author: Jared Tobin
+maintainer: jared@jtobin.ca
+-- copyright:
+category: Numeric
+build-type: Simple
+cabal-version: >=1.8
+
+library
+ exposed-modules: Numeric.MCMC.Hamiltonian
+ -- other-modules:
+ build-depends: base ==4.5.*, mtl ==2.1.*, primitive ==0.4.*, mwc-random ==0.12.*
