hbayes 0.4.1 → 0.5

raw patch · 20 files changed

+1450/−91 lines, 20 filesdep +gammadep +mwc-randomdep +statisticsdep ~QuickCheckdep ~arraydep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: gamma, mwc-random, statistics

Dependency ranges changed: QuickCheck, array, base, containers, mtl, pretty

API changes (from Hackage documentation)

Files

• Bayes.hs +39/−5
• Bayes/BayesianNetwork.hs +1/−1
• Bayes/Continuous.hs +330/−0
• Bayes/EM.hs +65/−0
• Bayes/Examples/ContinuousSampling.hs +118/−0
• Bayes/Examples/EMTest.hs +118/−0
• Bayes/Examples/Sampling.hs +98/−0
• Bayes/Examples/Tutorial.hs +23/−23
• Bayes/Factor.hs +2/−1
• Bayes/Factor/CPT.hs +17/−4
• Bayes/Factor/MaxCPT.hs +2/−3
• Bayes/Factor/PrivateCPT.hs +7/−5
• Bayes/FactorElimination.hs +9/−7
• Bayes/InfluenceDiagram.hs +9/−4
• Bayes/Network.hs +3/−0
• Bayes/PrivateTypes.hs +86/−14
• Bayes/Sampling.hs +490/−0
• Bayes/Test/CompareEliminations.hs +14/−14
• Bayes/Test/ReferencePatterns.hs +1/−1
• hbayes.cabal +18/−9

@@ -2,17 +2,23 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MultiParamTypeClasses #-}-{- | Discrete Bayesian Network Library.+{- | Bayesian Network Library.  It is a very preliminary version. It has only been tested on very simple examples where it worked. It should be considered as experimental and not used in any production work. -Look at the "Bayes.Examples" and "Bayes.Examples.Tutorial" in this package +* Look at the "Bayes.Examples" and "Bayes.Examples.Tutorial" in this package  to see how to use the library. -In "Bayes.Examples.Influence" you'll find additional examples about influence diagrams.+* In "Bayes.Examples.Influence" you'll find additional examples about influence diagrams. +* In "Bayes.Examples.Sampling" there are some explanations about the samplers for discrete networks.++* "Bayes.Examples.EMTest" is explaining learning with expectation / maximization.++* "Bayes.Examples.ContinuousSampling" is showing an example of sampling with a continuous network.+ -} module Bayes(   -- * Graph@@ -21,6 +27,7 @@   , UndirectedGraph(..)   , DirectedGraph(..)   , FoldableWithVertex(..)+  , FunctorWithVertex(..)   , NamedGraph(..)   -- ** Graph Monad   , GraphMonad@@ -35,6 +42,7 @@   , rootNode   , parentNodes   , childrenNodes+  , markovBlanket   -- ** Support functions for Graph constructions   , Vertex   , Edge @@ -86,9 +94,8 @@ -- | An implementation of the BayesianNetwork using the simple graph and no value for the edges type SBN f = DirectedSG () f --- | Bayesian network. g must be a directed graph and f a factor-type BayesianNetwork g f = g () f + instance Arbitrary (DirectedSG String String) where   arbitrary = do      let createVertex g i = do @@ -279,6 +286,16 @@     ingoing :: g a b -> Vertex -> Maybe [Edge]     outgoing :: g a b -> Vertex -> Maybe [Edge] ++-- | Return the Markov blanket of a node +markovBlanket :: DirectedGraph g => g a b -> Vertex -> [Vertex]+markovBlanket g c = +  let p = parentNodes g c+      children = childrenNodes g c +      parentOfChildrens = concatMap (parentNodes g) children +  in +  nub (c:p ++ children ++ parentOfChildrens)+ -- | Return the parents of a node parentNodes :: DirectedGraph g => g a b -> Vertex -> [Vertex] parentNodes g v = maybe [] id $ do @@ -418,6 +435,23 @@  instance Functor (SimpleGraph local edge) where    fmap f (SP em vm nm) = SP em (IM.map (\(l,d) -> (l, f d)) vm) nm++class FunctorWithVertex g where +   fmapWithVertex :: (Vertex -> a -> b) -> g c a -> g c b +   fmapWithVertexM :: Monad m => (Vertex -> a -> m b) -> g c a -> m (g c b) +++instance FunctorWithVertex (SimpleGraph local) where +    fmapWithVertex f (SP em vm nm) = SP em (IM.mapWithKey (\k (l,d) -> (l, f (Vertex k) d)) vm) nm+    fmapWithVertexM f (SP em vm nm) = do +      let l = IM.toList vm+          g f (k,(l,d)) = do +            r <- f (Vertex k) d +            return (k,(l,r))+      rl <- mapM (g f) $ l+      let vm' = IM.fromList rl+      return $ SP em vm' nm+  instance F.Foldable (SimpleGraph local edge) where   foldr f c (SP _ vm _) = IM.foldr (\(_,d) s -> f d s) c vm

@@ -142,7 +142,7 @@   -- an enumeration   (.==.) :: d -> v -> LE  -instance Instantiable d v => Testable d v where +instance Instantiable d v DVI => Testable d v where    (.==.) a b = LETest (a =: b)  infixl 8 .==.

@@ -0,0 +1,330 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{- Bayesian Networks with continuous variables++-}+module Bayes.Continuous(+    -- * Types +    -- ** For the graph description+      CNMonad(..)+    , CV+    , DN+    , VariableName(..)+    , BayesianVariable(..)+    -- ** For the sampling +    , Distri +    , ContinuousNetwork(..) +    , ContinuousSample(..)+    , InstantiationValue(..)+    , CVI +    -- * Network creation functions+    , uniform +    , normal+    , beta+    , beta'+    , exponential+    , execCN+    , runCN +    , evalCN+    -- * Sampling functions+    , runSampling +    , continuousMCMCSampler+    , (=:)+    -- * Result statistics+    , histogram+    , samplingHistograms +    ) where++import Bayes.Sampling+import Bayes+import Data.Monoid+import Control.Applicative+import Bayes.PrivateTypes+import Bayes.Network(NetworkMonad(..),runGraph,execGraph,evalGraph)+import Control.Monad.State.Strict+import Data.Maybe(fromJust)+import Control.Monad+import Control.Monad.Reader +import Bayes.Sampling(runSampling, samplingHistograms,gibbsMCMCSampler, ContinuousNetwork(..), ContinuousSample(..)) +import Math.Gamma(gamma)++-- | The Bayesian monad+type CNMonad a = GraphMonad DirectedSG () Distri a+++instance Monoid (DistributionSupport (Double,Double)) where +    mempty = BoundedSupport (0.0,0.0)+    (BoundedSupport (xa,ya)) `mappend` (BoundedSupport (xb,yb)) = BoundedSupport (min xa xb, max ya yb)+    (BoundedSupport _) `mappend` (Unbounded a) = Unbounded a +    (Unbounded a) `mappend` (BoundedSupport _) = Unbounded a +    (Unbounded a) `mappend` (Unbounded b) = Unbounded (a+b)++newtype RN = RN (Reader (Sample DirectedSG CVI) Double) ++instance Num RN where +   (RN a) + (RN b) = RN $ liftA2 (+) a b+   (RN a) - (RN b) = RN $ liftA2 (-) a b+   (RN a) * (RN b) = RN $ liftA2 (*) a b+   abs (RN a) = RN $ liftA abs a +   signum (RN a) = RN $ liftA signum a +   fromInteger a = RN (return (fromInteger a))++instance Fractional RN where +   (RN a) / (RN b) = RN $ liftA2 (/) a b+   fromRational a = RN (return (fromRational a))++-- | An expression which can be a constant, variable or formula.+-- In case it is a variable, it can be used as a 'BayesianVariable'+-- or instantiated as an 'Instantiable' type.+-- Otherwise you'll get an error+data DN = DN !([CV]) !RN++instance BayesianVariable DN where +    vertex (DN [] _) = error "No vertex for this expression"+    vertex (DN (v:[]) _) = vertex v +    vertex (DN l _) = error "Too many vertices for this expression"++dependencies :: DN -> [CV]+dependencies (DN l _) = l ++instance Instantiable DN Double CVI where +  (=:) (DN [] _) x = error "No variable is related to this expression"+  (=:) (DN (a:[]) _) x = a =: x+  (=:) (DN l _) x = error "Too many variables are related to this expression"++value :: DN -> Reader (Sample DirectedSG CVI) Double +value (DN _ (RN v)) = v++instance Num DN where +    (DN da na) + (DN db nb) = DN (da `mappend`db) (na + nb)+    (DN da na) - (DN db nb) = DN (da `mappend`db) (na - nb)+    (DN da na) * (DN db nb) = DN (da `mappend`db) (na * nb)+    abs (DN a b) = DN a (abs b)+    signum (DN a b) = DN a (signum b)+    fromInteger i = DN [] (fromInteger i)++instance Fractional DN where+    (DN da na) / (DN db nb) = DN (da `mappend`db) (na / nb)+    fromRational i = DN [] (fromRational i)+++{-++Graph creation++-}++-- | This class is used to simplify the network description. Variable+-- names can be optional. In that later case, () must be used instead of a name.+class VariableName m where +    mkVariable :: m -> CNMonad CV ++instance VariableName String where +    mkVariable = variable ++instance VariableName () where +    mkVariable _ = unamedVariable++-- | Create a new named Bayesian variable if not found.+-- Otherwise, return the found one.+addVariableIfNotFound :: String -> CNMonad Vertex+addVariableIfNotFound vertexName = graphNode vertexName (D (CV (Vertex 0)) (Distri (return mempty) (\s -> return 0.0)))++-- | Create a new unamed variable+unamedVariable :: CNMonad CV+unamedVariable = do +  ((namemap,count),g) <- get +  let vertexName = "unamed" ++ show count+  va <- addVariableIfNotFound vertexName+  return (CV va)+++whenJust Nothing _ = return ()+whenJust (Just i) f = f i >> return ()++-- | Set the node of a bayesian network under creation+setBayesianNode :: CV -> Distri -> CNMonad ()+setBayesianNode (CV v) newValue = do+  (aux,oldGraph) <- get+  let newGraph = changeVertexValue v newValue oldGraph+ +  whenJust newGraph $ \nvm -> do+     put $! (aux, nvm)++-- | Create an edge between two vertex of the Bayesian network+(<--) :: CV -> CV -> CNMonad ()+(CV va) <-- (CV vb) = newEdge vb va mempty++++node :: CV -> DN+node v = DN [v] (RN $ do +    r <- ask +    return . instantiationValue . fromJust . vertexValue r $ (vertex v)+    )++-- | Define a Bayesian variable (name and bounds)+variable :: String -- ^ Variable name+         -> CNMonad CV+variable name = do+  va <- addVariableIfNotFound name+  return (CV va)++cpt :: CV -> [CV] -> CNMonad CV+cpt node conditions = do+  mapM_ (node <--) (reverse conditions)+  return node++{-++Distributions++-}++-- | Exponential distribution +exponential :: VariableName s +            => s +            -> DN +            -> CNMonad DN +exponential s na = do +    v <- mkVariable s+    let la = dependencies na +    cpt v la +    let bound = do +           lambda <- value na+           return (Unbounded (1.0 / lambda))+        result = D v (Distri bound $ \inst -> do +                    lambda <- value na+                    let x = instantiationValue inst +                    if x < 0+                      then +                        return 0.0 +                      else +                        return $ lambda * exp(-lambda * x)+                    )+    setBayesianNode v result +    return (node v)++-- | Beta distribution +beta :: VariableName s +     => s +     -> DN +     -> DN +     -> CNMonad DN +beta s na nb = do +    v <- mkVariable s+    let la = dependencies na +        lb = dependencies nb +    let l = la ++ lb+    cpt v l+    let bound = return (BoundedSupport (0.0,0.1) )+        result = D v (Distri bound $ \inst -> do +                    a <- value na+                    b <- value nb+                    let x = instantiationValue inst +                    return $ gamma(a+b)/gamma(a)/gamma(b)*x**(a-1)*(1-x)**(b-1)+                    )+    setBayesianNode v result +    return (node v)++-- | Beta' distribution +beta' :: VariableName s +     => s +     -> DN +     -> DN +     -> CNMonad DN +beta' s na nb = do +    v <- mkVariable s+    let la = dependencies na +        lb = dependencies nb +    let l = la ++ lb+    cpt v l+    let bound = do +              a <- value na+              b <- value nb+              let std = a * (a + b - 1) -- Not really the standard deviation but related to it+                                        -- enough for generating new samples+              return (Unbounded std)+        result = D v (Distri bound $ \inst -> do +                    a <- value na+                    b <- value nb+                    let x = instantiationValue inst +                    if x <= 0 +                      then +                        return 0.0+                      else+                        return $ gamma(a+b)/gamma(a)/gamma(b)*x**(a-1)*(1+x)**(-a-b)+                    )+    setBayesianNode v result +    return (node v)++-- | Uniform dstribution+uniform :: VariableName s +        => s -- ^ Variable name+        -> DN -- ^ Min bound+        -> DN -- ^ Max bound+        -> CNMonad DN+uniform s na nb = do +    v <- mkVariable s+    let la = dependencies na +        lb = dependencies nb +    let l = la ++ lb+    cpt v l+    let bound = do +        a <- value na +        b <- value nb +        return (BoundedSupport (a,b) )+    let result = D v (Distri bound $ \inst -> do +                    a <- value na+                    b <- value nb+                    let d = 1.0 / (b - a)+                    if instantiationValue inst < a || instantiationValue inst > b then return 0.0 else return d)+    setBayesianNode v result +    return (node v)+++-- | Normal distribution+normal :: VariableName s +       => s -- ^ Variable name+       -> DN -- ^ Average+       -> DN -- ^ Standard deviation+       -> CNMonad DN+normal s mean std = do+    v <- mkVariable s+    let la = dependencies mean +        lb = dependencies std +    let l = la ++ lb+    cpt v l+    let bound = do +        s <- value std +        return (Unbounded s )+    let result = D v (Distri bound $ \inst -> do +                    m <- value mean+                    s <- value std+                    let x = instantiationValue inst +                    return (1.0 / (sqrt (2 * pi)*s) * exp(-(x-m)*(x-m)/(2*s*s))+                           ))+    setBayesianNode v result +    return (node v)++{-+ +Graph creation from the Monad.  ++-}++-- | Create a  network using the simple graph implementation+-- The initialized nodes are replaced by the value.+-- Returns the monad values and the built graph.+runCN :: CNMonad a -> (a,ContinuousNetwork)+runCN = runGraph++-- | Create a  network but only returns the monad value.+-- Mainly used for testing.+execCN :: CNMonad a -> ContinuousNetwork+execCN = execGraph++-- | Create a bayesian network but only returns the monad value.+-- Mainly used for testing.+evalCN :: CNMonad a -> a+evalCN = evalGraph

@@ -0,0 +1,65 @@+{- | Expectation / Maximization to learn Bayesian network values++-}+module Bayes.EM( +    -- * Learning function+      learnEM+    ) where ++import Bayes+import Bayes.Factor +import Bayes.Sampling +import Bayes.Factor.CPT+import Bayes.FactorElimination+import Data.Maybe(fromJust)+++sumG :: (FunctorWithVertex g, Graph g) +     => Sample g (CPT,CPT) +     -> Sample g (CPT,CPT) +     -> Sample g (CPT,CPT)+sumG ga gb = +  let sumNode vertex (xa,xb) = +        let (ya,yb) = fromJust . vertexValue gb $ vertex +        in +        (cptSum [xa,ya], cptSum [xb,yb])+  in +  fmapWithVertex sumNode ga ++divideG :: Vertex +        -> (CPT,CPT) +        -> CPT+divideG _ (a,b) = cptDivide a b ++-- | Learn network values from samples using the expectation / maximization algorithm.+learnEM :: (FunctorWithVertex g, NamedGraph g, FoldableWithVertex g, DirectedGraph g) +        => [[DVI]] -- ^ Samples+        -> BayesianNetwork g CPT -- ^ Start network+        -> BayesianNetwork g CPT -- ^ Network with new values learnt from the samples+learnEM samples startG = +    let jt = createJunctionTree nodeComparisonForTriangulation startG +        results = map (computeSample startG jt) samples+    in +    fmapWithVertex divideG (foldl1 sumG results)++computeSample :: (FunctorWithVertex g, Graph g) +              => BayesianNetwork g CPT +              -> JunctionTree CPT +              -> [DVI] +              -> Sample g (CPT,CPT)+computeSample startG jt s = fmapWithVertex (computeNode startG jt s) startG++computeNode :: Graph g +            => BayesianNetwork g CPT +            -> JunctionTree CPT +            -> [DVI] +            -> Vertex +            -> CPT +            -> (CPT,CPT)+computeNode g jt samples vertex _ = +    let jt' = changeEvidence samples jt+        f@(main:l) = factorVariables (fromJust . vertexValue g $ vertex)+        a = fromJust $ posterior jt' f +        b = if null l then factorFromScalar 1.0 else fromJust $ posterior jt' l +    in +    (a,b)

@@ -0,0 +1,118 @@+{- | Sampling example with continuous distributions++Continuous networks can't be handled by any of the functions defined for the discrete networks.+So, instead of using exact inference algorithms like the junction trees, sampling method have to be used.++In this example, we want to estimate a parameter which is measured by noisy sensors.++There are 'nbSensors' available. They are described with a 'normal' distribution centered on the value of+the unknown parameters and with a standard deviation of 0.1.++The unknown parameter is described with a 'uniform' distribution bounded by 1.0 and 2.0.++First, we describe the sensor:++@+sensor :: 'DN' -> 'CNMonad' 'DN' +sensor p = do +    'normal' () p 0.1 +@++It is just a 'normal' distribution. The mean of this distribution is the parameters p. This parameter has special type 'DN'.+All expressions used to build the continuous bayesian network are using values of type 'DN'. A value of type 'DN' can either+represent a constant, a variable or an expression.++If the sensor was biased, we could write:++@+    'normal' ()  (p + 0.2) 0.1+@++The Bayesian network describing the measurement process is given by:++@+    test = 'runCN' $ do+      a <- 'uniform' \"a\" 1.0 2.0 -- Unknown parameter+      sensors <- sequence (replicate 'nbSensors' (sensor a))+      return (a:sensors)+@++We are connecting 'nbSensors' nodes corresponding to the 'nbSensors' measurements.  In real life it can either be different sensors+or the same one used several times (assuming the value of the parameter is not dependent on time).++Now, as usual in all the examples of this package, we get the bayesian graph and a list of variables used+to compute some posterior or define some evidences++@+    debugcn = do +        let ((a:sensors), testG) = test+@++Then, we generate some random measurements and create the evidences++@+    g <- create +    measurements <- sequence . replicate  nbSensors $ (MWC.normal 1.5 0.1 g)+    let evidence = zipWith (=:) sensors measurements+@++Evidence has type 'CVI' and is created with the assigment operator '=:' .++Now, we generate some samples to estimate the posterior distributions.++@+    n <- 'runSampling' 10000 200 ('continuousMCMCSampler' testG evidence)+@++This function is generating a sequence of graphs ! We are not interested in the sensor values. They are known and fixed+since they have been measured. So, we extract the value of the parameter.++@+    let samples = map (\g -> 'instantiationValue' . fromJust . 'vertexValue' g $ ('vertex' a)) n+@++And with the samples for the parameters we can compute an histogram and get an approximation of the posterior.++@+    let samples = map (\g -> 'instantiationValue' . fromJust . 'vertexValue' g $ ('vertex' a)) n+        h = 'histogram' 6 samples +    print h+@++We see in the histogram that the estimated value is around 1.5.++-}+module Bayes.Examples.ContinuousSampling(+	  nbSensors +	, sensor +	, test +	, debugcn+	) where+import Bayes+import Bayes.Continuous +import qualified System.Random.MWC.Distributions as MWC(normal)+import System.Random.MWC(GenIO,create)+import Data.Maybe(fromJust)++nbSensors = 10 ++sensor :: DN -> CNMonad DN +sensor p = do +    normal () p 0.1 ++test = runCN $ do+  a <- uniform "a" 1.0 2.0 -- Unknown parameter+  sensors <- sequence (replicate nbSensors (sensor a))+  return (a:sensors)++debugcn = do +    let ((a:sensors), testG) = test+    g <- create +    measurements <- sequence . replicate  nbSensors $ (MWC.normal 1.5 0.1 g)+    let evidence = zipWith (=:) sensors measurements+    n <- runSampling 10000 200 (continuousMCMCSampler testG evidence)+    let samples = map (\g -> instantiationValue . fromJust . vertexValue g $ (vertex a)) n+        h = histogram 6 samples +    print h++

@@ -0,0 +1,118 @@+{- | Test of learning++In this example, two networks are used : 'simple' which is the reference and 'wrong' which is a wrong start.+The goal is to use test patterns to learn the right 'simple' network from 'wrong'. Only the values are learnt.+The topology of both networks is the same.++@+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]+@++and 'wrong' where the probability for a is wrong.++@+wrong :: (['TDV' Bool],'SBN' 'CPT')+wrong = 'runBN' $ do +    a <- 'variable' \"a\" ('t' :: Bool)+    b <- 'variable' \"b\" ('t' :: Bool) +--    +    'proba' a '~~' [0.2,0.8]+    'cpt' b [a] '~~' [0.8,0.2,0.2,0.8]+--+    return [a,b]+@++So, the first thing to do is generate test patterns. We are using the 'discreteAncestralSampler' for this. This function is+generating a sequence of graphs. We are just interested in the values. So, we get the values with 'allVertexValues'.++@+generatePatterns :: IO [[DVI]]+generatePatterns = do +    let (vars\@[a,b],exampleG) = simple+    r <- 'runSampling' 5000 0 ('discreteAncestralSampler' exampleG)+    return (map 'allVertexValues' r)+@++Once we have the data, we can try to learn the network:++@+emTest = do +  samples <- generatePatterns +  let (_,simpleG) = simple +      (_,wrongG) = wrong +  print simpleG +  'printGraphValues' simpleG+  'printGraphValues' wrongG+--+  'printGraphValues' ('learnEM' samples wrongG)+@++First, we display the topology of the graph and the values for the reference graph and the wrong one.+Then, we use the 'learnEM' function to learn a new network from the samples. And, we print the new values to check.++-}+module Bayes.Examples.EMTest(+    -- * Test function+    emTest+  ) 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+import Bayes.EM++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]++wrong :: ([TDV Bool],SBN CPT)+wrong = runBN $ do +    a <- variable "a" (t :: Bool)+    b <- variable "b" (t :: Bool) +    +    proba a ~~ [0.2,0.8]+    cpt b [a] ~~ [0.8,0.2,0.2,0.8]++    return [a,b]++setJust vertex v = Just v ++generatePatterns :: IO [[DVI]]+generatePatterns = do +    let (vars@[a,b],exampleG) = simple+    r <- runSampling 5000 0 (discreteAncestralSampler exampleG)+    return (map allVertexValues r)++emTest = do +  samples <- generatePatterns +  let (_,simpleG) = simple +      (_,wrongG) = wrong +  print simpleG +  printGraphValues simpleG+  printGraphValues wrongG++  +  printGraphValues (learnEM samples wrongG)+ 

@@ -0,0 +1,98 @@+{- | 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

@@ -67,7 +67,7 @@ Once the junction tree has been computd, it can be used to compute several marginals:  @-    'posterior' jt rain+    'posterior' jt [rain] @  The function is called posterior and will compute posterior only when solme evidence has@@ -77,7 +77,7 @@  @     let jt' = 'updateEvidence' [wet '=:' True] jt -    'posterior' jt' rain+    'posterior' jt' [rain] @  If you want to compute the posterior for a combination of variables, you have two possibilities : either going back to the@@ -99,12 +99,12 @@ @     let jt4 = 'changeEvidence' [wet '=:' True] jt      print \"Posterior Marginal : probability of rain and road slippery if grass wet\"-    let m = 'posterior' jt4 roadandrain+    let m = 'posterior' jt4 [roadandrain]     print m --     let jt5 = 'changeEvidence' [wet '=:' True, sprinkler '=:' False] jt      print \"Posterior Marginal : probability of rain and road slippery if grass wet and srinkler not used\"-    let m = 'posterior' jt5 roadandrain+    let m = 'posterior' jt5 [roadandrain] @  /Inferences with an imported network/@@ -146,11 +146,11 @@ @     let jtcancer = 'createJunctionTree' 'nodeComparisonForTriangulation' cancer ---    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer $ x)) [varA,varB,varC,varE]+    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer $ [x])) [varA,varB,varC,varE] --     print \"UPDATED EVIDENCE\"     let jtcancer' = 'updateEvidence' [varD '=:' Present] jtcancer -    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer' $ x)) [varA,varB,varC,varE]+    mapM_ (\x -> putStrLn (show x) >> (print . 'posterior' jtcancer' $ [x])) [varA,varB,varC,varE] @  The '=:' operator will check that the assignment is type compatible because varD is a typed discrete variable of type 'TDV' Coma.@@ -208,13 +208,13 @@  @     print \"Sensor 90%\"-    print $ posterior ('changeFactor' (theNewFactor 0.9) jt') a+    print $ posterior ('changeFactor' (theNewFactor 0.9) jt') [a] --     print \"Sensor 50%\"-    print $ posterior ('changeFactor' (theNewFactor 0.5) jt') a+    print $ posterior ('changeFactor' (theNewFactor 0.5) jt') [a] --     print \"Sensor 10%\"-    print $ posterior ('changeFactor' (theNewFactor 0.1) jt') a+    print $ posterior ('changeFactor' (theNewFactor 0.1) jt') [a] @  -}@@ -257,7 +257,7 @@   (varmap,perso) <- exampleDiabete   let jtperso = createJunctionTree nodeComparisonForTriangulation perso       cho0 = fromJust . Map.lookup "cho_0" $ varmap-  print $ posterior jtperso cho0+  print $ posterior jtperso [cho0] #endif  miscTest s = do @@ -270,7 +270,7 @@   print "FACTOR ELIMINATION"   let post (v,name) = do          putStrLn name -        print $ posterior jtperso v+        print $ posterior jtperso [v]   mapM_ post  (zip l names)    print "VARIABLE ELIMINATION"@@ -294,15 +294,15 @@     let [varA,varB,varC,varE] = fromJust $ mapM (flip Map.lookup varmap) ["A","B","C","E"]         jtcancer = createJunctionTree nodeComparisonForTriangulation cancer         varD = tdv (fromJust $ Map.lookup "D" varmap) :: TDV Coma-    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer $ x)) [varA,varB,varC,varE]+    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer $ [x])) [varA,varB,varC,varE]      print "UPDATED EVIDENCE : Coma present"     let jtcancer' = changeEvidence [varD =: Present] jtcancer -    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varE]+    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ [x])) [varA,varB,varC,varE]      print "UPDATED EVIDENCE : Coma absent"     let jtcancer' = changeEvidence [varD =: Absent] jtcancer -    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ x)) [varA,varB,varC,varE]+    mapM_ (\x -> putStrLn (show x) >> (print . posterior jtcancer' $ [x])) [varA,varB,varC,varE] #endif   -- | Display of factors generated by the logical keyword@@ -322,13 +322,13 @@         theNewFactor x = fromJust $ se seNode a x -- x % success for the sensor         jt' = changeEvidence [seNode =: True] jt     print "Sensor 90%"-    print $ posterior (changeFactor (theNewFactor 0.9) jt') a+    print $ posterior (changeFactor (theNewFactor 0.9) jt') [a]      print "Sensor 50%"-    print $ posterior (changeFactor (theNewFactor 0.5) jt') a+    print $ posterior (changeFactor (theNewFactor 0.5) jt') [a]      print "Sensor 10%"-    print $ posterior (changeFactor (theNewFactor 0.1) jt') a+    print $ posterior (changeFactor (theNewFactor 0.1) jt') [a]  -- | Inferences with the standard network inferencesOnStandardNetwork = do@@ -360,38 +360,38 @@     print "FACTOR ELIMINATION"     putStrLn ""     print "Prior Marginal : probability of rain"-    let m = posterior jt rain+    let m = posterior jt [rain]     print m     putStrLn ""      let jt' = changeEvidence [wet =: True] jt       print "Posterior Marginal : probability of rain if grass wet"-    let m = posterior jt' rain+    let m = posterior jt' [rain]     print m     putStrLn ""      let jt'' = changeEvidence [] jt'     print "Prior Marginal : probability of rain"-    let m = posterior jt rain+    let m = posterior jt [rain]     print m     putStrLn ""      let jt3 = changeEvidence [wet =: True, sprinkler =: True] jt'     print "Posterior Marginal : probability of rain if grass wet and sprinkler used"-    let m = posterior jt3 rain+    let m = posterior jt3 [rain]     print m     putStrLn ""      let jt4 = changeEvidence [wet =: True] jt      print "Posterior Marginal : probability of rain and road slippery if grass wet"-    let m = posterior jt4 roadandrain+    let m = posterior jt4 [roadandrain]     print m     putStrLn ""      let jt5 = changeEvidence [wet =: True, sprinkler =: False] jt      print "Posterior Marginal : probability of rain and road slippery if grass wet and srinkler not used"-    let m = posterior jt5 roadandrain+    let m = posterior jt5 [roadandrain]     print m     putStrLn "" 

@@ -16,6 +16,7 @@  -- * Set of variables   , Set(..)  , BayesianDiscreteVariable(..)+ , BayesianVariable(..)  -- * Implementation  , Vertex(..)  -- ** Discrete variables and instantiations@@ -24,10 +25,10 @@  , DVSet(..)  , DVI  , DVISet+ , InstantiationValue(..)  , tdvi  , tdv  , setDVValue- , instantiationValue  , instantiationVariable  , variableVertex  , (=:)

@@ -111,7 +111,6 @@     mkValue = id     scale = (*)     multiply = (*)-    elementUnit = 1.0     divide _ 0 = 0     divide a b = a / b      elementSum = (+)@@ -136,18 +135,32 @@          evidenceFrom = _evidenceFrom -    factorProduct = _factorProduct multiply+    factorProduct = _factorProduct (Op (factorFromScalar 1.0) 1.0 multiply)          factorProjectOut _ f@(Scalar v) = f     factorProjectOut s f= cptFactorProjectOutWith (sum . map fst) s f  --  Divie two CPT cptDivide :: CPT -> CPT -> CPT-cptDivide a b = _factorProduct divide [a, b]+cptDivide a b   | isScalarFactor a && isScalarFactor b = +                    let va = factorValue a [] +                        vb = factorValue b [] +                    in +                    factorFromScalar (divide va vb)+                | isScalarFactor a = +                    let va = factorValue a [] +                    in +                    factorScale va b +                | isScalarFactor b = +                    let vb = factorValue b [] +                    in +                    if vb == 0.0 then factorFromScalar 0.0 else factorScale (1.0 / vb) a +                | otherwise =+                    _factorProduct (Op (factorFromScalar 1.0) 1.0 divide) [b,a]   cptSum :: [CPT] -> CPT-cptSum = _factorProduct elementSum+cptSum = _factorProduct (Op (factorFromScalar 0.0) 0.0 elementSum)  instance IsBucketItem CPT where      scalarItem = isScalarFactor

@@ -29,7 +29,6 @@     multiply (da,[]) (db,lb) = (da*db,lb)     multiply (da,la) (db,[]) = (da*db,la)     multiply (da,la) (db,lb) = (da*db,[xa ++ xb | xa <- la, xb <- lb])-    elementUnit = (1.0,[])     divide _ _ = error "It does not make sense to divide MAXCPT elements"     elementSum _ _ = error "It does not make sense to sum MAXCPT elements" @@ -37,7 +36,7 @@ instantiationProba (a,b) = fst a   maximization :: [((Double,PossibleInstantiations), DVISet)] -> (Double,PossibleInstantiations)-maximization [] = elementUnit +maximization [] = (1.0, [])  maximization l =      let ((d,possible),newInst) = maximumBy (compare `on` instantiationProba) l         addTo [] l = l@@ -65,7 +64,7 @@      evidenceFrom = _evidenceFrom -    factorProduct = _factorProduct multiply+    factorProduct = _factorProduct (Op (factorFromScalar 1.0) (1.0,[]) multiply)          factorProjectOut _ f@(Scalar v) = f     factorProjectOut s f = cptFactorProjectOutWith maximization s f

@@ -11,6 +11,7 @@    -- ** Private functions    , FactorElement(..)    , DecisionFactor(..)+   , Operator(..)    , _emptyFactor    , _factorVariables    , _isScalarFactor@@ -130,7 +131,6 @@     multiply :: a -> a -> a     divide :: a -> a -> a     elementSum ::  a -> a -> a-    elementUnit :: a  _isUsingSameVariablesAs :: PrivateCPT v a -> PrivateCPT v a -> Bool _isUsingSameVariablesAs (Table va _ _) (Table vb _ _) = va == vb@@ -210,17 +210,19 @@ {-# INLINE _factorValue #-} _factorValue f d = doubleValue (privateFactorValue f d) +data Operator f a = Op f a (a -> a -> a)+ _factorProduct :: (FactorElement a, GV.Vector v a, Factor (PrivateCPT v a)) -               => (a -> a -> a)+               => Operator (PrivateCPT v a) a                -> [PrivateCPT v a]                 -> PrivateCPT v a-_factorProduct _ [] = factorFromScalar 1.0-_factorProduct op l = +_factorProduct (Op fUnit opUnit op) [] = fUnit+_factorProduct (Op fUnit opUnit op) l =      let nakedVars = L.foldr1 union . map factorVariables $ l         allVars = DVSet nakedVars         (scalars,cpts) = partition isScalarFactor l         stridesFromCPT (Table d _ _) = indexStridesFor (DVSet d) allVars-        elementProduct [] = elementUnit+        elementProduct [] = opUnit         elementProduct l = foldr1 op l         ps = elementProduct . map (flip privateFactorValue undefined) $ scalars         cptsStrides = map stridesFromCPT cpts

@@ -308,10 +308,12 @@   -- | Compute the marginal posterior (if some evidence is set on the junction tree)-  -- otherwise compute just the marginal prior.-posterior :: (BayesianDiscreteVariable dv, Factor f, IsBucketItem f) => JunctionTree f -> dv -> Maybe f+-- otherwise compute just the marginal prior. The set of variables must be included inside a cluster+-- for thr algorithm to work. So, most of the cases, it will be used to compute the posterior of just+-- one variable.+posterior :: (BayesianDiscreteVariable dv, Factor f, IsBucketItem f) => JunctionTree f -> [dv] -> Maybe f posterior t someDv = -  let v = dv someDv+  let v = map dv someDv   in   case snd $ traverseTree (findClusterFor v) Nothing t of      Nothing -> Nothing@@ -319,19 +321,19 @@                   d = maybe (factorFromScalar 1.0) id $ downMessage t =<< (nodeParent t c)                   u = map (upMessage t) (nodeChildren t c)                   allFactors = d:u ++ f ++ e-                  variablesToRemove = (nub (concatMap factorVariables allFactors)) \\ [v]-                  unNormalized = marginal allFactors variablesToRemove [v] []+                  variablesToRemove = (nub (concatMap factorVariables allFactors)) \\ v+                  unNormalized = marginal allFactors variablesToRemove v []               in                Just $ factorDivide unNormalized (factorNorm unNormalized)  -- | Find a cluster containing the variable-findClusterFor :: DV +findClusterFor :: [DV]                -> Maybe Cluster                -> Cluster -- ^ Current cluster                -> NodeValue f -- ^ Current value                -> Action (Maybe Cluster) (NodeValue f) findClusterFor dv s c@(Cluster sc) v = -  case Set.member dv sc of +  case Set.isSubsetOf (Set.fromList dv) sc of      False -> Skip s      True -> Stop (Just c) 

@@ -269,12 +269,14 @@ instance Show DEV where     show (DEV v d) = "D" ++ show v ++ "(" ++ show d ++ ")" +instance BayesianVariable DEV where +  vertex (DEV v _) = v+ instance BayesianDiscreteVariable DEV where    dimension (DEV _ d) = d    dv (DEV v d) = DV v d -  vertex (DEV v _) = v -instance Instantiable DEV Int where +instance Instantiable DEV Int DVI where    (=:) d@(DEV v dim) value = DVI (dv d) value  data PorD = P DV | D DEV deriving(Eq)@@ -288,13 +290,16 @@ instance ChanceVariable (TDV s) where    toDV = dv +instance BayesianVariable PorD where +    vertex (D d) = vertex d+    vertex (P p) = vertex p+ instance BayesianDiscreteVariable PorD where     dimension (D d) = dimension d     dimension (P p) = dimension p     dv (D x) = dv x     dv (P x) = dv x-    vertex (D d) = vertex d-    vertex (P p) = vertex p+  -- | Used to mix decision and chance variables and a same list p :: ChanceVariable c => c -> PorD

@@ -22,6 +22,9 @@ 	, runNetwork 	, execNetwork 	, evalNetwork+  , runGraph+  , execGraph+  , evalGraph 	, getCpt 	) where  

@@ -5,10 +5,13 @@ -} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-} module Bayes.PrivateTypes(   -- * Classes and types    BayesianDiscreteVariable(..)+ , BayesianVariable(..)  , Set(..)+ , BayesianNetwork(..)  -- * Variables  , DV(..)  , DVSet(..)@@ -16,11 +19,13 @@  , TDV  , tdv  , tdvi+ , CV(..)+ , CVI(..)  -- * Instantiations  , Instantiable(..)+ , InstantiationValue(..)  , DVI(..)  , setDVValue- , instantiationValue  , instantiationVariable  , fromDVSet  -- * Vertices, Graph@@ -38,6 +43,10 @@  , instantiationDetails  , instantiation  , allInstantiationsForOneVariable+ -- * For sampling + , Sample(..)+ , DistributionSupport(..)+ , DistributionF(..)  -- * Tests   , instantiationProp  ) where@@ -50,9 +59,32 @@ import System.Random(Random) import qualified Data.IntMap as IM import qualified Data.Map as M+import Control.Monad.Reader  +-- | Bayesian network. g must be a directed graph and f a factor+type BayesianNetwork g f = g () f +{-+For sampling+-} +++-- | A sample (graph of instantiations)+type Sample g a = BayesianNetwork g a++-- | Support of the distributiob p(x | abc ...) for x depending on the value of a,b,c ...+data DistributionSupport a = BoundedSupport !a +                           | Unbounded !Double -- ^ Standard deviation+                           deriving(Eq)++-- | Distribution function with a possible bound for the function target+-- The bound may depend on the current sample (because for unbounded we need a variance to be able+-- to cover enough of the support)+data DistributionF g bound inst = Distri !(Reader (Sample g inst) (DistributionSupport bound)) !(inst -> Reader (Sample g inst) Double)+++ {- 	Set -}@@ -128,24 +160,44 @@ instance Show Vertex where      show (Vertex v) = "v" ++ show v +-- | A Bayesian Variable is a variable part of Bayesian network and so which knows+-- its position : the vertex.+class BayesianVariable v where +  vertex :: v -> Vertex + -- | A discrete variable has a number of levels which is required to size the factors-class BayesianDiscreteVariable v where+class BayesianVariable v => BayesianDiscreteVariable v where     dimension :: v -> Int      dv :: v -> DV-    vertex :: v -> Vertex  -- | Get the minimum bound for a type getMinBound :: Bounded a => a -> a  getMinBound _ = minBound --+instance BayesianVariable Vertex where +  vertex v = v +   {-  Variables 	 -} +-- | A continuous variable+data CV = CV !Vertex deriving(Eq,Ord)++instance Show CV where +  show (CV v) = show v++instance BayesianVariable CV where +  vertex (CV v) = v ++-- | A continuous variable instantiation +data CVI = CVI !CV !Double deriving(Eq,Ord,Show)++instance BayesianVariable CVI where +  vertex (CVI v _) = vertex v+ -- | A discrete variable data DV = DV !Vertex !Int deriving(Eq,Ord) @@ -160,10 +212,12 @@ instance Show DV where     show (DV v d) = show v ++ "(" ++ show d ++ ")" +instance BayesianVariable DV where +    vertex (DV v _) = v+ instance BayesianDiscreteVariable DV where     dimension (DV _ d) = d     dv = id-    vertex (DV v _) = v  -- | A typed discrete variable data TDV s = TDV !Vertex !Int deriving(Eq,Ord)@@ -171,10 +225,12 @@ instance Show (TDV s) where     show (TDV v d) = show v +instance BayesianVariable (TDV s) where +    vertex (TDV v _) = v+ instance BayesianDiscreteVariable (TDV s) where     dimension (TDV _ d) = d     dv (TDV v nb) = DV v nb-    vertex (TDV v _) = v  -- | Typed discrete variable tdv :: DV -> TDV s @@ -228,16 +284,32 @@    -- | A set of variable instantiations type DVISet = [DVI] -class Instantiable d v where +class InstantiationValue i v | i -> v where +  instantiationValue :: i -> v +  toDouble :: i -> Double++instance InstantiationValue DVI Int where +  -- | Extract value of the instantiation+  instantiationValue (DVI _ v) = v+  toDouble (DVI _ v) = fromIntegral v++instance InstantiationValue CVI Double where +  -- | Extract value of the instantiation+  instantiationValue (CVI _ v) = v+  toDouble (CVI _ v) = v++class Instantiable d v r | d -> r where    -- | Create a variable instantiation using values from   -- an enumeration-  (=:) :: d -> v -> DVI +  (=:) :: d -> v -> r  +instance Instantiable CV Double CVI where +  (=:) c x = CVI c x -instance (Bounded b, Enum b) => Instantiable DV b where+instance (Bounded b, Enum b) => Instantiable DV b DVI where   (=:) a b = setDVValue a (fromEnum b - fromEnum (getMinBound b)) -instance (Bounded b, Enum b) => Instantiable (TDV b) b where+instance (Bounded b, Enum b) => Instantiable (TDV b) b DVI where   (=:) (TDV v nb) b = setDVValue (DV v nb) (fromEnum b - fromEnum (getMinBound b))  @@ -246,18 +318,18 @@ setDVValue v a = DVI v a  +instance BayesianVariable DVI where +    vertex (DVI dv _) = vertex dv  instance BayesianDiscreteVariable DVI where     dimension (DVI v _) = dimension v     dv = instantiationVariable-    vertex (DVI dv _) = vertex dv  -- | Get the variables and their values with a type constraint instantiationDetails :: [DVI] -> (DVSet s, MultiIndex s) instantiationDetails l = (DVSet $ map instantiationVariable l, MultiIndex . V.fromList . map (instantiationValue) $ l) --- | Extract value of the instantiation-instantiationValue (DVI _ v) = v+  -- | Discrete variable from the instantiation instantiationVariable (DVI dv _) = dv

@@ -0,0 +1,490 @@+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ExistentialQuantification #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{- | Sampling++Samplers for Bayesian network inferences++-}+module Bayes.Sampling(+  -- * Types +     Sampler(..)+   , Sample(..)+  -- * Sampling+   , runSampling+   , topologicalOrder+  -- * Sampling schemes+   , discreteAncestralSampler+   , gibbsSampler+   , gibbsMCMCSampler+  -- * Sampling results+   , samplingHistograms+   , histogram+  -- * For continuous distributions +   , ContinuousNetwork(..) +   , ContinuousSample(..)+   , Distri(..)+   , continuousMCMCSampler+) where +++import Bayes+import Bayes.Factor+import Control.Monad.Reader +import Data.Maybe(fromJust)+import Bayes.PrivateTypes+import System.Random.MWC.CondensedTable+import System.Random.MWC.Distributions(normal)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as UV+import qualified Data.Vector.Unboxed.Mutable as MV(write,read)++import System.Random.MWC(GenIO,create,withSystemRandom,asGenIO,uniformR,Variate)++--import Debug.Trace ++--debug s a = trace ("DEBUG:\n" ++ s ++ show a ++ "\n") a++-- | A graph of sampling functions (parental sampler or markov blanket sampler)+-- It is used to generate new samples.+-- The vertex is implicit in the function definition+-- Which means that if a function is for updating vertex V then it will be the value of vertex+-- V in the sampler graph.+type SamplerGraph g a = BayesianNetwork g (Sample g a -> IO a)++-- | Sampling function for a graph.+-- Generate a new sample and update the state. +type SamplingScheme g b a = GenIO -> SamplerGraph g a -> b -> Sample g a -> IO (b,Sample g a)++-- | Sampler defining the behavior of a sampling algorithms (init value,+-- sample generation, how to select nodes in the grapg)+data Sampler g a = forall b . Sampler !b !(GenIO -> IO (Sample g a)) !(GenIO -> SamplerGraph g a) !(SamplingScheme g b a)+++data Distri = D !CV !(DistributionF DirectedSG (Double,Double) CVI)++type ContinuousNetwork = SBN Distri++type ContinuousSample = SBN CVI++-- | Get current value of function using current network values+distributionOnNetwork ::(DistributionF DirectedSG (Double,Double) CVI) -> Vertex -> Reader (Sample DirectedSG CVI) Double +distributionOnNetwork (Distri _ f) v = do +  r <- ask +  let inst = fromJust . vertexValue r $ v +  f inst++instance Graph g => Show (Sample g [DVI]) where +    show = show . allNodes++instance Graph g => Show (Sample g CVI) where +    show = show . allNodes++instance Graph g => Show (Sample g [(Double,Double,Double)]) where +    show = show . allNodes++{-++Ancestral sampling ++-}++-- | Return the vertices in topological order+topologicalOrder :: DirectedGraph g => g a b -> [Vertex] +topologicalOrder g = _topologicalOrder g [] + where+  _topologicalOrder g current = case rootNode g of +    Nothing -> reverse current+    Just r -> _topologicalOrder (removeVertex r g) (r:current)++-- | Generate a new sample for vertex based upon the sampler functions+-- and the current sample+generateNewSample :: Graph g => SamplerGraph g a  -> Sample g a -> Vertex -> IO (Sample g a)+generateNewSample sampler sample theVertex = do+   let vertexSampler = fromJust $ vertexValue sampler theVertex +   newValue <- vertexSampler sample +   return . fromJust . changeVertexValue theVertex newValue $ sample ++-- | Update all vertices in the current sample and generate a new sample.+-- It is done using ancestral sampling+ancestralSampling :: DirectedGraph g => SamplingScheme g [Vertex] a+ancestralSampling gen sampler nodes s  = do +  r <- foldM (generateNewSample sampler) s nodes+  return (nodes,r)++-- | Select a vertex and generate a new sample using gibbs sampling+gibbsSampling :: DirectedGraph g => SamplingScheme g (CondensedTable V.Vector Vertex) a+gibbsSampling gen samplerGraph table s = do +  aVertex <- genFromTable table gen +  r <- generateNewSample samplerGraph s aVertex+  return (table,r)++-- | Sequence a list of the SAME graphs but with different node values+sequenceG :: (FunctorWithVertex g, Graph g) => [g a b] -> g a [b]+sequenceG [] = error "Can't sequence an empty list of graphs"+sequenceG l@(h:_) = +  let setVertexList vertex oldValue = map (\g -> fromJust . vertexValue g $ vertex) l+  in +  fmapWithVertex setVertexList h ++_range :: Double -> Double -> [Double] -> (Double,Double)+_range mi ma [] = (mi,ma) +_range mi ma (h:t) | h < mi = _range h ma t +                   | h > ma = _range mi h t +                   | otherwise = _range mi ma t ++range :: [Double] -> (Double,Double)+range = _range (1/0) (-1/0)++-- | Compute the histogram of values+histogram :: Int -- ^ Bins +          -> [Double] -- ^ Samples +          -> [(Double,Double,Double)] -- ^ Histogram  (minBound,maxBound, value)+histogram bins values | bins <= 0 = error "Bins must be >= 1" +                      | bins == 1 =+  let (bmin,bmax) = range values +  in +  [(bmin,bmax,1.0)]+                      | otherwise = +  let (bmin,bmax) = range values +  in+  if bmax == bmin +    then +      [(bmin,bmax,1.0)]+    else +      let delta = (bmax - bmin) / fromIntegral bins+          v = UV.replicate bins 0.0+          orMAX a | a >= bins = bins - 1 +                  | otherwise = a+          countSample value = UV.modify (\s -> do +            --when (isNaN value) $ error "NaN"+            let index = orMAX $ floor ((value - bmin) / delta) +            --when (index < 0) $ error ("index == " ++ show index ++ show " " ++ show value) +            --when (index >= bins) $ error ("index == " ++ show index ++ show " " ++ show value)+            d <- MV.read s index +            MV.write s index (d+1)) +          v' = foldr countSample v values+          t = UV.sum v'+          scaled = UV.map (/ t) v'+          addBounds value i = (bmin + delta * i ,bmin + delta + delta*i,value)+      in +      zipWith addBounds (UV.toList scaled) [0..]++-- | Generate a graph of sampling histogram for each variable+-- So, for a vertex v we have the posterior values p(v)+samplingHistograms :: (InstantiationValue i v,BayesianVariable i, FunctorWithVertex g, Graph g) +                   => Int -- ^ Number of bins+                   -> [Sample g i] -- ^ Samples+                   -> Sample g [(Double,Double,Double)] -- ^ Histogram with start bins and bin content+samplingHistograms bins g' = +  let g = sequenceG g' +  in+  fmapWithVertex (createHistogram bins g) g+ where +  createHistogram bins g vertex l =  histogram bins . map toDouble $ l++-- | Sample a bayesian network using a given sampling scheme+runSampling :: (DirectedGraph g, FunctorWithVertex g) +            => Int -- ^ Number of used samples+            -> Int -- ^ Burnin samples before the used samples+            -> Sampler g a -- ^ Sampler+            -> IO [Sample g a] -- ^ List of generated samples+runSampling n b (Sampler initState startSample sampler samplingScheme) = do +    withSystemRandom . asGenIO $ \gen -> do+       start <- startSample gen+       l <- _runAncestralSampling gen (n+b) initState start []+       return (drop b $ l)+  where +    _runAncestralSampling gen n istate i r | n <= 1 = return $ reverse (i:r)+                                           | otherwise = do+                                               (newState,i') <- samplingScheme gen (sampler gen)  istate i+                                               _runAncestralSampling gen (n-1) newState i' (i':r)  +++{-++Discrete variable sampling++-}  ++-- | Gibbs sampling+gibbsSampler :: (Factor f, FunctorWithVertex g, DirectedGraph g) +             => BayesianNetwork g f -- ^ Bayesian network+             -> [DVI] -- ^ Evidence+             -> Sampler g DVI -- ^ Sampler+gibbsSampler g evidence = +  let isEvidence v = v `elem` (map vertex evidence)+      v = filter (not . isEvidence) (allVertices g)+      t = tableFromWeights $ V.fromList $ zip v (repeat 1.0)+  in+  Sampler t (gibbsInitValue evidence g) (gibbsGraphSampler g)  gibbsSampling++-- | Gibbs sampling+gibbsMCMCSampler :: (Factor f, FunctorWithVertex g, DirectedGraph g) +                 => BayesianNetwork g f -- ^ Bayesian network+                 -> [DVI] -- ^ Evidence+                 -> Sampler g DVI -- ^ Sampler+gibbsMCMCSampler g evidence = +  let isEvidence v = v `elem` (map vertex evidence)+      v = filter (not . isEvidence) (allVertices g)+      t = tableFromWeights $ V.fromList $ zip v (repeat 1.0)+  in+  Sampler t (gibbsInitValue evidence g) (gibbsMCMCGraphSampler g)  gibbsSampling++-- | Gibbs sampling for continuous network +continuousMCMCSampler :: ContinuousNetwork -- ^ Bayesian network+                      -> [CVI] -- ^ Evidence+                      -> Sampler DirectedSG CVI-- ^ Sampler+continuousMCMCSampler g evidence = +  let isEvidence v = v `elem` (map vertex evidence)+      v = filter (not . isEvidence) (allVertices g)+      t = tableFromWeights $ V.fromList $ zip v (repeat 1.0)+  in+  Sampler t (continuousInitValue evidence g) (gibbsContinuousMCMCGraphSampler g)  gibbsSampling++-- | Ancestral sampler which does not support evidence+discreteAncestralSampler :: (Factor f, FunctorWithVertex g, DirectedGraph g) => BayesianNetwork g f -> Sampler g DVI+discreteAncestralSampler g = +  let nodes = topologicalOrder g+  in+  Sampler nodes (discreteInitValue [] g) (discreteSampler g)  ancestralSampling+++-- | Distribution function from a factor+-- When we compute the distribution P(x | abc...), abc are given by the graph+-- but x is varied so we must not take its value fromthe graph.+-- So, x is an argument. +distributionFromFactor :: (DirectedGraph g,Factor f) => f -> DistributionF g DV DVI +distributionFromFactor f = +    let vars@(main:_) = factorVariables f+    in +    Distri (return $ BoundedSupport main) $ \current -> do +      sample <- ask+      let currentV = vertex main+          sample' = fromJust . changeVertexValue currentV current $ sample+      return $ factorValue f (variableInstantiations sample' vars)++-- | Distribution function from  markov blanket but not normalized+distributionFromMarkovBlanket :: (DirectedGraph g,Factor f) => BayesianNetwork g f -> f -> DistributionF g DV DVI +distributionFromMarkovBlanket g f = +  let vars@(main:_) = factorVariables f+    in +    Distri (return $ BoundedSupport main) $ \current -> do +      sample <- ask+      let currentV = vertex main +          sample' = fromJust . changeVertexValue currentV current $ sample+          children = childrenNodes sample currentV+          childrenF = map (fromJust . vertexValue g) children +          values = map (\aFactor -> factorValue aFactor (variableInstantiations sample' (factorVariables aFactor))) (f:childrenF)+      return $ product values++continuousDistributionFromMarkovBlanket :: ContinuousNetwork +                                        -> Distri +                                        -> DistributionF DirectedSG (Double,Double) CVI+continuousDistributionFromMarkovBlanket g (D cv d@(Distri b distri)) = Distri b $ \current -> do +      sample <- ask+      let currentV = vertex cv +          children = childrenNodes sample currentV+          dist (D _ d) = d+          childrenF = map (dist . fromJust . vertexValue g) children +      local (fromJust . changeVertexValue currentV current) $ do+        s' <- ask+        values <- mapM (uncurry distributionOnNetwork) ((d,currentV):zip childrenF children)+        return $ product $ values++-- | Value of distribution functions for all range values+forAllRange :: Graph g => Sample g DVI -> DistributionF g DV DVI -> [Double] +forAllRange g (Distri bound r) = +  let range = allInstantiationsForOneVariable theBound +  in +  map (flip runReader g . r) range+ where +  BoundedSupport theBound = runReader bound g+++-- | Value of distribution on Sample+distributionValue :: Graph g => DistributionF g bounds inst -> Sample g inst -> inst -> Double +distributionValue (Distri bound r) g v = runReader (r v) g+++-- | Probability distribution of A given some values for the parents in P(A | parents)+distribution :: Graph g => Sample g DVI -> DistributionF g DV DVI  -> CondensedTableV Int+distribution g f = +  let notNullProba (i,p) = p /= 0.0+      values = forAllRange g f+  in+  tableFromWeights $ V.fromList  $ filter notNullProba $ zip [0..] values++-- | Get the values of factor variables on the current instantiation+variableInstantiations :: (DirectedGraph g, BayesianVariable v) => Sample g a -> [v] -> [a]+variableInstantiations g vars = +  let vertices = map vertex vars+  in+  map (fromJust . vertexValue g) vertices++-- | Force the value of the vertex if it is part of the evidence otherwise+-- use the function+withEvidence :: (Monad m, BayesianVariable v, BayesianVariable vb) => [v] -> vb -> (vb -> m v) -> m v+withEvidence e v f = +  if (vertex v) `elem` (map vertex e) +    then +      return . head . filter (\x -> vertex v == vertex x) $ e+    else +      f v+++-- | Null instantiations (min bound for each variable)+nullInstantiation :: Factor f => [DVI] -> Vertex -> f -> DVI+nullInstantiation evidence v f  = setDVValue (factorMainVariable f) 0 ++-- | Random instantiation for a variable +randomInstantiation :: Factor f => GenIO -> [DVI] -> Vertex -> f -> IO DVI +randomInstantiation gen evidence vertex f = withEvidence evidence vertex $ \v -> do +  let main = factorMainVariable f+      t = tableFromWeights $ V.fromList $ zip [0..dimension main - 1] (repeat 1.0)+  r <- genFromTable t gen +  return (setDVValue main r)++continuousRandomInstantiation :: GenIO -> [CVI] -> Vertex -> Distri -> IO CVI +continuousRandomInstantiation gen evidence vertex (D cv (Distri b _)) = withEvidence evidence vertex $ \v -> do +   value <- uniformR (-10.0,10.0) gen :: IO Double+   return (cv =: value)++-- | Update one vertex in the current sample+updateOneSample :: (DirectedGraph g,Factor f) => GenIO -> Vertex -> f -> (Sample g DVI -> IO DVI)+updateOneSample gen v f = \s -> do+    let d = distributionFromFactor f+    newVal <- genFromTable (distribution s d) gen +    let v = factorMainVariable f +    return (setDVValue v newVal)++-- | Update one vertex in the current sample using gibbs transitions+gibbsUpdateOneSample :: (DirectedGraph g,Factor f) => BayesianNetwork g f -> GenIO -> Vertex -> f -> (Sample g DVI -> IO DVI)+gibbsUpdateOneSample g gen v f = +    let d = distributionFromMarkovBlanket g f+        v = factorMainVariable f +    in +    \s -> do+        newVal <- genFromTable (distribution s d) gen +        return (setDVValue v newVal)+++class SamplingBounds bound inter | bound -> inter where +  samplingBounds :: Num inter => bound -> (inter,inter)++class SamplingBounds bound inter => SampleGeneration var inst bound inter | var -> inst, var -> bound, var -> inter where+  generateSample :: Instantiable var inter inst => Sample g inst -> DistributionF g bound inst -> GenIO -> var -> IO inst++instance SamplingBounds DV Int where +  samplingBounds dv = (0, dimension dv - 1)++instance SampleGeneration DV DVI DV Int where+  generateSample s (Distri b _) gen var = do +      let theBound = runReader b s+      case theBound of +        BoundedSupport s -> do+            let (sa,sb) = samplingBounds s+            u <- uniformR (sa,sb) gen+            return $ setDVValue var u+        Unbounded v -> error "DVI are bounded"++instance SamplingBounds (Double,Double) Double where +  samplingBounds (sa,sb) = (sa,sb)++instance SampleGeneration CV CVI (Double,Double) Double where+  generateSample s (Distri b _) gen var = do +     let theBound = runReader b s+     case theBound of +       BoundedSupport s -> do +                let (sa,sb) = samplingBounds s+                u <- uniformR (sa,sb) gen+                return $ var =: u+       Unbounded sigma -> do +                n <- normal 0.0 sigma gen+                return $ var =: n++gibbsMCMCUpdateOneSample :: (DirectedGraph g,Factor f) +                         => BayesianNetwork g f +                         -> GenIO +                         -> Vertex +                         -> f +                         -> (Sample g DVI -> IO DVI)+gibbsMCMCUpdateOneSample g gen v f = +  let d = distributionFromMarkovBlanket g f+      main = factorMainVariable f+      v = vertex main+  in +  _gibbsMCMCUpdateOneSample d gen v main++gibbsContinuousMCMCUpdateOneSample :: ContinuousNetwork +                                   -> GenIO +                                   -> Vertex +                                   -> Distri+                                   -> (ContinuousSample -> IO CVI)+gibbsContinuousMCMCUpdateOneSample g gen v distri@(D cv _) = +  let d = continuousDistributionFromMarkovBlanket g distri+      v = vertex cv+  in +  _gibbsMCMCUpdateOneSample d gen v cv++-- | Generate a sample using MCMC method.+-- | Update one vertex in the current sample using gibbs transitions+_gibbsMCMCUpdateOneSample :: (DirectedGraph g,BayesianVariable var, Instantiable var inter inst, SampleGeneration var inst bound inter) +                          => DistributionF g bound inst+                          -> GenIO+                          -> Vertex+                          -> var+                          -> (Sample g inst -> IO inst)+_gibbsMCMCUpdateOneSample d@(Distri _ _)  gen v var = \s -> do+        newSample <- generateSample s d gen var+        let oldSample = fromJust . vertexValue s $ v+            oldProbaValue = distributionValue d s oldSample+            newProbaValue = distributionValue d s newSample+        if oldProbaValue == 0.0+          then +            return newSample +          else +            let r = newProbaValue / oldProbaValue +            in+            if (r >= 1.0) +              then +                return newSample +              else do +                z <- uniformR (0.0, 1.0)  gen+                if z < r +                  then +                    return newSample +                  else +                    return oldSample++++-- | Sampler graph where each vertex know how to update its value+discreteSampler :: (FunctorWithVertex g, Factor f, DirectedGraph g) => BayesianNetwork g f -> GenIO -> SamplerGraph g DVI+discreteSampler g gen = fmapWithVertex (updateOneSample gen) g++-- | Sampler graph where each vertex know how to update its value+gibbsGraphSampler :: (FunctorWithVertex g, Factor f, DirectedGraph g) => BayesianNetwork g f -> GenIO -> SamplerGraph g DVI+gibbsGraphSampler g gen = fmapWithVertex (gibbsUpdateOneSample g gen) g++-- | Sampler graph where each vertex know how to update its value+gibbsMCMCGraphSampler :: (FunctorWithVertex g, Factor f, DirectedGraph g) => BayesianNetwork g f -> GenIO -> SamplerGraph g DVI+gibbsMCMCGraphSampler g gen = fmapWithVertex (gibbsMCMCUpdateOneSample g gen) g++gibbsContinuousMCMCGraphSampler :: ContinuousNetwork+                                -> GenIO +                                -> SamplerGraph DirectedSG CVI+gibbsContinuousMCMCGraphSampler g gen = fmapWithVertex (gibbsContinuousMCMCUpdateOneSample g gen) g++-- | Start value+discreteInitValue :: (FunctorWithVertex g, Factor f, DirectedGraph g) => [DVI] -> BayesianNetwork g f -> GenIO -> IO (Sample g DVI) +discreteInitValue evidences g gen = fmapWithVertexM (randomInstantiation gen evidences) g++-- | Start value+gibbsInitValue :: (FunctorWithVertex g, Factor f, DirectedGraph g) => [DVI] -> BayesianNetwork g f -> GenIO -> IO (Sample g DVI) +gibbsInitValue = discreteInitValue++continuousInitValue ::[CVI] -> ContinuousNetwork -> GenIO -> IO ContinuousSample+continuousInitValue evidences g gen = fmapWithVertexM (continuousRandomInstantiation gen evidences) g

@@ -46,20 +46,20 @@     -- We then update the values to the right ones by using the right factor         jt = changeFactor theNewFactor jt' -    compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain])+    compareFactors "PRIOR FOR RAIN" (posterior jt [rain]) (priorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain])      let jt1 = changeEvidence [wet =: True] jt          jt2 = changeEvidence [wet =: True, sprinkler =: True] jt1  -    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 rain) +    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 [rain])           (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True])-    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 rain) +    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 [rain])           (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True, sprinkler =: True]) -    compareFactors "PRIOR FOR WINTER" (posterior jt winter) (priorMarginal exampleG [sprinkler,wet,road,rain,roadandrain] [winter])-    compareFactors "PRIOR FOR SPRINKLER" (posterior jt sprinkler) (priorMarginal exampleG [winter,wet,road,rain,roadandrain] [sprinkler])-    compareFactors "PRIOR FOR WET" (posterior jt wet) (priorMarginal exampleG [winter,sprinkler,road,rain,roadandrain] [wet])-    compareFactors "PRIOR FOR ROAD" (posterior jt road) (priorMarginal exampleG [winter,sprinkler,wet,rain,roadandrain] [road])+    compareFactors "PRIOR FOR WINTER" (posterior jt [winter]) (priorMarginal exampleG [sprinkler,wet,road,rain,roadandrain] [winter])+    compareFactors "PRIOR FOR SPRINKLER" (posterior jt [sprinkler]) (priorMarginal exampleG [winter,wet,road,rain,roadandrain] [sprinkler])+    compareFactors "PRIOR FOR WET" (posterior jt [wet]) (priorMarginal exampleG [winter,sprinkler,road,rain,roadandrain] [wet])+    compareFactors "PRIOR FOR ROAD" (posterior jt [road]) (priorMarginal exampleG [winter,sprinkler,wet,rain,roadandrain] [road])   -- | Compare that variable elemination and factor elimination are giving@@ -68,18 +68,18 @@ compareVariableFactor = do      let ([winter,sprinkler,rain,wet,road,roadandrain],exampleG) = example         jt = createJunctionTree nodeComparisonForTriangulation exampleG-    compareFactors "PRIOR FOR RAIN" (posterior jt rain) (priorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain])+    compareFactors "PRIOR FOR RAIN" (posterior jt [rain]) (priorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain])      let jt1 = changeEvidence [wet =: True] jt          jt2 = changeEvidence [wet =: True, sprinkler =: True] jt1  -    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 rain) +    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt1 [rain])           (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True])-    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 rain) +    compareFactors "POSTERIOR RAIN FOR WET" (posterior jt2 [rain])           (posteriorMarginal exampleG [winter,sprinkler,wet,road,roadandrain] [rain]  [wet =: True, sprinkler =: True]) -    compareFactors "PRIOR FOR WINTER" (posterior jt winter) (priorMarginal exampleG [sprinkler,wet,road,rain,roadandrain] [winter])-    compareFactors "PRIOR FOR SPRINKLER" (posterior jt sprinkler) (priorMarginal exampleG [winter,wet,road,rain,roadandrain] [sprinkler])-    compareFactors "PRIOR FOR WET" (posterior jt wet) (priorMarginal exampleG [winter,sprinkler,road,rain,roadandrain] [wet])-    compareFactors "PRIOR FOR ROAD" (posterior jt road) (priorMarginal exampleG [winter,sprinkler,wet,rain,roadandrain] [road])+    compareFactors "PRIOR FOR WINTER" (posterior jt [winter]) (priorMarginal exampleG [sprinkler,wet,road,rain,roadandrain] [winter])+    compareFactors "PRIOR FOR SPRINKLER" (posterior jt [sprinkler]) (priorMarginal exampleG [winter,wet,road,rain,roadandrain] [sprinkler])+    compareFactors "PRIOR FOR WET" (posterior jt [wet]) (priorMarginal exampleG [winter,sprinkler,road,rain,roadandrain] [wet])+    compareFactors "PRIOR FOR ROAD" (posterior jt [road]) (priorMarginal exampleG [winter,sprinkler,wet,rain,roadandrain] [road]) 

@@ -31,7 +31,7 @@ value varmap jt s =    let v =  fromJust $ Map.lookup s varmap     in -    factorToList (fromJust $ posterior jt v) +    factorToList (fromJust $ posterior jt [v])   testWithRef varmap jt s l = assertBool s $ value varmap jt s ~=~ l testWithRefAndPrint varmap jt s l = do

@@ -7,13 +7,13 @@ -- The package version. See the Haskell package versioning policy -- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for -- standards guiding when and how versions should be incremented.-Version:             0.4.1+Version:             0.5  -- A short (one-line) description of the package.-Synopsis:            Inference with Discrete Bayesian Networks+Synopsis:            Bayesian Networks  -- A longer description of the package.-Description:  Algorithms for inference with Discrete Bayesian Networks.  +Description:  Algorithms for Bayesian Networks.    It is a very preliminary version. It has only been tested on very simple  examples where it worked. It should be considered as experimental and not used  in any production work.@@ -80,6 +80,12 @@     Bayes.InfluenceDiagram     Bayes.VariableElimination.Buckets     Bayes.Examples.Influence+    Bayes.Sampling+    Bayes.Examples.Sampling+    Bayes.EM+    Bayes.Examples.EMTest+    Bayes.Continuous+    Bayes.Examples.ContinuousSampling   other-modules:     Paths_hbayes     Bayes.ImportExport.HuginNet.Splitting@@ -100,11 +106,11 @@   -- Packages needed in order to build this package.   Build-depends:            base < 5,-    mtl == 2.0.1.0,-    containers == 0.4.2.1,-    array == 0.4.0.0,-    QuickCheck == 2.4.2,-    pretty == 1.1.1.0,+    mtl >= 2.0.1.0,+    containers >= 0.4.2.1,+    array >= 0.4.0.0,+    QuickCheck >= 2.4.2,+    pretty >= 1.1.1.0,     boxes,     vector,     random,@@ -116,4 +122,7 @@     test-framework-quickcheck2,     test-framework,     test-framework-hunit,-    HUnit+    HUnit,+    mwc-random >= 0.12,+    statistics >= 0.10.1,+    gamma >= 0.9