mighty-metropolis (empty) → 1.0.0

raw patch · 6 files changed

+233/−0 lines, 6 filesdep +basedep +containersdep +ghc-primsetup-changed

Dependencies added: base, containers, ghc-prim, mcmc-types, mighty-metropolis, mwc-probability, pipes, primitive, transformers

Files

LICENSE +19/−0
Numeric/MCMC/Metropolis.hs +113/−0
Setup.hs +2/−0
mighty-metropolis.cabal +72/−0
test/BNN.hs +15/−0
test/Rosenbrock.hs +12/−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.

+ Numeric/MCMC/Metropolis.hs view

@@ -0,0 +1,113 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE RecordWildCards #-}++-- |+-- Module: Numeric.MCMC.Metropolis+-- Copyright: (c) 2015 Jared Tobin+-- License: MIT+--+-- Maintainer: Jared Tobin <jared@jtobin.ca>+-- Stability: unstable+-- Portability: ghc+--+-- 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 benchmark other algorithms for 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-exported+  , module Data.Sampling.Types+  , MWC.create+  , MWC.createSystemRandom+  , MWC.withSystemRandom+  , MWC.asGenIO+  ) where++import Control.Monad (when)+import Control.Monad.Primitive (PrimMonad, PrimState)+import Control.Monad.Trans.Class (lift)+import Control.Monad.Trans.State.Strict (execStateT, get, put)+import Data.Sampling.Types (Target(..), Chain(..), Transition)+#if __GLASGOW_HASKELL__ < 710+import Data.Traversable (Traversable, traverse)+#endif+import GHC.Prim (RealWorld)+import Pipes (Producer, yield, (>->), runEffect)+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+  :: (PrimMonad m, Traversable f)+  => Double+  -> f Double+  -> Prob m (f Double)+propose radial = traverse perturb where+  perturb m = MWC.normal m radial++-- | A generic Metropolis transition operator.+metropolis+  :: (Traversable f, PrimMonad m)+  => Double+  -> Transition m (Chain (f Double) b)+metropolis radial = do+  Chain {..} <- get+  proposal <- lift (propose radial chainPosition)+  let proposalScore = lTarget chainTarget proposal+      acceptProb    = whenNaN 0 (exp (min 0 (proposalScore - chainScore)))++  accept <- lift (MWC.bernoulli acceptProb)+  when accept (put (Chain chainTarget proposalScore proposal chainTunables))++-- A Markov chain driven by the Metropolis transition operator.+chain+  :: (Traversable f, PrimMonad m)+  => Double+  -> Chain (f Double) b+  -> Gen (PrimState m)+  -> Producer (Chain (f Double) b) m ()+chain radial = loop where+  loop state prng = do+    next <- lift (MWC.sample (execStateT (metropolis radial) state) prng)+    yield next+    loop next prng++-- | 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+  -> Double+  -> f Double+  -> (f Double -> Double)+  -> Gen RealWorld+  -> IO ()+mcmc n radial chainPosition target gen = runEffect $+        chain radial Chain {..} gen+    >-> Pipes.take n+    >-> Pipes.mapM_ print+  where+    chainScore    = lTarget chainTarget chainPosition+    chainTunables = Nothing+    chainTarget   = Target target Nothing++-- Use a provided default value when the argument is NaN.+whenNaN :: RealFloat a => a -> a -> a+whenNaN val x+  | isNaN x   = val+  | otherwise = x

+ Setup.hs view

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

+ mighty-metropolis.cabal view

@@ -0,0 +1,72 @@+name:                mighty-metropolis+version:             1.0.0+synopsis:            The Metropolis algorithm.+homepage:            http://github.com/jtobin/mighty-metropolis+license:             MIT+license-file:        LICENSE+author:              Jared Tobin+maintainer:          jared@jtobin.ca+category:            Numeric+build-type:          Simple+cabal-version:       >= 1.10+description:+  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.+  .+  > 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++Source-repository head+  Type:     git+  Location: http://github.com/jtobin/mighty-metropolis.git++library+  default-language:    Haskell2010+  ghc-options:+      -Wall+  exposed-modules:+      Numeric.MCMC.Metropolis+  build-depends:+      base             < 5+    , ghc-prim+    , pipes+    , primitive+    , mcmc-types       >= 1.0.1+    , mwc-probability  >= 1.0.1+    , 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+    , mighty-metropolis+    , mwc-probability   >= 1.0.1++Test-suite bnn+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             BNN.hs+  default-language:    Haskell2010+  ghc-options:+    -rtsopts+  build-depends:+      base              < 5+    , containers+    , mighty-metropolis+    , mwc-probability   >= 1.0.1+

+ test/BNN.hs view

@@ -0,0 +1,15 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Main where++import Numeric.MCMC.Metropolis+import Data.Sequence (Seq, fromList, index)++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 (fromList [0, 0]) bnn+

+ test/Rosenbrock.hs view

@@ -0,0 +1,12 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Main where++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+