commit cd7b9e5f40201d3b3e066df83d87d8a9ffe2b914
parent 78372fbccb60397a57bc0b33dd66f03f5fc64f9f
Author: Jared Tobin <jared@jtobin.ca>
Date: Thu, 8 Oct 2015 21:04:34 +1300
1.1.0 update.
Diffstat:
10 files changed, 372 insertions(+), 199 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -19,3 +19,8 @@ tmp
*hp
*ps
*aux
+.stack-work
+debug
+test/Beale
+test/Booth
+test/Rosenbrock
diff --git a/LICENSE b/LICENSE
@@ -1,30 +1,19 @@
-Copyright (c) 2012, Jared Tobin
+Copyright (c) 2012-2015 Jared Tobin
-All rights reserved.
+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:
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions are met:
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
- * Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
-
- * Redistributions in binary form must reproduce the above
- copyright notice, this list of conditions and the following
- disclaimer in the documentation and/or other materials provided
- with the distribution.
-
- * Neither the name of Jared Tobin nor the names of other
- contributors may be used to endorse or promote products derived
- from this software without specific prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
-OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
-DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
-THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+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/Numeric/MCMC/Hamiltonian.hs b/Numeric/MCMC/Hamiltonian.hs
@@ -0,0 +1,218 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module: Numeric.MCMC.Hamiltonian
+-- Copyright: (c) 2015 Jared Tobin
+-- License: MIT
+--
+-- Maintainer: Jared Tobin <jared@jtobin.ca>
+-- Stability: unstable
+-- Portability: ghc
+--
+-- This implementation performs Hamiltonian Monte Carlo using an identity mass
+-- matrix.
+--
+-- The 'mcmc' function streams a trace to stdout to be processed elsewhere,
+-- while the `slice` transition can be used for more flexible purposes, such as
+-- working with samples in memory.
+--
+-- See <http://arxiv.org/pdf/1206.1901.pdf Neal, 2012> for the definitive
+-- reference of the algorithm.
+
+module Numeric.MCMC.Hamiltonian (
+ mcmc
+ , hamiltonian
+
+ -- * Re-exported
+ , Target(..)
+ , MWC.create
+ , MWC.createSystemRandom
+ , MWC.withSystemRandom
+ , MWC.asGenIO
+ ) where
+
+import Control.Lens hiding (index)
+import Control.Monad.Trans.State.Strict hiding (state)
+import Control.Monad.Primitive (PrimState, PrimMonad, RealWorld)
+import qualified Data.Foldable as Foldable (sum)
+import Data.Maybe (fromMaybe)
+import Data.Sampling.Types
+import Data.Traversable (for)
+import Pipes hiding (for, next)
+import qualified Pipes.Prelude as Pipes
+import System.Random.MWC.Probability (Prob, Gen)
+import qualified System.Random.MWC.Probability as MWC
+
+-- | Trace 'n' iterations of a Markov chain and stream them to stdout.
+--
+-- >>> withSystemRandom . asGenIO $ mcmc 3 1 [0, 0] target
+mcmc
+ :: (Num (IxValue (t Double)), Show (t Double), Traversable t
+ , FunctorWithIndex (Index (t Double)) t, Ixed (t Double)
+ , IxValue (t Double) ~ Double)
+ => Int
+ -> Double
+ -> Int
+ -> t Double
+ -> Target (t Double)
+ -> Gen RealWorld
+ -> IO ()
+mcmc n step leaps chainPosition chainTarget gen = runEffect $
+ chain step leaps Chain {..} gen
+ >-> Pipes.take n
+ >-> Pipes.mapM_ print
+ where
+ chainScore = lTarget chainTarget chainPosition
+ chainTunables = Nothing
+
+-- A Markov chain driven by the Metropolis transition operator.
+chain
+ :: (Num (IxValue (t Double)), Traversable t
+ , FunctorWithIndex (Index (t Double)) t, Ixed (t Double)
+ , PrimMonad m, IxValue (t Double) ~ Double)
+ => Double
+ -> Int
+ -> Chain (t Double) b
+ -> Gen (PrimState m)
+ -> Producer (Chain (t Double) b) m ()
+chain step leaps = loop where
+ loop state prng = do
+ next <- lift (MWC.sample (execStateT (hamiltonian step leaps) state) prng)
+ yield next
+ loop next prng
+
+-- | A Hamiltonian transition operator.
+hamiltonian
+ :: (Num (IxValue (t Double)), Traversable t
+ , FunctorWithIndex (Index (t Double)) t, Ixed (t Double), PrimMonad m
+ , IxValue (t Double) ~ Double)
+ => Double -> Int -> Transition m (Chain (t Double) b)
+hamiltonian e l = do
+ Chain {..} <- get
+ r0 <- lift (for chainPosition (const MWC.standard))
+ zc <- lift (MWC.uniform :: PrimMonad m => Prob m Double)
+ let (q, r) = leapfrogIntegrator chainTarget e l (chainPosition, r0)
+ perturbed = nextState chainTarget (chainPosition, q) (r0, r) zc
+ perturbedScore = lTarget chainTarget perturbed
+ put (Chain chainTarget perturbedScore perturbed chainTunables)
+
+-- Calculate the next state of the chain.
+nextState
+ :: (Foldable s, Foldable t, FunctorWithIndex (Index (t Double)) t
+ , FunctorWithIndex (Index (s Double)) s, Ixed (s Double)
+ , Ixed (t Double), IxValue (t Double) ~ Double
+ , IxValue (s Double) ~ Double)
+ => Target b
+ -> (b, b)
+ -> (s Double, t Double)
+ -> Double
+ -> b
+nextState target position momentum z
+ | z < pAccept = snd position
+ | otherwise = fst position
+ where
+ pAccept = acceptProb target position momentum
+
+-- Calculate the acceptance probability of a proposed moved.
+acceptProb
+ :: (Foldable t, Foldable s, FunctorWithIndex (Index (t Double)) t
+ , FunctorWithIndex (Index (s Double)) s, Ixed (t Double)
+ , Ixed (s Double), IxValue (t Double) ~ Double
+ , IxValue (s Double) ~ Double)
+ => Target a
+ -> (a, a)
+ -> (s Double, t Double)
+ -> Double
+acceptProb target (q0, q1) (r0, r1) = exp . min 0 $
+ auxilliaryTarget target (q1, r1) - auxilliaryTarget target (q0, r0)
+
+-- A momentum-augmented target.
+auxilliaryTarget
+ :: (Foldable t, FunctorWithIndex (Index (t Double)) t
+ , Ixed (t Double), IxValue (t Double) ~ Double)
+ => Target a
+ -> (a, t Double)
+ -> Double
+auxilliaryTarget target (t, r) = f t - 0.5 * innerProduct r r where
+ f = lTarget target
+
+innerProduct
+ :: (Num (IxValue s), Foldable t, FunctorWithIndex (Index s) t, Ixed s)
+ => t (IxValue s) -> s -> IxValue s
+innerProduct xs ys = Foldable.sum $ gzipWith (*) xs ys
+
+-- A container-generic zipwith.
+gzipWith
+ :: (FunctorWithIndex (Index s) f, Ixed s)
+ => (a -> IxValue s -> b) -> f a -> s -> f b
+gzipWith f xs ys = imap (\j x -> f x (fromMaybe err (ys ^? ix j))) xs where
+ err = error "gzipWith: invalid index"
+
+-- The leapfrog or Stormer-Verlet integrator.
+leapfrogIntegrator
+ :: (Num (IxValue (f Double)), Num (IxValue (t Double))
+ , FunctorWithIndex (Index (f Double)) t
+ , FunctorWithIndex (Index (t Double)) f
+ , Ixed (f Double), Ixed (t Double)
+ , IxValue (f Double) ~ Double
+ , IxValue (t Double) ~ Double)
+ => Target (f Double)
+ -> Double
+ -> Int
+ -> (f Double, t (IxValue (f Double)))
+ -> (f Double, t (IxValue (f Double)))
+leapfrogIntegrator target e l (q0, r0) = go q0 r0 l where
+ go q r 0 = (q, r)
+ go q r n = go q1 r1 (pred n) where
+ (q1, r1) = leapfrog target e (q, r)
+
+-- A single leapfrog step.
+leapfrog
+ :: (Num (IxValue (f Double)), Num (IxValue (t Double))
+ , FunctorWithIndex (Index (f Double)) t
+ , FunctorWithIndex (Index (t Double)) f
+ , Ixed (t Double), Ixed (f Double)
+ , IxValue (f Double) ~ Double, IxValue (t Double) ~ Double)
+ => Target (f Double)
+ -> Double
+ -> (f Double, t (IxValue (f Double)))
+ -> (f Double, t (IxValue (f Double)))
+leapfrog target e (q, r) = (qf, rf) where
+ rm = adjustMomentum target e (q, r)
+ qf = adjustPosition e (rm, q)
+ rf = adjustMomentum target e (qf, rm)
+
+adjustMomentum
+ :: (Functor f, Num (IxValue (f Double))
+ , FunctorWithIndex (Index (f Double)) t, Ixed (f Double))
+ => Target (f Double)
+ -> Double
+ -> (f Double, t (IxValue (f Double)))
+ -> t (IxValue (f Double))
+adjustMomentum target e (q, r) = r .+ ((0.5 * e) .* g q) where
+ g = fromMaybe err (glTarget target)
+ err = error "adjustMomentum: no gradient provided"
+
+adjustPosition
+ :: (Functor f, Num (IxValue (f Double))
+ , FunctorWithIndex (Index (f Double)) t, Ixed (f Double))
+ => Double
+ -> (f Double, t (IxValue (f Double)))
+ -> t (IxValue (f Double))
+adjustPosition e (r, q) = q .+ (e .* r)
+
+-- Scalar-vector product.
+(.*) :: (Num a, Functor f) => a -> f a -> f a
+z .* xs = fmap (* z) xs
+
+-- Vector addition.
+(.+)
+ :: (Num (IxValue t), FunctorWithIndex (Index t) f, Ixed t)
+ => f (IxValue t)
+ -> t
+ -> f (IxValue t)
+(.+) = gzipWith (+)
+
diff --git a/Numeric/Sampling/Hamiltonian.hs b/Numeric/Sampling/Hamiltonian.hs
@@ -1,143 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-
--- | Hamiltonian Monte Carlo. See, ex: Neal (2012)
--- http://arxiv.org/pdf/1206.1901.pdf.
-
-module Numeric.Sampling.Hamiltonian (
- Target(..)
- , Tunables(..)
- , Options(..)
- , hamiltonian
- , hmc
- ) where
-
-import Control.Applicative
-import Control.Monad
-import Control.Monad.Primitive
-import Control.Monad.Trans.Class
-import Control.Monad.Trans.State.Strict
-import qualified Data.Foldable as Foldable
-import Data.Sampling.Types
-import Data.Vector (Vector)
-import qualified Data.Vector as V
-import System.Random.MWC
-import System.Random.MWC.Distributions
-
-data Tunables = Tunables {
- stepSize :: !Double -- ^ step size for a given proposal
- , leapfrogs :: !Int -- ^ number of leapfrog steps to take
- } deriving Show
-
-data Options = Options {
- epochs :: !Int -- ^ number of epochs to iterate the chain
- , initial :: !Parameters -- ^ start location
- } deriving Show
-
-type Particle = (Parameters, Parameters)
-
-type Transition m = StateT Parameters m Parameters
-
--- | The Hamiltonian transition operator.
-hamiltonian
- :: PrimMonad m
- => Target
- -> Tunables
- -> Gen (PrimState m)
- -> Transition m
-hamiltonian target tunables g = do
- q0 <- get
- r0 <- V.replicateM (V.length q0) (lift $ standard g)
- zc <- lift $ uniform g
- let (q, r) = leapfrogIntegrator target tunables (q0, r0)
- next = nextState target (q0, q) (r0, r) zc
- put next
- return next
-{-# INLINE hamiltonian #-}
-
--- | The leapfrog or Stormer-Verlet integrator.
-leapfrogIntegrator :: Target -> Tunables -> Particle -> Particle
-leapfrogIntegrator target tunables (q0, r0) = go q0 r0 l where
- l = leapfrogs tunables
- go q r 0 = (q, r)
- go q r n =
- let (q1, r1) = leapfrog target tunables (q, r)
- in go q1 r1 (pred n)
-{-# INLINE leapfrogIntegrator #-}
-
--- | A single iteration of the leapfrog integrator.
-leapfrog :: Target -> Tunables -> Particle -> Particle
-leapfrog target tunables (q, r) = (qf, rf) where
- rm = adjustMomentum target tunables (q, r)
- qf = adjustPosition tunables (rm, q)
- rf = adjustMomentum target tunables (qf, rm)
-{-# INLINE leapfrog #-}
-
--- | Adjust momentum according to a half-leapfrog step.
-adjustMomentum :: Target -> Tunables -> Particle -> Parameters
-adjustMomentum target tunables (q, r) = r .+ ((0.5 * e) .* g q) where
- e = stepSize tunables
- g = glTarget target
-{-# INLINE adjustMomentum #-}
-
--- | Adjust position according to a half-leapfrog step.
-adjustPosition :: Tunables -> Particle -> Parameters
-adjustPosition tunables (r, q) = q .+ (e .* r) where
- e = stepSize tunables
-{-# INLINE adjustPosition #-}
-
--- | Scalar-vector multiplication.
-(.*) :: Double -> Parameters -> Parameters
-z .* xs = (* z) <$> xs
-{-# INLINE (.*) #-}
-
--- | Scalar-vector addition.
-(.+) :: Parameters -> Parameters -> Parameters
-xs .+ ys = V.zipWith (+) xs ys
-{-# INLINE (.+) #-}
-
--- | The next state of the Markov chain.
-nextState
- :: Target
- -> Particle
- -> Particle
- -> Double
- -> Parameters
-nextState target position momentum z
- | z < pAccept = snd position
- | otherwise = fst position
- where
- pAccept = acceptProb target position momentum
-{-# INLINE nextState #-}
-
--- | A target augmented by momentum auxilliary variables.
-auxilliaryTarget :: Target -> Particle -> Double
-auxilliaryTarget target (t, r) = f t - 0.5 * innerProduct r r where
- f = lTarget target
-{-# INLINE auxilliaryTarget #-}
-
--- | The acceptance probability of a move.
-acceptProb :: Target -> Particle -> Particle -> Double
-acceptProb target (q0, q1) (r0, r1) = exp . min 0 $
- auxilliaryTarget target (q1, r1) - auxilliaryTarget target (q0, r0)
-{-# INLINE acceptProb #-}
-
--- | Simple inner product.
-innerProduct :: Parameters -> Parameters -> Double
-innerProduct xs ys = V.sum $ V.zipWith (*) xs ys
-{-# INLINE innerProduct #-}
-
--- | Run a chain of HMC.
-hmc
- :: PrimMonad m
- => Target
- -> Tunables
- -> Options
- -> Gen (PrimState m)
- -> m [Parameters]
-hmc target tunables options g = do
- let h = hamiltonian target tunables g
- n = epochs options
- q = initial options
- evalStateT (replicateM n h) q
-{-# INLINE hmc #-}
-
diff --git a/README.md b/README.md
@@ -1,16 +1,31 @@
-# hasty-hamiltonian [](http://travis-ci.org/jtobin/hasty-hamiltonian)
+# hasty-hamiltonian
-Speedy, gradient-based traversal through parameter space. See the *examples* folder for example usage.
+[](http://travis-ci.org/jtobin/hasty-hamiltonian)
+[](http://hackage.haskell.org/package/hasty-hamiltonian)
-Install notes:
+Speedy, gradient-based traversal through parameter space.
- cabal sandbox init
- cabal install -j --only-dep
+Exports a `mcmc` function that prints a trace to stdout, as well as a
+`hamiltonian` transition operator that can be used more generally.
-To build the example,
+If you don't want to calculate your gradients by hand you can use the handy
+[ad](https://hackage.haskell.org/package/ad) library for automatic
+differentiation.
- cabal install -j --only-dep --enable-tests
- cabal build -j hasty-examples
+ import Numeric.AD (grad)
+ import Numeric.MCMC.Hamiltonian
-You can then use `./dist/build/hasty-examples/hasty-examples` to see an example trace of a chain wandering over the [Rosenbrock density](http://en.wikipedia.org/wiki/Rosenbrock_function).
+ target :: RealFloat a => [a] -> a
+ target [x0, x1] = negate ((x0 + 2 * x1 - 7) ^ 2 + (2 * x0 + x1 - 5) ^ 2)
+
+ gTarget :: [Double] -> [Double]
+ gTarget = grad target
+
+ booth :: Target [Double]
+ booth = Target target (Just gTarget)
+
+ main :: IO ()
+ main = withSystemRandom . asGenIO $ mcmc 10000 0.05 20 [0, 0] booth
+
+
diff --git a/hasty-hamiltonian.cabal b/hasty-hamiltonian.cabal
@@ -1,17 +1,42 @@
name: hasty-hamiltonian
-version: 0.3.0.0
+version: 1.1.0
synopsis: Speedy traversal through parameter space.
homepage: http://jtobin.github.com/hasty-hamiltonian
-license: BSD3
+license: MIT
license-file: LICENSE
author: Jared Tobin
maintainer: jared@jtobin.ca
category: Numeric
build-type: Simple
-cabal-version: >=1.18
+cabal-version: >=1.10
Description:
-
- Speedy gradient-based traversal through parameter space.
+ Gradient-based traversal through parameter space.
+ .
+ This implementation of HMC algorithm uses 'lens' as a means to operate over
+ generic indexed traversable functors, so you can expect it to work if your
+ target function takes a list, vector, map, sequence, etc. as its argument.
+ .
+ If you don't want to calculate your gradients by hand you can use the
+ handy <https://hackage.haskell.org/package/ad ad> library for automatic
+ differentiation.
+ .
+ Exports a 'mcmc' function that prints a trace to stdout, as well as a
+ 'hamiltonian' transition operator that can be used more generally.
+ .
+ > import Numeric.AD (grad)
+ > import Numeric.MCMC.Hamiltonian
+ >
+ > target :: RealFloat a => [a] -> a
+ > target [x0, x1] = negate ((x0 + 2 * x1 - 7) ^ 2 + (2 * x0 + x1 - 5) ^ 2)
+ >
+ > gTarget :: [Double] -> [Double]
+ > gTarget = grad target
+ >
+ > booth :: Target [Double]
+ > booth = Target target (Just gTarget)
+ >
+ > main :: IO ()
+ > main = withSystemRandom . asGenIO $ mcmc 10000 0.05 20 [0, 0] booth
Source-repository head
Type: git
@@ -19,30 +44,30 @@ Source-repository head
library
default-language: Haskell2010
+ ghc-options:
+ -Wall
exposed-modules:
- Numeric.Sampling.Hamiltonian
+ Numeric.MCMC.Hamiltonian
build-depends:
- base
- , mtl
- , mcmc-types
- , mwc-random
+ base < 5
+ , ghc-prim
+ , mcmc-types >= 1.0.1
+ , mwc-probability >= 1.0.1
+ , lens
+ , pipes
, primitive
, transformers
- , vector
- ghc-options:
- -Wall
-Test-Suite hasty-examples
+Test-suite booth
type: exitcode-stdio-1.0
- hs-source-dirs: examples
- main-is: Example.hs
+ hs-source-dirs: test
+ main-is: Booth.hs
default-language: Haskell2010
ghc-options:
- -threaded -rtsopts -Wall -O2
+ -rtsopts
build-depends:
ad
- , base
+ , base < 5
+ , mwc-probability >= 1.0.1
, hasty-hamiltonian
- , mwc-random >= 0.13.2.0
- , vector
diff --git a/stack.yaml b/stack.yaml
@@ -0,0 +1,7 @@
+flags: {}
+packages:
+- '.'
+extra-deps:
+ - mcmc-types-1.0.1
+ - mwc-probability-1.0.1
+resolver: lts-3.3
diff --git a/test/Beale.hs b/test/Beale.hs
@@ -0,0 +1,23 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+import Numeric.AD (grad)
+import Numeric.MCMC.Hamiltonian
+
+target :: RealFloat a => [a] -> a
+target [x0, x1]
+ | and [x0 >= -4.5, x0 <= 4.5, x1 >= -4.5, x1 <= 4.5]
+ = negate $ (1.5 - x0 + x0 * x1) ^ 2
+ + (2.25 - x0 + x0 * x1 ^ 2) ^ 2
+ + (2.625 - x0 + x0 * x1 ^ 3) ^ 2
+ | otherwise = - (1 / 0)
+
+gTarget :: [Double] -> [Double]
+gTarget = grad target
+
+beale :: Target [Double]
+beale = Target target (Just gTarget)
+
+main :: IO ()
+main = withSystemRandom . asGenIO $ mcmc 10000 0.05 20 [0, 0] beale
+
+
diff --git a/test/Booth.hs b/test/Booth.hs
@@ -0,0 +1,17 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+import Numeric.AD (grad)
+import Numeric.MCMC.Hamiltonian
+
+target :: RealFloat a => [a] -> a
+target [x0, x1] = negate ((x0 + 2 * x1 - 7) ^ 2 + (2 * x0 + x1 - 5) ^ 2)
+
+gTarget :: [Double] -> [Double]
+gTarget = grad target
+
+booth :: Target [Double]
+booth = Target target (Just gTarget)
+
+main :: IO ()
+main = withSystemRandom . asGenIO $ mcmc 10000 0.05 20 [0, 0] booth
+
diff --git a/test/Rosenbrock.hs b/test/Rosenbrock.hs
@@ -0,0 +1,17 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+import Numeric.AD (grad)
+import Numeric.MCMC.Hamiltonian
+
+target :: RealFloat a => [a] -> a
+target [x0, x1] = negate (5 *(x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+
+gTarget :: [Double] -> [Double]
+gTarget = grad target
+
+rosenbrock :: Target [Double]
+rosenbrock = Target target (Just gTarget)
+
+main :: IO ()
+main = withSystemRandom . asGenIO $ mcmc 3000 0.05 20 [0, 0] rosenbrock
+
