commit 2a8520195be383c0239b8d5eff99f31635f7c502
parent c7d11968abdf85bcc01be17aaf2902ee19dd5e10
Author: Jared Tobin <jared@jtobin.ca>
Date: Tue, 6 Oct 2015 21:37:24 +1300
1.0.0 update.
Diffstat:
4 files changed, 54 insertions(+), 12 deletions(-)
diff --git a/Numeric/MCMC/Metropolis.hs b/Numeric/MCMC/Metropolis.hs
@@ -1,11 +1,30 @@
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE RecordWildCards #-}
+-- |
+-- Module: Numeric.MCMC.Metropolis
+-- Copyright: (c) 2015 Jared Tobin
+-- License: MIT
+--
+-- Maintainer: Jared Tobin <jared@jtobin.ca>
+-- Stability: unstable
+-- Portability: ghc
+--
+-- You have to wake up pretty early to beat the Metropolis algorithm.
+--
+-- This implementation uses spherical Gaussian proposals to implement a
+-- reliable and computationally inexpensive sampling routine. It can be used
+-- as a baseline from which to benchmarks other algorithms on a given problem.
+--
+-- The 'mcmc' function streams a trace to stdout to be processed elsewhere,
+-- while the `metropolis` transition can be used for more flexible purposes,
+-- such as working with samples in memory.
+
module Numeric.MCMC.Metropolis (
mcmc
, metropolis
- -- * re-export
+ -- * re-exported
, module Data.Sampling.Types
, MWC.create
, MWC.createSystemRandom
@@ -24,8 +43,8 @@ import qualified Pipes.Prelude as Pipes (mapM_, take)
import System.Random.MWC.Probability (Gen, Prob)
import qualified System.Random.MWC.Probability as MWC
--- | Propose a state transition according to a Gaussian proposal distribution
--- with the specified standard deviation.
+-- Propose a state transition according to a Gaussian proposal distribution
+-- with the specified standard deviation.
propose
:: (PrimMonad m, Traversable f)
=> Double
@@ -34,7 +53,7 @@ propose
propose radial = traverse perturb where
perturb m = MWC.normal m radial
--- | A Metropolis transition operator.
+-- | A generic Metropolis transition operator.
metropolis
:: (Traversable f, PrimMonad m)
=> Double
@@ -48,7 +67,7 @@ metropolis radial = do
accept <- lift (MWC.bernoulli acceptProb)
when accept (put (Chain chainTarget proposalScore proposal chainTunables))
--- | A Markov chain.
+-- A Markov chain driven by the Metropolis transition operator.
chain
:: (Traversable f, PrimMonad m)
=> Double
@@ -61,7 +80,13 @@ chain radial = loop where
yield next
loop next prng
--- | Trace 'n' iterations of a Markov chain.
+-- | Trace 'n' iterations of a Markov chain and stream them to stdout.
+--
+-- >>> let rosenbrock [x0, x1] = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+-- >>> withSystemRandom . asGenIO $ mcmc 3 1 [0, 0] rosenbrock
+-- 0.5000462419822702,0.5693944056267897
+-- 0.5000462419822702,0.5693944056267897
+-- -0.7525995304580824,1.2240725505283248
mcmc
:: (Traversable f, Show (f Double))
=> Int
@@ -79,7 +104,7 @@ mcmc n radial chainPosition target gen = runEffect $
chainTunables = Nothing
chainTarget = Target target Nothing
--- | Use a provided default value when the argument is NaN.
+-- Use a provided default value when the argument is NaN.
whenNaN :: RealFloat a => a -> a -> a
whenNaN val x
| isNaN x = val
diff --git a/README.md b/README.md
@@ -1,6 +1,19 @@
# mighty-metropolis [](http://travis-ci.org/jtobin/mighty-metropolis)
-The classic Metropolis sampling algorithm.
+The classic Metropolis algorithm. Wander around parameter space according to a
+simple spherical Gaussian distribution.
+
+Exports a `mcmc` function that prints a trace to stdout, as well as a
+`metropolis` transition operator that can be used more generally.
See the *test* folder for example usage.
+ import Numeric.MCMC.Metropolis
+
+ rosenbrock :: [Double] -> Double
+ rosenbrock [x0, x1] = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+
+ main :: IO ()
+ main = withSystemRandom . asGenIO $ mcmc 10000 1 [0, 0] rosenbrock
+
+
diff --git a/mighty-metropolis.cabal b/mighty-metropolis.cabal
@@ -1,5 +1,5 @@
name: mighty-metropolis
-version: 0.2.0.0
+version: 1.0.0
synopsis: The Metropolis algorithm.
homepage: http://github.com/jtobin/mighty-metropolis
license: MIT
@@ -37,5 +37,6 @@ Test-suite rosenbrock
-rtsopts -Wall
build-depends:
base
+ , containers
, mwc-probability
diff --git a/test/BNN.hs b/test/BNN.hs
@@ -2,10 +2,13 @@
module Main where
import Numeric.MCMC.Metropolis
+import Data.Sequence (Seq, fromList, index)
-bnn :: [Double] -> Double
-bnn [x0, x1] = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1)
+bnn :: Seq Double -> Double
+bnn xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1) where
+ x0 = index xs 0
+ x1 = index xs 1
main :: IO ()
-main = withSystemRandom . asGenIO $ mcmc 10000 1 [0, 0] bnn
+main = withSystemRandom . asGenIO $ mcmc 10000 1 (fromList [0, 0]) bnn
