Sampling.hs

{- | Example of sampling

Two samplers are availables : the 'discreteAncestralSampler' and the 'gibbsSampler'.
Only the 'gibbsSampler' can be used with evidence.

In this example, we have a very simple network.

@
    simple :: (['TDV' Bool],'SBN' 'CPT')
    simple = 'runBN' $ do 
        a <- 'variable' \"a\" ('t' :: Bool)
        b <- 'variable' \"b\" ('t' :: Bool) 
--        
        'proba' a '~~' [0.4,0.6]
        'cpt' b [a] '~~' [0.8,0.2,0.2,0.8]
--    
        return [a,b]
@

This network is representing a sensor b. We observe the value of b and we want to infer the value of a.

We use the 'gibbsSampler' for this with an initial period of 200 samples which are dropped. The 'gibbsSampler' is
generate a stream of samples. From this stream, we need to compute a probability distribution. For this, we use
the 'samplingHistograms' histogram function which is generating a list : the probability values of each vertex.

@
    let (vars\@[a,b],exampleG) = simple
    n <- 'runSampling' 5000 200 ('gibbsSampler' exampleG [b '=:' True])
    let h = 'samplingHistograms' n
    print $ h
@

Then, we compare this result with the exact one we get with a junction tree.

@
    let jt = 'createJunctionTree' 'nodeComparisonForTriangulation' exampleG
        jt' = 'changeEvidence' [b '=:' True] jt
    mapM_ (\x -> print . 'posterior' jt' $ [x]) vars
@

We can also use the 'discreteAncestralSampler' to compute the posterior but it is not supporting the use of evidence in this
version. The syntax is similar.

@
    n <- 'runSampling' 500 ('discreteAncestralSampler' exampleG)
@

-}
module Bayes.Examples.Sampling(
    -- * Test function
	  testSampling
) where

import Bayes.Sampling 
import System.Random.MWC.CondensedTable
import qualified Data.Vector as V
import Bayes.Factor
import Bayes
import Bayes.FactorElimination
import Data.Maybe(fromJust)
import Bayes.BayesianNetwork
import Bayes.Factor.CPT
import Statistics.Sample.Histogram


simple :: ([TDV Bool],SBN CPT)
simple = runBN $ do 
    a <- variable "a" (t :: Bool)
    b <- variable "b" (t :: Bool) 
    
    proba a ~~ [0.4,0.6]
    cpt b [a] ~~ [0.8,0.2,0.2,0.8]

    return [a,b]


testSampling = do
    let (vars@[a,b],exampleG) = simple
        jt = createJunctionTree nodeComparisonForTriangulation exampleG
    --n <- runSampling 500 (discreteAncestralSampler exampleG)
    putStrLn "The bayesian network"
    print exampleG
    putStrLn "\nThe values of the Bayesian network"
    printGraphValues exampleG

    n <- runSampling 5000 200 (gibbsSampler exampleG [b =: True])
    let h = samplingHistograms 2 n
    putStrLn "\nThe result of the sampling"
    print $ h

    n <- runSampling 5000 200 (gibbsMCMCSampler exampleG [b =: True])
    let h = samplingHistograms 2 n
    putStrLn "\nThe result of the sampling with MCMC"
    print $ h

    let jt' = changeEvidence [b =: True] jt
    putStrLn "\nThe result of the junction tree inference"
    mapM_ (\x -> print . posterior jt' $ [x]) vars