flat-mcmc

Painless, efficient, general-purpose sampling from continuous distributions.
Log | Files | Refs | README | LICENSE

commit 9dc746eb8a341a58b45b677280d7e6be91149d94
parent f5e22006f136e705e999535dd1ecf2cd7a03c9a1
Author: Jared Tobin <jared@jtobin.ca>
Date:   Thu,  1 Dec 2016 10:25:08 +1300

Version bump (1.4.0).

Add a couple of convenience aliases for creating ensembles.

Diffstat:
ACHANGELOG | 7+++++++
MREADME.md | 22+++++++++++-----------
Mflat-mcmc.cabal | 35+++++++++++++++++++++--------------
Mlib/Numeric/MCMC/Flat.hs | 31+++++++++++++++++++++++--------
Msrc/Main.hs | 17++++++++---------
Mtest/BNN.hs | 17++++++++---------
Mtest/Rosenbrock.hs | 17++++++++---------
7 files changed, 86 insertions(+), 60 deletions(-)

diff --git a/CHANGELOG b/CHANGELOG @@ -0,0 +1,7 @@ +# Changelog + +- 1.4.0 (2016-12-01) + * Export the 'particle' and 'ensemble' helper aliases to make it a little + less burdensome to create an origin for a Markov chain. + + diff --git a/README.md b/README.md @@ -28,23 +28,23 @@ as a 'flat' transition operator that can be used more generally. ``` haskell import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (unsafeIndex) rosenbrock :: Particle -> Double rosenbrock xs = negate (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where - [x0, x1] = U.toList xs - -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] + x0 = U.unsafeIndex xs 0 + x1 = U.unsafeIndex xs 1 + +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 12500 ensemble rosenbrock +main = withSystemRandom . asGenIO $ mcmc 12500 origin rosenbrock ``` ![trace](http://jtobin.ca/flat-mcmc/img/Rosenbrock_AIE.png) diff --git a/flat-mcmc.cabal b/flat-mcmc.cabal @@ -1,5 +1,5 @@ name: flat-mcmc -version: 1.3.0 +version: 1.4.0 synopsis: Painless general-purpose sampling. homepage: https://github.com/jtobin/flat-mcmc license: MIT @@ -8,7 +8,7 @@ author: Jared Tobin maintainer: jared@jtobin.ca category: Math build-type: Simple -cabal-version: >=1.10 +cabal-version: >= 1.18 description: flat-mcmc is a Haskell library for painless, efficient, general-purpose sampling from continuous distributions. @@ -25,23 +25,23 @@ description: as a 'flat' transition operator that can be used more generally. . > import Numeric.MCMC.Flat - > import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) - > import qualified Data.Vector as V (fromList) + > import qualified Data.Vector.Unboxed as U (unsafeIndex) > > rosenbrock :: Particle -> Double > rosenbrock xs = negate (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where - > [x0, x1] = U.toList xs + > x0 = U.unsafeIndex xs 0 + > x1 = U.unsafeIndex xs 1 > - > ensemble :: Ensemble - > ensemble = V.fromList [ - > U.fromList [negate 1.0, negate 1.0] - > , U.fromList [negate 1.0, 1.0] - > , U.fromList [1.0, negate 1.0] - > , U.fromList [1.0, 1.0] + > origin :: Ensemble + > origin = ensemble [ + > particle [negate 1.0, negate 1.0] + > , particle [negate 1.0, 1.0] + > , particle [1.0, negate 1.0] + > , particle [1.0, 1.0] > ] > > main :: IO () - > main = withSystemRandom . asGenIO $ mcmc 12500 ensemble rosenbrock + > main = withSystemRandom . asGenIO $ mcmc 12500 origin rosenbrock Source-repository head Type: git @@ -49,9 +49,16 @@ Source-repository head library default-language: Haskell2010 - ghc-options: -Wall hs-source-dirs: lib - exposed-modules: Numeric.MCMC.Flat + ghc-options: + -Wall + + exposed-modules: + Numeric.MCMC.Flat + + other-modules: + Data.Vector.Extended + build-depends: base > 4 && < 6 , formatting >= 6 && < 7 diff --git a/lib/Numeric/MCMC/Flat.hs b/lib/Numeric/MCMC/Flat.hs @@ -36,6 +36,9 @@ module Numeric.MCMC.Flat ( , MWC.createSystemRandom , MWC.withSystemRandom , MWC.asGenIO + + , VE.ensemble + , VE.particle ) where import Control.Monad (replicateM) @@ -50,6 +53,7 @@ import qualified Data.Text as T import qualified Data.Text.IO as T (putStrLn) import Data.Vector (Vector) import qualified Data.Vector as V +import qualified Data.Vector.Extended as VE (ensemble, particle) import qualified Data.Vector.Unboxed as U import Formatting ((%)) import qualified Formatting as F @@ -83,8 +87,20 @@ renderEnsemble = glue a b = a <> "\n" <> renderParticle b {-# INLINE renderEnsemble #-} +-- | A particle is an n-dimensional point in Euclidean space. +-- +-- You can create a particle by using the 'particle' helper function, or just +-- use Data.Vector.Unboxed.fromList. type Particle = U.Vector Double +-- | An ensemble is a collection of particles. +-- +-- The Markov chain we're interested in will run over the space of ensembles, +-- so you'll want to build an ensemble out of a reasonable number of +-- particles to kick off the chain. +-- +-- You can create an ensemble by using the 'ensemble' helper function, or just +-- use Data.Vector.fromList. type Ensemble = Vector Particle symmetric :: PrimMonad m => Prob m Double @@ -130,7 +146,6 @@ execute target e0 e1 n = do w0 k = e0 `V.unsafeIndex` pred k w1 k ks = e1 `V.unsafeIndex` pred (ks `U.unsafeIndex` pred k) - worker (k, z, zc) = move target (w0 k) (w1 k js) z zc !result = runPar $ parMapChunk granularity worker (zip3 [1..n] zs zcs) @@ -169,20 +184,20 @@ chain = loop where -- you'll need to provide an ensemble of particles for the start location. -- -- >>> import Numeric.MCMC.Flat --- >>> import Data.Vector (Vector, toList, fromList) +-- >>> import Data.Vector.Unboxed (toList) -- >>> :{ -- >>> let rosenbrock xs = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) -- where [x0, x1] = toList xs -- >>> :} -- >>> :{ --- >>> let ensemble = fromList [ --- >>> fromList [negate 1.0, negate 1.0] --- >>> , fromList [negate 1.0, 1.0] --- >>> , fromList [1.0, negate 1.0] --- >>> , fromList [1.0, 1.0] +-- >>> let origin = ensemble [ +-- >>> particle [negate 1.0, negate 1.0] +-- >>> , particle [negate 1.0, 1.0] +-- >>> , particle [1.0, negate 1.0] +-- >>> , particle [1.0, 1.0] -- >>> ] -- >>> :} --- >>> withSystemRandom . asGenIO $ mcmc 2 ensemble rosenbrock +-- >>> withSystemRandom . asGenIO $ mcmc 2 origin rosenbrock -- -1.0,-1.0 -- -1.0,1.0 -- 1.0,-1.0 diff --git a/src/Main.hs b/src/Main.hs @@ -3,21 +3,20 @@ module Main where import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (toList) bnn :: Particle -> Double bnn xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1) where [x0, x1] = U.toList xs -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 10000 ensemble bnn +main = withSystemRandom . asGenIO $ mcmc 10000 origin bnn diff --git a/test/BNN.hs b/test/BNN.hs @@ -3,21 +3,20 @@ module Main where import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (toList) bnn :: Particle -> Double bnn xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1) where [x0, x1] = U.toList xs -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 100 ensemble bnn +main = withSystemRandom . asGenIO $ mcmc 100 origin bnn diff --git a/test/Rosenbrock.hs b/test/Rosenbrock.hs @@ -3,21 +3,20 @@ module Main where import Numeric.MCMC.Flat -import qualified Data.Vector.Unboxed as U (Vector, toList, fromList) -import qualified Data.Vector as V (fromList) +import qualified Data.Vector.Unboxed as U (toList) rosenbrock :: Particle -> Double rosenbrock xs = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2) where [x0, x1] = U.toList xs -ensemble :: Ensemble -ensemble = V.fromList [ - U.fromList [negate 1.0, negate 1.0] - , U.fromList [negate 1.0, 1.0] - , U.fromList [1.0, negate 1.0] - , U.fromList [1.0, 1.0] +origin :: Ensemble +origin = ensemble [ + particle [negate 1.0, negate 1.0] + , particle [negate 1.0, 1.0] + , particle [1.0, negate 1.0] + , particle [1.0, 1.0] ] main :: IO () -main = withSystemRandom . asGenIO $ mcmc 100 ensemble rosenbrock +main = withSystemRandom . asGenIO $ mcmc 100 origin rosenbrock