A Metropolis sampler.

Metropolis.hs

module Metropolis where

import Control.Monad
import Control.Monad.Primitive
import System.Random.MWC as MWC
import System.Random.MWC.Distributions as MWC

propose :: [Double] -> Gen RealWorld -> IO [Double]
propose location gen = traverse (perturb gen) location where
  perturb rng m = MWC.normal m 1 rng

moveProbability :: ([Double] -> Double) -> [Double] -> [Double] -> Double
moveProbability altitude current proposed =
    whenNaN 0 (exp (min 0 (altitude proposed - altitude current)))
  where
    whenNaN val x
      | isNaN x   = val
      | otherwise = x

decide :: [Double] -> [Double] -> Double -> Gen RealWorld -> IO [Double]
decide current proposed prob rng = do
  accept <- MWC.bernoulli prob rng
  return $
    if   accept
    then proposed
    else current

metropolis :: ([Double] -> Double) -> [Double] -> Gen RealWorld -> IO [Double]
metropolis altitude current rng = do
  proposed <- propose current rng
  let prob = moveProbability altitude current proposed
  decide current proposed prob rng

chain :: Int -> ([Double] -> Double) -> [Double] -> Gen RealWorld -> IO [[Double]]
chain n altitude origin rng = loop n [origin] where
  loop j history@(current:_)
    | j <= 0    = return history
    | otherwise = do
        next <- metropolis altitude current rng
        loop (j - 1) (next:history)

---

landscape :: [Double] -> Double
landscape xs = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1)
  where
    x0 = xs !! 0
    x1 = xs !! 1

main :: IO ()
main = do
  rng <- MWC.createSystemRandom
  let origin = [-0.2, 0.3]
  trace <- chain 1000 landscape origin rng
  mapM_ print trace