1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

{-# 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)