commit c5a0a5fa541d4a29bd9b3c77373349e9dc56ac8d
Author: Jared Tobin <jared@jtobin.ca>
Date: Tue, 6 Nov 2012 20:55:26 +1300
Initial commit.
Diffstat:
6 files changed, 213 insertions(+), 0 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,14 @@
+.git
+*swp
+.DS_Store
+dist
+Setup.hs
+*.hi
+*.o
+demos
+data
+_cache/
+_site/
+build
+cabal-dev
+site
diff --git a/Examples/code/Rosenbrock_MH.hs b/Examples/code/Rosenbrock_MH.hs
@@ -0,0 +1,39 @@
+import Numeric.MCMC.Metropolis
+import System.Random.MWC
+import System.Environment
+import System.Exit
+import System.IO
+import Control.Monad
+
+target :: RealFloat a => [a] -> a
+target [x0, x1] = (-1)*(5*(x1 - x0^2)^2 + 0.05*(1 - x0)^2)
+
+main = do
+ args <- getArgs
+ when (args == []) $ do
+ putStrLn "(mighty-metropolis) Rosenbrock density "
+ putStrLn "Usage: ./Rosenbrock_MH <numSteps> <stepSize> <inits> "
+ putStrLn " "
+ putStrLn "numSteps : Number of Markov chain iterations to run."
+ 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 !! 2)) :: IO [Double]
+
+ let nepochs = read (head args) :: Int
+ eps = read (args !! 1) :: Double
+ params = Options target eps
+ config = MarkovChain inits 0
+
+ g <- create
+ results <- runChain params nepochs 1 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/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/Metropolis.hs b/Numeric/MCMC/Metropolis.hs
@@ -0,0 +1,95 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE BangPatterns #-}
+
+module Numeric.MCMC.Metropolis (
+ MarkovChain(..), Options(..)
+ , runChain
+ ) where
+
+import Control.Monad.Trans
+import Control.Monad.Reader
+import Control.Monad.Primitive
+import Control.Arrow
+import System.Random.MWC
+import System.Random.MWC.Distributions
+import Data.List
+import Statistics.Distribution
+import Statistics.Distribution.Normal hiding (standard)
+import GHC.Float
+
+-- | 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) and
+-- a step size tuning parameter.
+data Options = Options { _target :: [Double] -> 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 (map double2Float (theta config))
+
+-- | Density function for an isotropic Gaussian. The (identity) covariance
+-- matrix is multiplied by the scalar 'sig'.
+isoGauss :: [Double] -> [Double] -> Double -> Double
+isoGauss xs mu sig = foldl1' (*) (zipWith density nds xs)
+ where nds = map (`normalDistr` sig) mu
+{-# INLINE isoGauss #-}
+
+-- | Perturb the state, creating a new proposal.
+perturb :: PrimMonad m
+ => [Double] -- Current state
+ -> Gen (PrimState m) -- MWC PRNG
+ -> ViewsOptions m [Double] -- Resulting perturbation.
+perturb t g = do
+ Options _ e <- ask
+ mapM (\m -> lift $ normal m e g) t
+{-# INLINE perturb #-}
+
+-- | Perform a Metropolis accept/reject step.
+metropolisStep :: PrimMonad m
+ => MarkovChain -- Current state
+ -> Gen (PrimState m) -- MWC PRNG
+ -> ViewsOptions m MarkovChain -- New state
+metropolisStep state g = do
+ Options target e <- ask
+ let (t0, nacc) = (theta &&& accepts) state
+ zc <- lift $ uniformR (0, 1) g
+ proposal <- perturb t0 g
+ let mc = if zc < acceptProb
+ then (proposal, 1)
+ else (t0, 0)
+
+ acceptProb = if isNaN val then 0 else val where val = arRatio
+
+ arRatio = exp . min 0 $
+ target proposal + log (isoGauss t0 proposal e)
+ - target t0 - log (isoGauss proposal t0 e)
+
+ return $! MarkovChain (fst mc) (nacc + snd mc)
+{-# INLINE metropolisStep #-}
+
+-- | Diffuse through states.
+runChain :: Options -- Options of the Markov chain.
+ -> Int -- Number of epochs to iterate the chain.
+ -> Int -- Print every nth iteration
+ -> MarkovChain -- Initial state of the Markov chain.
+ -> Gen RealWorld -- MWC PRNG
+ -> IO MarkovChain -- End state of the Markov chain, wrapped in IO.
+runChain = go
+ where go o n t !c g | n == 0 = return c
+ | n `rem` t /= 0 = do
+ r <- runReaderT (metropolisStep c g) o
+ go o (n - 1) t r g
+ | otherwise = do
+ r <- runReaderT (metropolisStep c g) o
+ print r
+ go o (n - 1) t r g
+{-# INLINE runChain #-}
+
+
diff --git a/README.md b/README.md
@@ -0,0 +1,6 @@
+# mighty-metropolis
+
+The classic Metropolis-Hastings sampling algorithm.
+
+See the *Examples* folder for example usage.
+
diff --git a/mighty-metropolis.cabal b/mighty-metropolis.cabal
@@ -0,0 +1,29 @@
+-- Initial mighty-metropolis.cabal generated by cabal init. For further
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name: mighty-metropolis
+version: 0.1.0.0
+synopsis: The classic Metropolis-Hastings sampling algorithm.
+-- description:
+homepage: http://github.com/jtobin/mighty-metropolis
+license: BSD3
+license-file: LICENSE
+author: Jared Tobin
+maintainer: jared@jtobin.ca
+-- copyright:
+category: Numeric
+build-type: Simple
+cabal-version: >=1.8
+Description:
+
+ Sampling via Gaussian perturbations.
+
+Source-repository head
+ Type: git
+ Location: http://github.com/jtobin/mighty-metropolis.git
+
+library
+ exposed-modules: Numeric.MCMC.Metropolis
+ -- other-modules:
+ build-depends: base ==4.*, mtl ==2.1.*, primitive ==0.4.*, mwc-random ==0.12.*, statistics ==0.10.*
+ ghc-options: -Wall
