commit cc012bc576452288968201ab48fe61f0ab204682
Author: Jared Tobin <jared@jtobin.ca>
Date: Sat, 10 Oct 2015 00:02:46 +1300
Initial commit.
Diffstat:
9 files changed, 323 insertions(+), 0 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,6 @@
+*swp
+.stack-work
+debug
+*o
+*hi
+
diff --git a/.travis.yml b/.travis.yml
@@ -0,0 +1,2 @@
+language: haskell
+sudo: false
diff --git a/LICENSE b/LICENSE
@@ -0,0 +1,19 @@
+Copyright (c) 2012-2015 Jared Tobin
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff --git a/README.md b/README.md
@@ -0,0 +1,53 @@
+# declarative
+
+[](http://travis-ci.org/jtobin/declarative)
+[](http://hackage.haskell.org/package/declarative)
+
+DIY Markov Chains.
+
+This package presents a simple combinator language for Markov transition
+operators that are useful in MCMC.
+
+Any transition operators sharing the same stationary distribution and obeying
+the Markov and reversibility properties can be combined in a couple of ways,
+such that the resulting operator preserves the stationary distribution and
+desirable properties amenable for MCMC.
+
+We can deterministically concatenate operators end-to-end, or sample from
+a collection of them according to some probability distribution. See
+[Geyer, 2005](www.stat.umn.edu/geyer/f05/8931/n1998.pdf) for details.
+
+The result is a simple grammar for building composite, property-preserving
+transition operators from existing ones:
+
+ transition ::= primitive <transition>
+ | concatT transition transition
+ | sampleT transition transition
+
+In addition to the above, this module provides a number of combinators for
+building composite transition operators. It re-exports a number of
+production-quality transition operators from the *mighty-metropolis*,
+*speedy-slice*, and *hasty-hamiltonian* libraries.
+
+Markov chains can then be run over arbitrary `Target`s using whatever
+transition operator is desired.
+
+ import Numeric.MCMC
+ import Data.Sampling.Types
+
+ target :: [Double] -> Double
+ target [x0, x1] = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+
+ rosenbrock :: Target [Double]
+ rosenbrock = Target target Nothing
+
+ transition :: Transition IO (Chain [Double] b)
+ transition =
+ concatT
+ (sampleT (metropolis 0.5) (metropolis 1.0))
+ (sampleT (slice 2.0) (slice 3.0))
+
+ main :: IO ()
+ main = withSystemRandom . asGenIO $ mcmc 10000 [0, 0] rosenbrock transition
+
+
diff --git a/Setup.hs b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/declarative.cabal b/declarative.cabal
@@ -0,0 +1,31 @@
+name: declarative
+version: 0.1.0.0
+synopsis: DIY Markov Chains.
+homepage: http://github.com/jtobin/declarative
+license: MIT
+license-file: LICENSE
+author: Jared Tobin
+maintainer: jared@jtobin.ca
+category: Math
+build-type: Simple
+cabal-version: >=1.10
+description:
+ DIY Markov Chains.
+
+library
+ default-language: Haskell2010
+ hs-source-dirs: lib
+ exposed-modules:
+ Numeric.MCMC
+ build-depends:
+ base < 5
+ , mcmc-types >= 1.0.1
+ , mwc-probability >= 1.0.1
+ , mighty-metropolis >= 1.0.1
+ , lens >= 4 && < 5
+ , primitive
+ , pipes >= 4 && < 5
+ , hasty-hamiltonian >= 1.1.1
+ , speedy-slice >= 0.1.2
+ , transformers
+
diff --git a/lib/Numeric/MCMC.hs b/lib/Numeric/MCMC.hs
@@ -0,0 +1,178 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module: Numeric.MCMC
+-- Copyright: (c) 2015 Jared Tobin
+-- License: MIT
+--
+-- Maintainer: Jared Tobin <jared@jtobin.ca>
+-- Stability: unstable
+-- Portability: ghc
+--
+-- This module presents a simple combinator language for Markov transition
+-- operators that are useful in MCMC.
+--
+-- Any transition operators sharing the same stationary distribution and
+-- obeying the Markov and reversibility properties can be combined in a couple
+-- of ways, such that the resulting operator preserves the stationary
+-- distribution and desirable properties amenable for MCMC.
+--
+-- We can deterministically concatenate operators end-to-end, or sample from
+-- a collection of them according to some probability distribution. See
+-- <www.stat.umn.edu/geyer/f05/8931/n1998.pdf Geyer, 2005> for details.
+--
+-- The result is a simple grammar for building composite, property-preserving
+-- transition operators from existing ones:
+--
+-- @
+-- transition ::= primitive <transition>
+-- | concatT transition transition
+-- | sampleT transition transition
+-- @
+--
+-- In addition to the above, this module provides a number of combinators for
+-- building composite transition operators. It re-exports a number of
+-- production-quality transition operators from the /mighty-metropolis/,
+-- /speedy-slice/, and /hasty-hamiltonian/ libraries.
+--
+-- Markov chains can then be run over arbitrary 'Target's using whatever
+-- transition operator is desired.
+--
+-- > import Numeric.MCMC
+-- > import Data.Sampling.Types
+-- >
+-- > target :: [Double] -> Double
+-- > target [x0, x1] = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+-- >
+-- > rosenbrock :: Target [Double]
+-- > rosenbrock = Target target Nothing
+-- >
+-- > transition :: Transition IO (Chain [Double] b)
+-- > transition =
+-- > concatT
+-- > (sampleT (metropolis 0.5) (metropolis 1.0))
+-- > (sampleT (slice 2.0) (slice 3.0))
+-- >
+-- > main :: IO ()
+-- > main = withSystemRandom . asGenIO $ mcmc 10000 [0, 0] rosenbrock transition
+--
+-- See the attached test suite for other examples.
+
+module Numeric.MCMC (
+ concatT
+ , concatAllT
+ , sampleT
+ , sampleAllT
+ , bernoulliT
+ , frequency
+ , mcmc
+
+ -- * Re-exported
+ , module Data.Sampling.Types
+
+ , metropolis
+ , hamiltonian
+ , slice
+
+ , MWC.create
+ , MWC.createSystemRandom
+ , MWC.withSystemRandom
+ , MWC.asGenIO
+ ) where
+
+import Control.Lens hiding (index)
+import Control.Monad.Primitive (PrimMonad, PrimState, RealWorld)
+import Control.Monad.Trans.State.Strict (execStateT)
+import Data.Sampling.Types
+import Numeric.MCMC.Metropolis hiding (mcmc)
+import Numeric.MCMC.Hamiltonian hiding (mcmc)
+import Numeric.MCMC.Slice hiding (mcmc)
+import Pipes hiding (next)
+import qualified Pipes.Prelude as Pipes
+import System.Random.MWC.Probability (Gen, Variate)
+import qualified System.Random.MWC.Probability as MWC
+
+-- | Deterministically concat transition operators together.
+concatT :: Monad m => Transition m a -> Transition m a -> Transition m a
+concatT = (>>)
+
+-- | Deterministically concat a list of transition operators together.
+concatAllT :: Monad m => [Transition m a] -> Transition m a
+concatAllT = foldl1 (>>)
+
+-- | Probabilistically concat transition operators together.
+sampleT :: PrimMonad m => Transition m a -> Transition m a -> Transition m a
+sampleT = bernoulliT 0.5
+
+-- | Probabilistically concat transition operators together using a Bernoulli
+-- distribution with the supplied success probability.
+--
+-- This is just a generalization of sampleT.
+bernoulliT
+ :: PrimMonad m
+ => Double
+ -> Transition m a
+ -> Transition m a
+ -> Transition m a
+bernoulliT p t0 t1 = do
+ heads <- lift (MWC.bernoulli p)
+ if heads then t0 else t1
+
+-- | Probabilistically concat transition operators together via a uniform
+-- distribution.
+sampleAllT :: PrimMonad m => [Transition m a] -> Transition m a
+sampleAllT ts = do
+ j <- lift (MWC.uniformR (0, length ts - 1))
+ ts !! j
+
+-- | Probabilistically concat transition operators together using the supplied
+-- frequency distribution.
+--
+-- This function is more-or-less an exact copy of 'QuickCheck.frequency',
+-- except here applied to transition operators.
+frequency :: PrimMonad m => [(Int, Transition m a)] -> Transition m a
+frequency xs = lift (MWC.uniformR (1, tot)) >>= (`pick` xs) where
+ tot = sum . map fst $ xs
+ pick n ((k, v):vs)
+ | n <= k = v
+ | otherwise = pick (n - k) vs
+ pick _ _ = error "frequency: no distribution specified"
+
+-- | Trace 'n' iterations of a Markov chain and stream them to stdout.
+--
+-- >>> withSystemRandom . asGenIO $ mcmc 3 [0, 0] (metropolis 0.5) rosenbrock
+-- -0.48939312153007863,0.13290702689491818
+-- 1.4541485365128892e-2,-0.4859905564050404
+-- 0.22487398491619448,-0.29769783186855125
+mcmc
+ :: (Show (t a), FoldableWithIndex (Index (t a)) t, Ixed (t a),
+ Num (IxValue (t a)), Variate (IxValue (t a)))
+ => Int
+ -> t a
+ -> Transition IO (Chain (t a) b)
+ -> Target (t a)
+ -> Gen RealWorld
+ -> IO ()
+mcmc n chainPosition transition chainTarget gen = runEffect $
+ chain transition Chain {..} gen
+ >-> Pipes.take n
+ >-> Pipes.mapM_ print
+ where
+ chainScore = lTarget chainTarget chainPosition
+ chainTunables = Nothing
+
+-- A Markov chain driven by an arbitrary transition operator.
+chain
+ :: PrimMonad m
+ => Transition m b
+ -> b
+ -> Gen (PrimState m)
+ -> Producer b m a
+chain transition = loop where
+ loop state prng = do
+ next <- lift (MWC.sample (execStateT transition state) prng)
+ yield next
+ loop next prng
+
diff --git a/stack.yaml b/stack.yaml
@@ -0,0 +1,10 @@
+resolver: lts-3.8
+packages:
+- '.'
+extra-deps:
+ - mighty-metropolis-1.0.2
+ - hasty-hamiltonian-1.1.2
+ - speedy-slice-0.1.3
+ - mcmc-types-1.0.1
+ - mwc-probability-1.0.2
+flags: {}
diff --git a/test/Rosenbrock.hs b/test/Rosenbrock.hs
@@ -0,0 +1,22 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Main where
+
+import Numeric.MCMC
+import Data.Sampling.Types
+
+target :: [Double] -> Double
+target [x0, x1] = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+
+rosenbrock :: Target [Double]
+rosenbrock = Target target Nothing
+
+transition :: Transition IO (Chain [Double] b)
+transition =
+ concatT
+ (sampleT (metropolis 0.5) (metropolis 1.0))
+ (sampleT (slice 2.0) (slice 3.0))
+
+main :: IO ()
+main = withSystemRandom . asGenIO $ mcmc 10000 [0, 0] rosenbrock transition
+
