declarative (empty) → 0.1.0.0

raw patch · 5 files changed

+299/−0 lines, 5 filesdep +basedep +declarativedep +hasty-hamiltoniansetup-changed

Dependencies added: base, declarative, hasty-hamiltonian, lens, mcmc-types, mighty-metropolis, mwc-probability, pipes, primitive, speedy-slice, transformers

Files

LICENSE +19/−0
Setup.hs +2/−0
declarative.cabal +78/−0
lib/Numeric/MCMC.hs +178/−0
test/Rosenbrock.hs +22/−0

+ LICENSE view

@@ -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.

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

+ declarative.cabal view

@@ -0,0 +1,78 @@+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.+  .+  Build composite Markov transition operators from existing ones for fun and+  profit.+  .+  A useful strategy is to hedge one's sampling risk by occasionally+  interleaving a computationally-expensive transition (such as a gradient-based+  algorithm like Hamiltonian Monte Carlo or NUTS) with cheap Metropolis+  transitions.+  .+  > transition = frequency [+  >     (9, metropolis 1.0)+  >   , (1, hamiltonian 0.05 20)+  >   ]+  .+  Alternatively: sample consecutively using the same algorithm, but over a+  range of different proposal distributions.+  .+  > transition = concatAllT [+  >     slice 0.5+  >   , slice 1.0+  >   , slice 2.0+  >   ]+  .+  Or just mix and match and see what happens!+  .+  > transition =+  >   sampleT+  >     (sampleT (metropolis 0.5) (slice 0.1))+  >     (sampleT (hamiltonian 0.01 20) (metropolis 2.0))+  .+  Check the test suite for example usage.++Source-repository head+  Type:     git+  Location: http://github.com/jtobin/declarative.git++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++Test-suite rosenbrock+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             Rosenbrock.hs+  default-language:    Haskell2010+  ghc-options:+    -rtsopts+  build-depends:+      base              < 5+    , mwc-probability   >= 1.0.1+    , declarative+

+ lib/Numeric/MCMC.hs view

@@ -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+

+ test/Rosenbrock.hs view

@@ -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+