Newcomb's problem in LazyPPL

This is an implementation of Newcombs problem using LazyPPL, a Haskell library for probabilistic programming. It uses the Metropolis-Hasting algorithm to obtain a numeric approximation of probability.

{-# LANGUAGE ExtendedDefaultRules #-}
 
import LazyPPL
import Distr
import Data.List
 
-- Lazy PPL example. We avoid plotting libraries.
model :: Meas Bool
model = do
  -- Prior: 
  -- 1:6 is Sunday.
  sunday <- sample $ bernoulli (1/7)
  -- Update: 
  -- On Sunday buses are distributed as Poisson(3) in one hour.
  -- On other days, buses are distributed as Poisson(10) in one hour.
  -- I have observed 4 buses in one hour.
  score $ poissonPdf (if sunday then 3 else 10) 4
  -- Posterior: 
  -- What is the probability of Sunday?
  return sunday
 
probSunday :: IO Double
probSunday = do
  xws <- mh 1 model
  let resolution = 3000
  let trues = length $ filter (== True) $ map fst $ take resolution xws
  return $ (fromIntegral trues :: Double) / (fromIntegral resolution :: Double)
 
 
 
-- PARTIAL MARKOV
observe :: Bool -> Meas ()
observe fact = score $ if fact then 0.999 else 0.001
 
-- NEWCOMB'S PROBLEM
newcomb :: Bool -> Meas Double
newcomb x = do
   -- Prior
   actor     <- sample $ bernoulli (1/2)
   predictor <- sample $ bernoulli (1/2)
   -- Observations
   observe $ actor == x
   observe $ predictor == actor
   -- Posterior
   return $ outcome actor predictor
 where
    outcome :: Bool -> Bool -> Double
    -- oneboxing
    outcome True  True  = 1000
    outcome True  False = 0
    -- twoboxing
    outcome False True  = 1001
    outcome False False = 1
 
 
execDouble :: Meas Double -> IO Double
execDouble model = do
  xws <- mh 1 model
  let resolution = 10000
  return 
    $ (\x -> x /(fromIntegral resolution)) 
    $ sum 
    $ map fst (take resolution xws)
 
-- probNewcomb :: IO 
probNewcomb x = do
  xws <- mh 1 (newcomb x)
  let resolution = 10000
  return 
    $ (\x -> x /(fromIntegral resolution)) 
    $ sum 
    $ map fst (take resolution xws)