diff --git a/Bayes.hs b/Bayes.hs
--- a/Bayes.hs
+++ b/Bayes.hs
@@ -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
diff --git a/Bayes/BayesianNetwork.hs b/Bayes/BayesianNetwork.hs
--- a/Bayes/BayesianNetwork.hs
+++ b/Bayes/BayesianNetwork.hs
@@ -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 .==.
diff --git a/Bayes/Continuous.hs b/Bayes/Continuous.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Continuous.hs
@@ -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
diff --git a/Bayes/EM.hs b/Bayes/EM.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/EM.hs
@@ -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)
diff --git a/Bayes/Examples/ContinuousSampling.hs b/Bayes/Examples/ContinuousSampling.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Examples/ContinuousSampling.hs
@@ -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
+
+
diff --git a/Bayes/Examples/EMTest.hs b/Bayes/Examples/EMTest.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Examples/EMTest.hs
@@ -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)
+ 
diff --git a/Bayes/Examples/Sampling.hs b/Bayes/Examples/Sampling.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Examples/Sampling.hs
@@ -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
diff --git a/Bayes/Examples/Tutorial.hs b/Bayes/Examples/Tutorial.hs
--- a/Bayes/Examples/Tutorial.hs
+++ b/Bayes/Examples/Tutorial.hs
@@ -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 ""
 
diff --git a/Bayes/Factor.hs b/Bayes/Factor.hs
--- a/Bayes/Factor.hs
+++ b/Bayes/Factor.hs
@@ -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
  , (=:)
diff --git a/Bayes/Factor/CPT.hs b/Bayes/Factor/CPT.hs
--- a/Bayes/Factor/CPT.hs
+++ b/Bayes/Factor/CPT.hs
@@ -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
diff --git a/Bayes/Factor/MaxCPT.hs b/Bayes/Factor/MaxCPT.hs
--- a/Bayes/Factor/MaxCPT.hs
+++ b/Bayes/Factor/MaxCPT.hs
@@ -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
diff --git a/Bayes/Factor/PrivateCPT.hs b/Bayes/Factor/PrivateCPT.hs
--- a/Bayes/Factor/PrivateCPT.hs
+++ b/Bayes/Factor/PrivateCPT.hs
@@ -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
diff --git a/Bayes/FactorElimination.hs b/Bayes/FactorElimination.hs
--- a/Bayes/FactorElimination.hs
+++ b/Bayes/FactorElimination.hs
@@ -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)
 
diff --git a/Bayes/InfluenceDiagram.hs b/Bayes/InfluenceDiagram.hs
--- a/Bayes/InfluenceDiagram.hs
+++ b/Bayes/InfluenceDiagram.hs
@@ -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
diff --git a/Bayes/Network.hs b/Bayes/Network.hs
--- a/Bayes/Network.hs
+++ b/Bayes/Network.hs
@@ -22,6 +22,9 @@
 	, runNetwork
 	, execNetwork
 	, evalNetwork
+  , runGraph
+  , execGraph
+  , evalGraph
 	, getCpt
 	) where 
 
diff --git a/Bayes/PrivateTypes.hs b/Bayes/PrivateTypes.hs
--- a/Bayes/PrivateTypes.hs
+++ b/Bayes/PrivateTypes.hs
@@ -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
diff --git a/Bayes/Sampling.hs b/Bayes/Sampling.hs
new file mode 100644
--- /dev/null
+++ b/Bayes/Sampling.hs
@@ -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
diff --git a/Bayes/Test/CompareEliminations.hs b/Bayes/Test/CompareEliminations.hs
--- a/Bayes/Test/CompareEliminations.hs
+++ b/Bayes/Test/CompareEliminations.hs
@@ -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])
 
diff --git a/Bayes/Test/ReferencePatterns.hs b/Bayes/Test/ReferencePatterns.hs
--- a/Bayes/Test/ReferencePatterns.hs
+++ b/Bayes/Test/ReferencePatterns.hs
@@ -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
diff --git a/hbayes.cabal b/hbayes.cabal
--- a/hbayes.cabal
+++ b/hbayes.cabal
@@ -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