diff --git a/Data/Distribution.hs b/Data/Distribution.hs
--- a/Data/Distribution.hs
+++ b/Data/Distribution.hs
@@ -43,6 +43,10 @@
       -- | The "Data.Distribution.Sample" module provides ways to efficiently
       --   randomly sample values from distributions.
 
+    , module Data.Distribution.Monadic
+      -- | The "Data.Distribution.Monadic" module provides a monadic interface
+      --   to build distributions.
+
       -- * Other modules
 
       -- | Not all modules related to distributions are exported by default.
@@ -169,7 +173,9 @@
 
 import System.Random (mkStdGen)  -- For doctest.
 
+import Data.Distribution.Aggregator
 import Data.Distribution.Core
 import Data.Distribution.Measure
-import Data.Distribution.Aggregator
+import Data.Distribution.Monadic
 import Data.Distribution.Sample
+
diff --git a/Data/Distribution/Core.hs b/Data/Distribution/Core.hs
--- a/Data/Distribution/Core.hs
+++ b/Data/Distribution/Core.hs
@@ -35,15 +35,18 @@
       -- ** Transformation
     , select
     , assuming
+    , observing
       -- ** Combination
-    , combine
+    , combineWith
       -- ** Sequences
       -- *** Independant experiments
     , trials
     , times
+    , iid
       -- *** Dependant experiments
     , andThen
-    , on
+      -- ** Utilities
+    , isValid
     ) where
 
 import Control.Arrow (second)
@@ -134,15 +137,15 @@
     abs = select abs
     signum = select signum
     negate = select negate
-    d1 + d2 = d1 `andThen` (+) `on` d2
-    d1 - d2 = d1 `andThen` (-) `on` d2
-    d1 * d2 = d1 `andThen` (*) `on` d2
+    d1 + d2 = combineWith (+) d1 d2
+    d1 - d2 = combineWith (-) d1 d2
+    d1 * d2 = combineWith (*) d1 d2
 
 -- Binary operations on distributions are defined to
 -- be the binary operation on each pair of elements.
 instance (Ord a, Fractional a) => Fractional (Distribution a) where
     fromRational = always . fromRational
-    d1 / d2 = d1 `andThen` (/) `on` d2
+    d1 / d2 = combineWith (+) d1 d2
     recip = select recip
 
 -- Binary operations on distributions are defined to
@@ -152,8 +155,8 @@
     exp = select exp
     sqrt = select sqrt
     log = select log
-    d1 ** d2 = d1 `andThen` (**) `on` d2
-    d1 `logBase` d2 = d1 `andThen` logBase `on` d2
+    d1 ** d2 = combineWith (**) d1 d2
+    d1 `logBase` d2 = combineWith logBase d1 d2
     sin = select sin
     tan = select tan
     cos = select cos
@@ -169,7 +172,7 @@
 
 instance (Ord a, Monoid a) => Monoid (Distribution a) where
     mempty = always mempty
-    mappend d1 d2 = d1 `andThen` mappend `on` d2
+    mappend d1 d2 = combineWith mappend d1 d2
 
 -- | Converts the distribution to a list of increasing values whose probability
 --   is greater than @0@. To each value is associated its probability.
@@ -257,7 +260,7 @@
 --   >>> assuming (> 2) $ uniform [1 .. 6]
 --   fromList [(3,1 % 4),(4,1 % 4),(5,1 % 4),(6,1 % 4)]
 --
---   Note that the resulting distribution will be empty
+--   Note that the resulting distribution will be invalid
 --   if the predicate does not hold on any of the values.
 --
 --   >>> assuming (> 7) $ uniform [1 .. 6]
@@ -270,28 +273,48 @@
     total = sum $ Map.elems filtered
 
 
--- Combination
-
-
--- | Combines multiple weighted distributions into a single distribution.
+-- | Returns a new distribution using the Bayesian update rule.
 --
---   The probability of each element is the weighted sum of the element's
---   probability in every distribution.
+-- Using this example:
+-- https://en.wikipedia.org/wiki/Bayesian_inference#Probability_of_a_hypothesis
 --
---   >>> combine [(always 2, 1 / 3), (uniform [1..6], 2 / 3)]
---   fromList [(1,1 % 9),(2,4 % 9),(3,1 % 9),(4,1 % 9),(5,1 % 9),(6,1 % 9)]
+-- > data CookieBowl = Bowl1 | Bowl2 deriving (Eq,Ord)
+-- > data CookieType = Plain | ChocolateChip deriving (Eq,Ord)
+-- >
+-- > assumption :: Distribution CookieBowl
+-- > assumption = uniform [Bowl1,Bowl2]
+-- >
+-- > update :: Cookie -> Distribution CookieBowl -> Distribution CookieBowl
+-- > update c = observing f where
+-- >   f b = case b of
+-- >     -- Bowl #1 contains 10 chocolate chip cookies and 30 plain cookies
+-- >     Bowl1 -> fromList [(c == ChocolateChip,10),(c == Plain,30)]
+-- >     -- Bowl #2 contains 20 of each flavour of cookie
+-- >     Bowl2 -> fromList [(c == ChocolateChip,20),(c == Plain,20)]
 --
---   Note that the weights do not have to sum up to @1@. Distributions with
---   negative or null weight will be ignored.
-combine :: (Ord a, Real p) => [(Distribution a, p)] -> Distribution a
-combine dws = Distribution $ Map.unionsWith (+) $ zipWith go ds ps
+-- The "update" function in this example can be used to update the probability
+-- distribution of which bowl you have based on observing a random cookie inside
+-- the bowl.
+observing :: (a -> Distribution Bool) -> Distribution a -> Distribution a
+observing f (Distribution xs) = Distribution $ fmap adjust filtered
   where
-    (ds, ws) = unzip $ filter ((> 0) . snd) $ map (second toRational) dws
-    w = sum ws
-    ps = map (/ w) ws
-    go (Distribution xs) p = fmap (* p) xs
+    filtered = Map.filter (/= 0) $ Map.mapWithKey tweak xs
+    tweak x p = let
+      Distribution px = f x
+      pt = fromMaybe 0 $ Map.lookup True px
+      in pt * p
+    adjust x = x * (1 / total)
+    total = sum $ Map.elems filtered
 
+-- Combination
 
+combineWith :: (Ord b) => (a -> a -> b) -> Distribution a -> Distribution a -> Distribution b
+combineWith f (Distribution xs) (Distribution ys) = Distribution $ Map.unionsWith (+) $ do
+    (x, p) <- Map.toList xs
+    return $ Map.fromListWith (+) $ do
+        (y, q) <- Map.toList ys
+        return (f x y, p * q)
+
 -- Sequences
 
 
@@ -363,13 +386,29 @@
             let p' = 2 * p * q
             return (y + x, p')
 
+iid :: (Ord a) => (a -> a -> a) -> Int -> Distribution a -> Distribution a
+iid f n d
+    | n <= 0 = error "Called iid with a non-positive number of trials."
+    | otherwise = go n
+    where
+      go 1 = d
+      go m =
+        let (i, j) = quotRem m 2 
+            sub = go i
+            combined = combineWith f sub sub
+        in if j == 0
+          then combined
+          else combineWith f combined d
+
+
 -- | Computes for each value in the distribution a new distribution, and then
 --   combines those distributions, giving each the weight of the original value.
 --
 --   >>> uniform [1 .. 3] `andThen` (\ n -> uniform [1 .. n])
 --   fromList [(1,11 % 18),(2,5 % 18),(3,1 % 9)]
 --
---   See the 'on' function for a convenient way to chain distributions.
+--   See the 'Experiment' data type in the 'Data.Distribution.Monadic' module
+--   for a more "natural" monadic interface. 
 infixl 7 `andThen`
 andThen :: Ord b => Distribution a -> (a -> Distribution b) -> Distribution b
 andThen (Distribution xs) f = Distribution $
@@ -377,12 +416,10 @@
   where
     go (x, p) = fmap (* p) $ toMap $ f x
 
--- | Utility to partially apply a function on a distribution.
---   A use case for 'on' is to use it in conjunction with 'andThen'
---   to combine distributions.
+
+-- | Determines if a distribution is valid.
 --
---   >>> uniform [1 .. 3] `andThen` (+) `on` uniform [1 .. 2]
---   fromList [(2,1 % 6),(3,1 % 3),(4,1 % 3),(5,1 % 6)]
-infixl 8 `on`
-on :: Ord c => (a -> b -> c) -> Distribution b -> a -> Distribution c
-on f d x = select (f x) d
+--   A distribution is valid if and only if its domain is non-empty.
+--   Invalid distributions may arise from the use of 'assuming' for instance.
+isValid :: Distribution a -> Bool
+isValid (Distribution xs) = not $ Map.null xs
diff --git a/Data/Distribution/Monadic.hs b/Data/Distribution/Monadic.hs
new file mode 100644
--- /dev/null
+++ b/Data/Distribution/Monadic.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE GADTs #-}
+
+{- Copyright 2016 Romain Edelmann
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License. -}
+
+-- | This modules provides a monadic interface to build distributions.
+module Data.Distribution.Monadic 
+    ( Experiment, from, run ) where
+
+import Data.Distribution.Core
+
+-- | Monadic description of distributions.
+data Experiment a where
+  Return :: a -> Experiment a
+  Bind :: Experiment b -> (b -> Experiment a) -> Experiment a
+  Prim :: (Ord a) => Distribution a -> Experiment a
+
+instance Functor Experiment where
+  fmap f d = Bind d (\ x -> Return (f x))
+
+instance Applicative Experiment where
+  pure x = Return x
+  df <*> d = Bind df (\ f -> Bind d (\ x -> Return (f x)))
+
+instance Monad Experiment where
+  return x = Return x
+  d >>= f = Bind d f
+
+-- | Converts a concrete distribution into its
+--   monadic representation.
+from :: (Ord a) => Distribution a -> Experiment a
+from d = Prim d
+
+-- | Converts the monadic description of the distribution 
+--   to a concrete distribution.
+run :: (Ord a) => Experiment a -> Distribution a
+run (Return x) = always x
+run (Bind (Return x) f) = run (f x)
+run (Bind (Bind i g) f) = run (Bind i (\ x -> Bind (g x) f))
+run (Bind (Prim d) f) = d `andThen` \ x -> run (f x)
+run (Prim d) = d
diff --git a/Data/Distribution/Sample.hs b/Data/Distribution/Sample.hs
--- a/Data/Distribution/Sample.hs
+++ b/Data/Distribution/Sample.hs
@@ -176,8 +176,8 @@
 getSample :: MonadRandom m => Generator a -> m a
 getSample g = do
     let n = capacity g
-    u <- getRandom
     j <- getRandomR (0, n - 1)
+    u <- getRandom
     let i = if u < probabilities g ! j
                 then j
                 else indexes g ! j
diff --git a/distribution.cabal b/distribution.cabal
--- a/distribution.cabal
+++ b/distribution.cabal
@@ -8,7 +8,7 @@
 -- PVP summary:      +-+------- breaking API changes
 --                   | | +----- non-breaking API additions
 --                   | | | +--- code changes with no API change
-version:             1.0.1.0
+version:             1.1.1.0
 
 -- A short (one-line) description of the package.
 synopsis:             Finite discrete probability distributions.
@@ -51,8 +51,10 @@
                        Data.Distribution.Domain.Coin,
                        Data.Distribution.Domain.Dice,
                        Data.Distribution.Measure,
+                       Data.Distribution.Monadic,
                        Data.Distribution.Sample
 
+
   -- Modules included in this library but not exported.
   -- other-modules:
 
@@ -60,6 +62,6 @@
   build-depends:       array >=0.4,
                        base >=4.5 && <5,
                        containers ==0.5.*,
-                       MonadRandom ==0.4.*,
+                       MonadRandom >=0.4,
                        random ==1.1.*