diff --git a/Bayes.hs b/Bayes.hs
--- a/Bayes.hs
+++ b/Bayes.hs
@@ -777,17 +777,17 @@
 
 instance (Show b, Show e) => Show (UndirectedSG e b)where
   show g@(SP em vm nm) = execWriter $ do
-  tell "graph dot {\n"
-  mapM_ (addVertexToUndirectedGraphviz nm) $ IM.toList vm
-  tell "\n"
-  mapM_ (addEdgeToGraphviz) $ M.toList em
-  tell "}\n"
-   where
-     addEdgeToGraphviz (e@(Edge (Vertex vs) (Vertex ve)),l) = do
-       tell $ show vs 
-       tell " -- "
-       tell $ show ve
-       tell "\n"
+    tell "graph dot {\n"
+    mapM_ (addVertexToUndirectedGraphviz nm) $ IM.toList vm
+    tell "\n"
+    mapM_ (addEdgeToGraphviz) $ M.toList em
+    tell "}\n"
+     where
+       addEdgeToGraphviz (e@(Edge (Vertex vs) (Vertex ve)),l) = do
+         tell $ show vs 
+         tell " -- "
+         tell $ show ve
+         tell "\n"
 
 
 displayFactors :: (NeighborhoodStructure n, Show f, Factor f, Graph (SimpleGraph n)) => SimpleGraph n a f -> String
@@ -821,7 +821,7 @@
 -- | Graph monad.
 -- The monad used to simplify the description of a new graph
 -- g is the graph type. e the edge type. f the node type (generally a 'Factor')
-newtype GraphMonad g e f a = GM {runGraphMonad :: State (GMState g e f) a} deriving(Monad, MonadState (GMState g e f))
+newtype GraphMonad g e f a = GM {runGraphMonad :: State (GMState g e f) a} deriving(Functor,Applicative,Monad, MonadState (GMState g e f))
 
 -- | Get a named vertex from the graph monad
 getVertex :: Graph g => String -> GraphMonad g e f (Maybe Vertex)
diff --git a/Bayes/Continuous.hs b/Bayes/Continuous.hs
--- a/Bayes/Continuous.hs
+++ b/Bayes/Continuous.hs
@@ -23,6 +23,7 @@
     , beta
     , beta'
     , exponential
+    , gammaD
     , execCN
     , runCN 
     , evalCN
@@ -73,6 +74,26 @@
    (RN a) / (RN b) = RN $ liftA2 (/) a b
    fromRational a = RN (return (fromRational a))
 
+instance Floating RN where 
+  pi = RN (return pi)
+  exp (RN b) = RN $ liftA exp b
+  sqrt (RN b) = RN $ liftA sqrt b 
+  log (RN b) = RN $ liftA log b
+  sin (RN b) = RN $ liftA sin b
+  cos (RN b) = RN $ liftA cos b
+  tan (RN b) = RN $ liftA tan b
+  asin (RN b) = RN $ liftA asin b 
+  atan (RN b) = RN $ liftA atan b 
+  acos (RN b) = RN $ liftA acos b 
+  sinh (RN b) = RN $ liftA sinh b 
+  tanh (RN b) = RN $ liftA tanh b 
+  cosh (RN b) = RN $ liftA cosh b 
+  asinh (RN b) = RN $ liftA asinh b 
+  atanh (RN b) = RN $ liftA atanh b 
+  acosh (RN b) = RN $ liftA acosh b 
+  (**) (RN xb) (RN yb) = RN $ liftA2 (**) xb yb
+  (logBase) (RN xb) (RN yb) = RN $ liftA2 logBase xb yb
+
 -- | 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.
@@ -107,7 +128,28 @@
     (DN da na) / (DN db nb) = DN (da `mappend`db) (na / nb)
     fromRational i = DN [] (fromRational i)
 
+instance Floating DN where 
+  pi = DN [] pi
+  exp (DN a b) = DN a (exp b) 
+  sqrt (DN a b) = DN a (sqrt b) 
+  log (DN a b) = DN a (log b) 
+  sin (DN a b) = DN a (sin b) 
+  cos (DN a b) = DN a (cos b) 
+  tan (DN a b) = DN a (tan b) 
+  asin (DN a b) = DN a (asin b) 
+  atan (DN a b) = DN a (atan b) 
+  acos (DN a b) = DN a (acos b) 
+  sinh (DN a b) = DN a (sinh b) 
+  tanh (DN a b) = DN a (tanh b) 
+  cosh (DN a b) = DN a (cosh b) 
+  asinh (DN a b) = DN a (asinh b) 
+  atanh (DN a b) = DN a (atanh b) 
+  acosh (DN a b) = DN a (acosh b) 
+  (**) (DN xa xb) (DN ya yb) = DN (xa `mappend` ya) (xb ** yb)
+  (logBase) (DN xa xb) (DN ya yb) = DN (xa `mappend` ya) (xb ** yb)
 
+
+
 {-
 
 Graph creation
@@ -205,6 +247,35 @@
     setBayesianNode v result 
     return (node v)
 
+-- | Gamma distribution 
+gammaD :: VariableName s 
+       => s 
+       -> DN -- ^ r
+       -> DN -- ^ lambda
+       -> CNMonad DN 
+gammaD s rm lambdam = do 
+    v <- mkVariable s
+    let la = dependencies rm 
+        lb = dependencies lambdam 
+    let l = la ++ lb
+    cpt v l
+    let bound = do 
+            lambda <- value lambdam 
+            r <- value rm
+            return . Unbounded $ r/lambda/lambda
+        result = D v (Distri bound $ \inst -> do 
+                    lambda <- value lambdam
+                    r <- value rm
+                    let x = instantiationValue inst 
+                    if x < 0 
+                      then 
+                        return 0.0 
+                      else 
+                        return $ lambda**r * x**(r-1)*exp(-lambda*x)/gamma r
+                    )
+    setBayesianNode v result 
+    return (node v)
+
 -- | Beta distribution 
 beta :: VariableName s 
      => s 
@@ -270,10 +341,10 @@
         lb = dependencies nb 
     let l = la ++ lb
     cpt v l
-    let bound = do 
-        a <- value na 
-        b <- value nb 
-        return (BoundedSupport (a,b) )
+    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
@@ -296,8 +367,8 @@
     let l = la ++ lb
     cpt v l
     let bound = do 
-        s <- value std 
-        return (Unbounded s )
+         s <- value std 
+         return (Unbounded s )
     let result = D v (Distri bound $ \inst -> do 
                     m <- value mean
                     s <- value std
diff --git a/Bayes/Examples/ContinuousSampling.hs b/Bayes/Examples/ContinuousSampling.hs
--- a/Bayes/Examples/ContinuousSampling.hs
+++ b/Bayes/Examples/ContinuousSampling.hs
@@ -81,18 +81,75 @@
 
 We see in the histogram that the estimated value is around 1.5.
 
+The example 'complexsamples' will create three files alpha.txt, beta.txt and tau.txt. It is corresponding to the following
+bugs model
+
+@
+model {
+    for (i in 1:n) {
+          mu[i] <- alpha + beta*i/n;
+          y[i]   ~ dnorm(mu[i],tau);
+    }
+    alpha    ~ dnorm(0.0,1.0E-4);
+    beta     ~ dnorm(0.0,1.0E-4);
+    tau      ~ dgamma(1.0E-3,1.0E-3);
+    sigma   <- 1.0/sqrt(tau);
+}
+@
+
+with alpha = 0, beta = 0 and tau = 1
+
+The Haskell code for this model is
+
+@
+complexMeasures = 100 
+--
+complex = 'runCN' $ do 
+  let n = complexMeasures
+  alpha <- 'normal' \"alpha\" 0.0 1e-4 
+  beta <- 'normal' \"beta\" 0.0 1e-4 
+  tau <- 'gammaD' \"tau\" 1e-3 1e-3 
+  let sigma = 1.0 / sqrt(tau)
+      sample i = do 
+        let mu = alpha + beta * fromIntegral i / fromIntegral n
+        y <- 'normal' () mu tau 
+        return y
+  l <- mapM sample [0..n-1]
+  return (alpha:beta:tau:l)
+@
+
+And the generation of the samples is done with
+
+@
+complexsamples = do 
+  let n = complexMeasures
+      ((alpha:beta:tau:obs),complexg) = complex 
+      alphat = 0.0 
+      betat = 0.0 
+      taut = 1.0 
+      aMeasure g i = do 
+        let mu = alphat + betat * fromIntegral i / fromIntegral n
+        MWC.normal mu taut g
+  g <- create
+  measurements <- mapM (aMeasure g) [0..n-1] 
+  let evidence = zipWith ('=:') obs measurements
+@
+
 -}
 module Bayes.Examples.ContinuousSampling(
-	  nbSensors 
-	, sensor 
-	, test 
-	, debugcn
-	) where
+      nbSensors 
+    , sensor 
+    , test 
+    , debugcn
+  , writesamples
+  , complexsamples
+    ) 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)
+import System.IO(withFile,IOMode(..),hPrint)
 
 nbSensors = 10 
 
@@ -105,6 +162,41 @@
   sensors <- sequence (replicate nbSensors (sensor a))
   return (a:sensors)
 
+complexMeasures = 100 
+
+complex = runCN $ do 
+  let n = complexMeasures
+  alpha <- normal "alpha" 0.0 1e-4 
+  beta <- normal "beta" 0.0 1e-4 
+  tau <- gammaD "tau" 1e-3 1e-3 
+  let sigma = 1.0 / sqrt(tau)
+      sample i = do 
+        let mu = alpha + beta * fromIntegral i / fromIntegral n
+        y <- normal () mu tau 
+        return y
+  l <- mapM sample [0..n-1]
+  return (alpha:beta:tau:l)
+
+complexsamples = do 
+  let n = complexMeasures
+      ((alpha:beta:tau:obs),complexg) = complex 
+      alphat = 0.0 
+      betat = 0.0 
+      taut = 1.0 
+      aMeasure g i = do 
+        let mu = alphat + betat * fromIntegral i / fromIntegral n
+        MWC.normal mu taut g
+  g <- create
+  measurements <- mapM (aMeasure g) [0..n-1] 
+  let evidence = zipWith (=:) obs measurements
+  n <- runSampling 10000 200 (continuousMCMCSampler complexg evidence)
+  let samples s a = do 
+       let v = map (\g -> instantiationValue . fromJust . vertexValue g $ (vertex a)) n
+       writeList s v
+  samples "alpha.txt" alpha 
+  samples "beta.txt" beta
+  samples "tau.txt" tau
+
 debugcn = do 
     let ((a:sensors), testG) = test
     g <- create 
@@ -115,4 +207,16 @@
         h = histogram 6 samples 
     print h
 
+writeList :: Show a => FilePath -> [a] -> IO ()
+writeList s l = do 
+  withFile s WriteMode $ \h -> do 
+    mapM_ (hPrint h) l
 
+writesamples s = 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
+    writeList s samples
diff --git a/Bayes/ImportExport.hs b/Bayes/ImportExport.hs
--- a/Bayes/ImportExport.hs
+++ b/Bayes/ImportExport.hs
@@ -4,13 +4,13 @@
 
 -}
 module Bayes.ImportExport (
-	-- * Networks
-	  writeNetworkToFile
-	, readNetworkFromFile
-	-- * Junction Tree 
-	, writeVariableMapAndJunctionTreeToFile
-	, readVariableMapAndJunctionTreeToFile
-	) where 
+    -- * Networks
+      writeNetworkToFile
+    , readNetworkFromFile
+    -- * Junction Tree 
+    , writeVariableMapAndJunctionTreeToFile
+    , readVariableMapAndJunctionTreeToFile
+    ) where 
 
 import Data.Binary
 import Bayes
@@ -46,130 +46,130 @@
 readVariableMapAndJunctionTreeToFile f = decodeFile f
 
 instance Binary Cluster where 
-	put (Cluster s) = put s 
-	get = get >>= return . Cluster
+    put (Cluster s) = put s 
+    get = get >>= return . Cluster
 
 instance (Ord c, Binary c, Binary f) => Binary (JTree c f) where 
-	put (JTree r ls cm pm spm scm nvm svm sck sclm) = do 
-		put r 
-		put ls 
-		put cm 
-		put pm 
-		put spm 
-		put scm 
-		put nvm 
-		put svm 
-		put sck 
-		put sclm 
-	get = do 
-		r <- get
-		ls <- get 
-		cm <- get 
-		pm <- get 
-		spm <- get 
-		scm <- get 
-		nvm <- get 
-		svm <- get 
-		sck <- get 
-		sclm <- get
-		return $ JTree r ls cm pm spm scm nvm svm sck sclm
+    put (JTree r ls cm pm spm scm nvm svm sck sclm) = do 
+        put r 
+        put ls 
+        put cm 
+        put pm 
+        put spm 
+        put scm 
+        put nvm 
+        put svm 
+        put sck 
+        put sclm 
+    get = do 
+        r <- get
+        ls <- get 
+        cm <- get 
+        pm <- get 
+        spm <- get 
+        scm <- get 
+        nvm <- get 
+        svm <- get 
+        sck <- get 
+        sclm <- get
+        return $ JTree r ls cm pm spm scm nvm svm sck sclm
 
 instance Binary a => Binary (NodeValue a) where 
-	put (NodeValue v f e) = do 
-		put v 
-		put f 
-		put e 
-	get = do 
-		v <- get 
-		f <- get 
-		e <- get
-		return $ NodeValue v f e
+    put (NodeValue v f e) = do 
+        put v 
+        put f 
+        put e 
+    get = do 
+        v <- get 
+        f <- get 
+        e <- get
+        return $ NodeValue v f e
 
 instance Binary a => Binary (SeparatorValue a) where
-	put (EmptySeparator) = do 
-		putWord8 0 
-	put (SeparatorValue a b) = do 
-		putWord8 1 
-		put a 
-		put b 
-	get = do 
-		tag <- getWord8 
-		case tag of 
-			0 -> return EmptySeparator
-			_ -> do 
-				a <- get 
-				b <- get 
-				return $ SeparatorValue a b 
+    put (EmptySeparator) = do 
+        putWord8 0 
+    put (SeparatorValue a b) = do 
+        putWord8 1 
+        put a 
+        put b 
+    get = do 
+        tag <- getWord8 
+        case tag of 
+            0 -> return EmptySeparator
+            _ -> do 
+                a <- get 
+                b <- get 
+                return $ SeparatorValue a b 
 
 instance Binary (V.Vector Double) where 
-	put = put . V.toList 
-	get = get >>= return . V.fromList 
+    put = put . V.toList 
+    get = get >>= return . V.fromList 
 
 instance Binary (NV.Vector Double) where 
-	put = put . NV.toList 
-	get = get >>= return . NV.fromList 
+    put = put . NV.toList 
+    get = get >>= return . NV.fromList 
 
 instance Binary DV where 
-	put (DV v i) = do 
-		put v 
-		put i 
-	get = do 
-		v <- get 
-		i <- get 
-		return $ DV v i 
+    put (DV v i) = do 
+        put v 
+        put i 
+    get = do 
+        v <- get 
+        i <- get 
+        return $ DV v i 
 
 instance Binary (v Double) => Binary (PrivateCPT v Double) where 
     put (Table d m v) = do 
-    	putWord8 0 
-    	put d 
-    	put m 
-    	put v 
+        putWord8 0 
+        put d 
+        put m 
+        put v 
     put (Scalar v) = do 
-    	putWord8 1 
-    	put v 
+        putWord8 1 
+        put v 
     get = do 
-    	tag <- getWord8 
-    	case tag of 
-    		0 -> do 
-    			d <- get 
-    			m <- get 
-    			v <- get 
-    			return $ Table d m v 
-    		_ -> get >>= return . Scalar 
+        tag <- getWord8 
+        case tag of 
+            0 -> do 
+                d <- get 
+                m <- get 
+                v <- get 
+                return $ Table d m v 
+            _ -> get >>= return . Scalar 
 
 instance Binary Vertex where 
   put (Vertex v) = put v 
   get = get >>= return . Vertex 
 
 instance Binary Edge where 
-	put (Edge va vb) = do 
-		put va 
-		put vb 
-	get = do 
-		va <- get
-		vb <- get 
-		return $ Edge va vb
+    put (Edge va vb) = do 
+        put va 
+        put vb 
+    get = do 
+        va <- get
+        vb <- get 
+        return $ Edge va vb
 
 instance (Binary l, Binary e, Binary v) => Binary (SimpleGraph l e v) where 
-	put (SP e v n) = do 
-		put e 
-		put v 
-		put n 
-	get = do 
-		e <- get 
-		v <- get 
-		n <- get 
-		return $ SP e v n
+    put (SP e v n) = do 
+        put e 
+        put v 
+        put n 
+    get = do 
+        e <- get 
+        v <- get 
+        n <- get 
+        return $ SP e v n
 
 instance Binary DE where 
     put (DE a b) = do 
-    	put a 
-    	put b 
+        put a 
+        put b 
     get = do 
-    	a <- get 
-    	b <- get 
-    	return $ DE a b 
+        a <- get 
+        b <- get 
+        return $ DE a b 
 
 instance Binary UE where 
-	put (UE a) = put a 
-	get = get >>= return . UE
+    put (UE a) = put a 
+    get = get >>= return . UE
diff --git a/Bayes/Network.hs b/Bayes/Network.hs
--- a/Bayes/Network.hs
+++ b/Bayes/Network.hs
@@ -3,30 +3,30 @@
 
 -}
 module Bayes.Network(
-	-- * Types 
-	  MaybeNode(..)
-	, NetworkMonad(..)
-	-- * Functions
-	, factorVariable
-	, (<--)
-	, getBayesianNode 
-	, setBayesianNode
-	, initializeNodeWithValue
-	, setVariableBoundWithSize
-	, setVariableBound
-	, addVariableIfNotFound
-	, unamedVariable
-	, variable
-	, variableWithSize
-	, unNamedVariableWithSize
-	, runNetwork
-	, execNetwork
-	, evalNetwork
+  -- * Types 
+    MaybeNode(..)
+  , NetworkMonad(..)
+  -- * Functions
+  , factorVariable
+  , (<--)
+  , getBayesianNode 
+  , setBayesianNode
+  , initializeNodeWithValue
+  , setVariableBoundWithSize
+  , setVariableBound
+  , addVariableIfNotFound
+  , unamedVariable
+  , variable
+  , variableWithSize
+  , unNamedVariableWithSize
+  , runNetwork
+  , execNetwork
+  , evalNetwork
   , runGraph
   , execGraph
   , evalGraph
-	, getCpt
-	) where 
+  , getCpt
+  ) where 
 
 import Bayes.PrivateTypes
 import Bayes 
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c)2012, alpheccar
+Copyright (c)2012-2016, alpheccar
 
 All rights reserved.
 
diff --git a/ModuleTest.hs b/ModuleTest.hs
new file mode 100644
--- /dev/null
+++ b/ModuleTest.hs
@@ -0,0 +1,3 @@
+import Bayes.Test
+
+main = runTests
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,6 @@
+[![Hackage](https://img.shields.io/hackage/v/hbayes.svg)](https://hackage.haskell.org/package/hbayes)
+
+hbayes
+======
+
+WARNING : This package is not yet using log arithmetic !
diff --git a/hbayes.cabal b/hbayes.cabal
--- a/hbayes.cabal
+++ b/hbayes.cabal
@@ -7,7 +7,7 @@
 -- 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.5
+Version:             0.5.2
 
 -- A short (one-line) description of the package.
 Synopsis:            Bayesian Networks
@@ -35,7 +35,7 @@
 Maintainer:          misc@alpheccar.org
 
 -- A copyright notice.
-Copyright: Copyright (c) 2012, alpheccar 
+Copyright: Copyright (c) 2012-2016, alpheccar 
 
 Stability: experimental      
 
@@ -45,8 +45,11 @@
 tested-with: GHC==7.4.1 
 
 -- Constraint on the version of Cabal needed to build this package.
-Cabal-version:       >=1.8
+Cabal-version:       >=1.10
 
+extra-source-files:
+  README.md
+
 data-files: cancer.net
 
 Flag local {
@@ -59,6 +62,42 @@
   Default:     False
 }
 
+source-repository head
+  type:     git
+  location: https://github.com/alpheccar/hbayes.git
+
+Test-Suite hbayes-Tests
+  Type:              exitcode-stdio-1.0
+  Main-is:           ModuleTest.hs
+  hs-source-dirs:    .
+  Build-depends:     base >= 4, 
+                     hbayes,
+                     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,
+                     boxes,
+                     vector,
+                     random,
+                     split,
+                     parsec,
+                     filepath,
+                     directory,
+                     binary >= 0.5,
+                     test-framework-quickcheck2,
+                     test-framework,
+                     test-framework-hunit,
+                     HUnit,
+                     mwc-random >= 0.12,
+                     statistics >= 0.10.1,
+                     gamma >= 0.9
+  Default-language:  Haskell2010
+  default-extensions: CPP
+  if flag(local)
+    cpp-options: -DLOCAL
+
 Library
   -- Modules exported by the library.
   Exposed-modules:
@@ -95,8 +134,9 @@
     Bayes.Factor.PrivateCPT
     Bayes.Network
 
+  Default-language:  Haskell2010
   GHC-Options: -funbox-strict-fields
-  Extensions: CPP
+  default-extensions: CPP
   if flag(local)
     cpp-options: -DLOCAL
   if flag(profile)