Skip to content

Instantly share code, notes, and snippets.

@mstksg
Created May 13, 2013 16:07
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mstksg/5569453 to your computer and use it in GitHub Desktop.
Save mstksg/5569453 to your computer and use it in GitHub Desktop.
import System.Random
import Control.Monad.Random
import Control.Monad
metroStep :: (RandomGen g) => (Double -> Double) -> Double -> Double -> Rand g Double
metroStep p step x = do
dx <- getRandomR (-step, step)
let
newx = x + dx
p0 = p x
p1 = p newx
weightselect = do
let odds = exp (p1 - p0)
select <- getRandom
return (if (select < odds) then newx else x)
if (p1 > p0)
then return newx
else weightselect
metroChain :: (RandomGen g) =>
(Double -> Double)
-> Double
-> Double
-> Int
-> Rand g Double
metroChain p step start n = (iterate mstep (return start)) !! n
where
mstep a = do
xval <- a
metroStep p step xval
main = do
let
p x = exp (-1*x) * (if (x>0) then 1 else 0)
res <- evalRandIO (replicateM 100 (metroChain p 0.1 0.5 10000))
mapM_ print res
@tel
Copy link

tel commented May 13, 2013

Here's a refactoring that uses mtl transformer stacks.

{-# LANGUAGE FlexibleContexts #-}

import System.Random

import Control.Monad
import Control.Monad.Random
import Control.Monad.State

type World a = StateT Double (Rand StdGen) a

pureMetro :: (Double -> Double) -> (Double, Double) -> Double -> Double
pureMetro p (x, newx) select =
  if (p1 > p0) || (select < exp (p1 - p0)) then newx else x
  where (p0, p1) = (p x, p newx)

metroStep :: (MonadRandom m, MonadState Double m)
             => (Double -> Double) -> Double -> m ()
metroStep p step = do
  dx <- getRandomR (-step, step)
  select <- getRandom
  modify $ \x -> (pureMetro p (x, x + dx) select)

runWorld :: Double -> World a -> IO a
runWorld x m = evalRandIO . (`evalStateT` x) $ m 

-- | Iterates a state monad, collecting the sequence of states
chain :: MonadState s m => Int -> m a -> m [s]
chain n step = replicateM n (step >> gets id)

main = do
  let p x = exp (-1*x) * (if (x>0) then 1 else 0)
  chain <- runWorld 0.5 (chain 1000 $ metroStep p 0.1)
  print (drop 200 chain)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment