diff --git a/.travis.yml b/.travis.yml
index b7b23e694..25750bac9 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -22,6 +22,7 @@ install:
   - pip install tensorflow
   - pip install opencv-python
 
+
 script:
   - py.test --cov=./
   - python -m doctest -v *.py
diff --git a/DeepNeuralNet4e.py b/DeepNeuralNet4e.py
index a353df95c..b68192ba8 100644
--- a/DeepNeuralNet4e.py
+++ b/DeepNeuralNet4e.py
@@ -1,6 +1,7 @@
 import math
 import statistics
-from utils4e import sigmoid, dotproduct, softmax1D, conv1D, GaussianKernel, element_wise_product, \
+
+from utils4e import sigmoid, dotproduct, softmax1D, conv1D, gaussian_kernel_2d, GaussianKernel, element_wise_product, \
     vector_add, random_weights, scalar_vector_product, matrix_multiplication, map_vector
 import random
 
diff --git a/probability4e.py b/probability4e.py
new file mode 100644
index 000000000..94429f2dd
--- /dev/null
+++ b/probability4e.py
@@ -0,0 +1,758 @@
+"""Probability models.
+"""
+
+from utils import product, argmax, isclose, probability
+from logic import extend
+from math import sqrt, pi, exp
+import copy
+import random
+from collections import defaultdict
+from functools import reduce
+
+# ______________________________________________________________________________
+# Chapter 12 Qualifying Uncertainty
+# 12.1 Acting Under Uncertainty
+
+
+def DTAgentProgram(belief_state):
+    """A decision-theoretic agent. [Figure 12.1]"""
+    def program(percept):
+        belief_state.observe(program.action, percept)
+        program.action = argmax(belief_state.actions(),
+                                key=belief_state.expected_outcome_utility)
+        return program.action
+    program.action = None
+    return program
+
+# ______________________________________________________________________________
+# 12.2 Basic Probability Notation
+
+
+class ProbDist:
+    """A discrete probability distribution. You name the random variable
+    in the constructor, then assign and query probability of values.
+    >>> P = ProbDist('Flip'); P['H'], P['T'] = 0.25, 0.75; P['H']
+    0.25
+    >>> P = ProbDist('X', {'lo': 125, 'med': 375, 'hi': 500})
+    >>> P['lo'], P['med'], P['hi']
+    (0.125, 0.375, 0.5)
+    """
+
+    def __init__(self, varname='?', freqs=None):
+        """If freqs is given, it is a dictionary of values - frequency pairs,
+        then ProbDist is normalized."""
+        self.prob = {}
+        self.varname = varname
+        self.values = []
+        if freqs:
+            for (v, p) in freqs.items():
+                self[v] = p
+            self.normalize()
+
+    def __getitem__(self, val):
+        """Given a value, return P(value)."""
+        try:
+            return self.prob[val]
+        except KeyError:
+            return 0
+
+    def __setitem__(self, val, p):
+        """Set P(val) = p."""
+        if val not in self.values:
+            self.values.append(val)
+        self.prob[val] = p
+
+    def normalize(self):
+        """Make sure the probabilities of all values sum to 1.
+        Returns the normalized distribution.
+        Raises a ZeroDivisionError if the sum of the values is 0."""
+        total = sum(self.prob.values())
+        if not isclose(total, 1.0):
+            for val in self.prob:
+                self.prob[val] /= total
+        return self
+
+    def show_approx(self, numfmt='{:.3g}'):
+        """Show the probabilities rounded and sorted by key, for the
+        sake of portable doctests."""
+        return ', '.join([('{}: ' + numfmt).format(v, p)
+                          for (v, p) in sorted(self.prob.items())])
+
+    def __repr__(self):
+        return "P({})".format(self.varname)
+
+# ______________________________________________________________________________
+# 12.3 Inference Using Full Joint Distributions
+
+
+class JointProbDist(ProbDist):
+    """A discrete probability distribute over a set of variables.
+    >>> P = JointProbDist(['X', 'Y']); P[1, 1] = 0.25
+    >>> P[1, 1]
+    0.25
+    >>> P[dict(X=0, Y=1)] = 0.5
+    >>> P[dict(X=0, Y=1)]
+    0.5"""
+
+    def __init__(self, variables):
+        self.prob = {}
+        self.variables = variables
+        self.vals = defaultdict(list)
+
+    def __getitem__(self, values):
+        """Given a tuple or dict of values, return P(values)."""
+        values = event_values(values, self.variables)
+        return ProbDist.__getitem__(self, values)
+
+    def __setitem__(self, values, p):
+        """Set P(values) = p.  Values can be a tuple or a dict; it must
+        have a value for each of the variables in the joint. Also keep track
+        of the values we have seen so far for each variable."""
+        values = event_values(values, self.variables)
+        self.prob[values] = p
+        for var, val in zip(self.variables, values):
+            if val not in self.vals[var]:
+                self.vals[var].append(val)
+
+    def values(self, var):
+        """Return the set of possible values for a variable."""
+        return self.vals[var]
+
+    def __repr__(self):
+        return "P({})".format(self.variables)
+
+
+def event_values(event, variables):
+    """Return a tuple of the values of variables in event.
+    >>> event_values ({'A': 10, 'B': 9, 'C': 8}, ['C', 'A'])
+    (8, 10)
+    >>> event_values ((1, 2), ['C', 'A'])
+    (1, 2)
+    """
+    if isinstance(event, tuple) and len(event) == len(variables):
+        return event
+    else:
+        return tuple([event[var] for var in variables])
+
+
+def enumerate_joint_ask(X, e, P):
+    """Return a probability distribution over the values of the variable X,
+    given the {var:val} observations e, in the JointProbDist P. [Section 12.3]
+    >>> P = JointProbDist(['X', 'Y'])
+    >>> P[0,0] = 0.25; P[0,1] = 0.5; P[1,1] = P[2,1] = 0.125
+    >>> enumerate_joint_ask('X', dict(Y=1), P).show_approx()
+    '0: 0.667, 1: 0.167, 2: 0.167'
+    """
+    assert X not in e, "Query variable must be distinct from evidence"
+    Q = ProbDist(X)  # probability distribution for X, initially empty
+    Y = [v for v in P.variables if v != X and v not in e]  # hidden variables.
+    for xi in P.values(X):
+        Q[xi] = enumerate_joint(Y, extend(e, X, xi), P)
+    return Q.normalize()
+
+
+def enumerate_joint(variables, e, P):
+    """Return the sum of those entries in P consistent with e,
+    provided variables is P's remaining variables (the ones not in e)."""
+    if not variables:
+        return P[e]
+    Y, rest = variables[0], variables[1:]
+    return sum([enumerate_joint(rest, extend(e, Y, y), P)
+                for y in P.values(Y)])
+
+# ______________________________________________________________________________
+# 12.4 Independence
+
+
+def is_independent(variables, P):
+    """
+    Return whether a list of variables are independent given their distribution P
+    P is an instance of JoinProbDist
+    >>> P = JointProbDist(['X', 'Y'])
+    >>> P[0,0] = 0.25; P[0,1] = 0.5; P[1,1] = P[1,0] = 0.125
+    >>> is_independent(['X', 'Y'], P)
+    False
+    """
+    for var in variables:
+        event_vars = variables[:]
+        event_vars.remove(var)
+        event = {}
+        distribution = enumerate_joint_ask(var, event, P)
+        events = gen_possible_events(event_vars, P)
+        for e in events:
+            conditional_distr = enumerate_joint_ask(var, e, P)
+            if conditional_distr.prob != distribution.prob:
+                return False
+    return True
+
+
+def gen_possible_events(vars, P):
+    """Generate all possible events of a collection of vars according to distribution of P"""
+    events = []
+
+    def backtrack(vars, P, temp):
+        if not vars:
+            events.append(temp)
+            return
+        var = vars[0]
+        for val in P.values(var):
+            temp[var] = val
+            backtrack([v for v in vars if v != var], P, copy.copy(temp))
+    backtrack(vars, P, {})
+    return events
+
+# ______________________________________________________________________________
+# Chapter 13 Probabilistic Reasoning
+# 13.1 Representing Knowledge in an Uncertain Domain
+
+
+class BayesNet:
+    """Bayesian network containing only boolean-variable nodes."""
+
+    def __init__(self, node_specs=None):
+        """
+        Nodes must be ordered with parents before children.
+        :param node_specs: an nested iterable object, each element contains (variable name, parents name, cpt)
+                           for each node
+        """
+
+        self.nodes = []
+        self.variables = []
+        node_specs = node_specs or []
+        for node_spec in node_specs:
+            self.add(node_spec)
+
+    def add(self, node_spec):
+        """
+        Add a node to the net. Its parents must already be in the
+        net, and its variable must not.
+        Initialize Bayes nodes by detecting the length of input node specs
+        """
+        if len(node_spec)>=5:
+            node = ContinuousBayesNode(*node_spec)
+        else:
+            node = BayesNode(*node_spec)
+        assert node.variable not in self.variables
+        assert all((parent in self.variables) for parent in node.parents)
+        self.nodes.append(node)
+        self.variables.append(node.variable)
+        for parent in node.parents:
+            self.variable_node(parent).children.append(node)
+
+    def variable_node(self, var):
+        """
+        Return the node for the variable named var.
+        >>> burglary.variable_node('Burglary').variable
+        'Burglary'
+        """
+        for n in self.nodes:
+            if n.variable == var:
+                return n
+        raise Exception("No such variable: {}".format(var))
+
+    def variable_values(self, var):
+        """Return the domain of var."""
+        return [True, False]
+
+    def __repr__(self):
+        return 'BayesNet({0!r})'.format(self.nodes)
+
+
+class BayesNode:
+    """
+    A conditional probability distribution for a boolean variable,
+    P(X | parents). Part of a BayesNet.
+    """
+
+    def __init__(self, X, parents, cpt):
+        """
+        :param X: variable name,
+        :param parents: a sequence of variable names or a space-separated string. Representing the names of parent nodes.
+        :param cpt: the conditional probability table, takes one of these forms:
+
+        * A number, the unconditional probability P(X=true). You can
+          use this form when there are no parents.
+
+        * A dict {v: p, ...}, the conditional probability distribution
+          P(X=true | parent=v) = p. When there's just one parent.
+
+        * A dict {(v1, v2, ...): p, ...}, the distribution P(X=true |
+          parent1=v1, parent2=v2, ...) = p. Each key must have as many
+          values as there are parents. You can use this form always;
+          the first two are just conveniences.
+
+        In all cases the probability of X being false is left implicit,
+        since it follows from P(X=true).
+
+        >>> X = BayesNode('X', '', 0.2)
+        >>> Y = BayesNode('Y', 'P', {T: 0.2, F: 0.7})
+        >>> Z = BayesNode('Z', 'P Q',
+        ...    {(T, T): 0.2, (T, F): 0.3, (F, T): 0.5, (F, F): 0.7})
+        """
+        if isinstance(parents, str):
+            parents = parents.split()
+
+        # We store the table always in the third form above.
+        if isinstance(cpt, (float, int)):  # no parents, 0-tuple
+            cpt = {(): cpt}
+        elif isinstance(cpt, dict):
+            # one parent, 1-tuple
+            if cpt and isinstance(list(cpt.keys())[0], bool):
+                cpt = {(v,): p for v, p in cpt.items()}
+
+        assert isinstance(cpt, dict)
+        for vs, p in cpt.items():
+            assert isinstance(vs, tuple) and len(vs) == len(parents)
+            assert all(isinstance(v, bool) for v in vs)
+            assert 0 <= p <= 1
+
+        self.variable = X
+        self.parents = parents
+        self.cpt = cpt
+        self.children = []
+
+    def p(self, value, event):
+        """
+        Return the conditional probability
+        P(X=value | parents=parent_values), where parent_values
+        are the values of parents in event. (event must assign each
+        parent a value.)
+        >>> bn = BayesNode('X', 'Burglary', {T: 0.2, F: 0.625})
+        >>> bn.p(False, {'Burglary': False, 'Earthquake': True})
+        0.375
+        """
+        assert isinstance(value, bool)
+        ptrue = self.cpt[event_values(event, self.parents)]
+        return ptrue if value else 1 - ptrue
+
+    def sample(self, event):
+        """
+        Sample from the distribution for this variable conditioned
+        on event's values for parent_variables. That is, return True/False
+        at random according with the conditional probability given the
+        parents.
+        """
+        return probability(self.p(True, event))
+
+    def __repr__(self):
+        return repr((self.variable, ' '.join(self.parents)))
+
+# Burglary example [Figure 13 .2]
+
+
+T, F = True, False
+
+burglary = BayesNet([
+    ('Burglary', '', 0.001),
+    ('Earthquake', '', 0.002),
+    ('Alarm', 'Burglary Earthquake',
+     {(T, T): 0.95, (T, F): 0.94, (F, T): 0.29, (F, F): 0.001}),
+    ('JohnCalls', 'Alarm', {T: 0.90, F: 0.05}),
+    ('MaryCalls', 'Alarm', {T: 0.70, F: 0.01})
+])
+
+# ______________________________________________________________________________
+# Section 13.2. The Semantics of Bayesian Networks
+# Bayesian nets with continuous variables
+
+
+def gaussian_probability(param, event, value):
+    """
+    Gaussian probability of a continuous Bayesian network node on condition of
+    certain event and the parameters determined by the event
+    :param param: parameters determined by discrete parent events of current node
+    :param event: a dict, continuous event of current node, the values are used
+                  as parameters in calculating distribution
+    :param value: float, the value of current continuous node
+    :return: float, the calculated probability
+    >>> param = {'sigma':0.5, 'b':1, 'a':{'h1':0.5, 'h2': 1.5}}
+    >>> event = {'h1':0.6, 'h2': 0.3}
+    >>> gaussian_probability(param, event, 1)
+    0.2590351913317835
+    """
+
+    assert isinstance(event, dict)
+    assert isinstance(param, dict)
+    buff = 0
+    for k, v in event.items():
+        # buffer varianle to calculate h1*a_h1 + h2*a_h2
+        buff += param['a'][k] * v
+    res = 1/(param['sigma']*sqrt(2*pi)) * exp(-0.5*((value-buff-param['b'])/param['sigma'])**2)
+    return res
+
+
+def logistic_probability(param, event, value):
+    """
+    Logistic probability of a discrete node in Bayesian network with continuous parents,
+    :param param: a dict, parameters determined by discrete parents of current node
+    :param event: a dict, names and values of continuous parent variables of current node
+    :param value: boolean, True or False
+    :return: int, probability
+    """
+
+    buff = 1
+    for _,v in event.items():
+        # buffer variable to calculate (value-mu)/sigma
+
+        buff *= (v-param['mu'])/param['sigma']
+    p = 1 - 1/(1+exp(-4/sqrt(2*pi)*buff))
+    return p if value else 1-p
+
+
+class ContinuousBayesNode:
+    """ A Bayesian network node with continuous distribution or with continuous distributed parents """
+
+    def __init__(self, name, d_parents, c_parents, parameters, type):
+        """
+        A continuous Bayesian node has two types of parents: discrete and continuous.
+        :param d_parents: str, name of discrete parents, value of which determines distribution parameters
+        :param c_parents: str, name of continuous parents, value of which is used to calculate distribution
+        :param parameters: a dict, parameters for distribution of current node, keys corresponds to discrete parents
+        :param type: str, type of current node's value, either 'd' (discrete) or 'c'(continuous)
+        """
+
+        self.parameters = parameters
+        self.type = type
+        self.d_parents = d_parents.split()
+        self.c_parents = c_parents.split()
+        self.parents = self.d_parents + self.c_parents
+        self.variable = name
+        self.children = []
+
+    def continuous_p(self, value, c_event, d_event):
+        """
+        Probability given the value of current node and its parents
+        :param c_event: event of continuous nodes
+        :param d_event: event of discrete nodes
+        """
+        assert isinstance(c_event, dict)
+        assert isinstance(d_event, dict)
+
+        d_event_vals = event_values(d_event, self.d_parents)
+        if len(d_event_vals) == 1:
+            d_event_vals = d_event_vals[0]
+        param = self.parameters[d_event_vals]
+        if self.type == "c":
+            p = gaussian_probability(param, c_event, value)
+        if self.type == "d":
+            p = logistic_probability(param, c_event, value)
+        return p
+
+# harvest-buy example. Figure 13.5
+
+
+harvest_buy = BayesNet([
+    ('Subsidy', '', 0.001),
+    ('Harvest', '', 0.002),
+    ('Cost', 'Subsidy', 'Harvest',
+     {True: {'sigma': 0.5, 'b': 1, 'a': {'Harvest': 0.5}},
+      False: {'sigma': 0.6, 'b': 1, 'a': {'Harvest': 0.5}}}, 'c'),
+    ('Buys', '', 'Cost', {T: {'mu':0.5, 'sigma':0.5}, F: {'mu': 0.6, 'sigma':0.6}}, 'd'),
+])
+
+
+# ______________________________________________________________________________
+# 13.3 Exact Inference in Bayesian Networks
+# 13.3.1 Inference by enumeration
+
+
+def enumeration_ask(X, e, bn):
+    """
+    Return the conditional probability distribution of variable X
+    given evidence e, from BayesNet bn. [Figure 13.10]
+    >>> enumeration_ask('Burglary', dict(JohnCalls=T, MaryCalls=T), burglary
+    ...  ).show_approx()
+    'False: 0.716, True: 0.284'
+    """
+
+    assert X not in e, "Query variable must be distinct from evidence"
+    Q = ProbDist(X)
+    for xi in bn.variable_values(X):
+        Q[xi] = enumerate_all(bn.variables, extend(e, X, xi), bn)
+    return Q.normalize()
+
+
+def enumerate_all(variables, e, bn):
+    """
+    Return the sum of those entries in P(variables | e{others})
+    consistent with e, where P is the joint distribution represented
+    by bn, and e{others} means e restricted to bn's other variables
+    (the ones other than variables). Parents must precede children in variables.
+    """
+
+    if not variables:
+        return 1.0
+    Y, rest = variables[0], variables[1:]
+    Ynode = bn.variable_node(Y)
+    if Y in e:
+        return Ynode.p(e[Y], e) * enumerate_all(rest, e, bn)
+    else:
+        return sum(Ynode.p(y, e) * enumerate_all(rest, extend(e, Y, y), bn)
+                   for y in bn.variable_values(Y))
+
+# ______________________________________________________________________________
+# 13.3.2 The variable elimination algorithm
+
+
+def elimination_ask(X, e, bn):
+    """
+    Compute bn's P(X|e) by variable elimination. [Figure 13.12]
+    >>> elimination_ask('Burglary', dict(JohnCalls=T, MaryCalls=T), burglary
+    ...  ).show_approx()
+    'False: 0.716, True: 0.284'
+    """
+    assert X not in e, "Query variable must be distinct from evidence"
+    factors = []
+    for var in reversed(bn.variables):
+        factors.append(make_factor(var, e, bn))
+        if is_hidden(var, X, e):
+            factors = sum_out(var, factors, bn)
+    return pointwise_product(factors, bn).normalize()
+
+
+def is_hidden(var, X, e):
+    """Is var a hidden variable when querying P(X|e)?"""
+    return var != X and var not in e
+
+
+def make_factor(var, e, bn):
+    """
+    Return the factor for var in bn's joint distribution given e.
+    That is, bn's full joint distribution, projected to accord with e,
+    is the pointwise product of these factors for bn's variables.
+    """
+    node = bn.variable_node(var)
+    variables = [X for X in [var] + node.parents if X not in e]
+    cpt = {event_values(e1, variables): node.p(e1[var], e1)
+           for e1 in all_events(variables, bn, e)}
+    return Factor(variables, cpt)
+
+
+def pointwise_product(factors, bn):
+    return reduce(lambda f, g: f.pointwise_product(g, bn), factors)
+
+
+def sum_out(var, factors, bn):
+    """Eliminate var from all factors by summing over its values."""
+    result, var_factors = [], []
+    for f in factors:
+        (var_factors if var in f.variables else result).append(f)
+    result.append(pointwise_product(var_factors, bn).sum_out(var, bn))
+    return result
+
+
+class Factor:
+    """A factor in a joint distribution."""
+
+    def __init__(self, variables, cpt):
+        self.variables = variables
+        self.cpt = cpt
+
+    def pointwise_product(self, other, bn):
+        """Multiply two factors, combining their variables."""
+        variables = list(set(self.variables) | set(other.variables))
+        cpt = {event_values(e, variables): self.p(e) * other.p(e)
+               for e in all_events(variables, bn, {})}
+        return Factor(variables, cpt)
+
+    def sum_out(self, var, bn):
+        """Make a factor eliminating var by summing over its values."""
+        variables = [X for X in self.variables if X != var]
+        cpt = {event_values(e, variables): sum(self.p(extend(e, var, val))
+                                               for val in bn.variable_values(var))
+               for e in all_events(variables, bn, {})}
+        return Factor(variables, cpt)
+
+    def normalize(self):
+        """Return my probabilities; must be down to one variable."""
+        assert len(self.variables) == 1
+        return ProbDist(self.variables[0],
+                        {k: v for ((k,), v) in self.cpt.items()})
+
+    def p(self, e):
+        """Look up my value tabulated for e."""
+        return self.cpt[event_values(e, self.variables)]
+
+
+def all_events(variables, bn, e):
+    """Yield every way of extending e with values for all variables."""
+    if not variables:
+        yield e
+    else:
+        X, rest = variables[0], variables[1:]
+        for e1 in all_events(rest, bn, e):
+            for x in bn.variable_values(X):
+                yield extend(e1, X, x)
+
+# ______________________________________________________________________________
+# 13.3.4 Clustering algorithms
+# [Figure 13.14a]: sprinkler network
+
+
+sprinkler = BayesNet([
+    ('Cloudy', '', 0.5),
+    ('Sprinkler', 'Cloudy', {T: 0.10, F: 0.50}),
+    ('Rain', 'Cloudy', {T: 0.80, F: 0.20}),
+    ('WetGrass', 'Sprinkler Rain',
+     {(T, T): 0.99, (T, F): 0.90, (F, T): 0.90, (F, F): 0.00})])
+
+# ______________________________________________________________________________
+# 13.4 Approximate Inference for Bayesian Networks
+# 13.4.1 Direct sampling methods
+
+
+def prior_sample(bn):
+    """
+    Randomly sample from bn's full joint distribution. The result
+    is a {variable: value} dict. [Figure 13.15]
+    """
+    event = {}
+    for node in bn.nodes:
+        event[node.variable] = node.sample(event)
+    return event
+
+# _________________________________________________________________________
+
+
+def rejection_sampling(X, e, bn, N=10000):
+    """
+    Estimate the probability distribution of variable X given
+    evidence e in BayesNet bn, using N samples.  [Figure 13.16]
+    Raises a ZeroDivisionError if all the N samples are rejected,
+    i.e., inconsistent with e.
+    >>> random.seed(47)
+    >>> rejection_sampling('Burglary', dict(JohnCalls=T, MaryCalls=T),
+    ...   burglary, 10000).show_approx()
+    'False: 0.7, True: 0.3'
+    """
+    counts = {x: 0 for x in bn.variable_values(X)}  # bold N in [Figure 13.16]
+    for j in range(N):
+        sample = prior_sample(bn)  # boldface x in [Figure 13.16]
+        if consistent_with(sample, e):
+            counts[sample[X]] += 1
+    return ProbDist(X, counts)
+
+
+def consistent_with(event, evidence):
+    """Is event consistent with the given evidence?"""
+    return all(evidence.get(k, v) == v
+               for k, v in event.items())
+
+# _________________________________________________________________________
+
+
+def likelihood_weighting(X, e, bn, N=10000):
+    """
+    Estimate the probability distribution of variable X given
+    evidence e in BayesNet bn.  [Figure 13.17]
+    >>> random.seed(1017)
+    >>> likelihood_weighting('Burglary', dict(JohnCalls=T, MaryCalls=T),
+    ...   burglary, 10000).show_approx()
+    'False: 0.702, True: 0.298'
+    """
+
+    W = {x: 0 for x in bn.variable_values(X)}
+    for j in range(N):
+        sample, weight = weighted_sample(bn, e)  # boldface x, w in [Figure 14.15]
+        W[sample[X]] += weight
+    return ProbDist(X, W)
+
+
+def weighted_sample(bn, e):
+    """
+    Sample an event from bn that's consistent with the evidence e;
+    return the event and its weight, the likelihood that the event
+    accords to the evidence.
+    """
+
+    w = 1
+    event = dict(e)  # boldface x in [Figure 13.17]
+    for node in bn.nodes:
+        Xi = node.variable
+        if Xi in e:
+            w *= node.p(e[Xi], event)
+        else:
+            event[Xi] = node.sample(event)
+    return event, w
+
+# _________________________________________________________________________
+# 13.4.2 Inference by Markov chain simulation
+
+
+def gibbs_ask(X, e, bn, N=1000):
+    """[Figure 13.19]"""
+    assert X not in e, "Query variable must be distinct from evidence"
+    counts = {x: 0 for x in bn.variable_values(X)}  # bold N in [Figure 14.16]
+    Z = [var for var in bn.variables if var not in e]
+    state = dict(e)  # boldface x in [Figure 14.16]
+    for Zi in Z:
+        state[Zi] = random.choice(bn.variable_values(Zi))
+    for j in range(N):
+        for Zi in Z:
+            state[Zi] = markov_blanket_sample(Zi, state, bn)
+            counts[state[X]] += 1
+    return ProbDist(X, counts)
+
+
+def markov_blanket_sample(X, e, bn):
+    """
+    Return a sample from P(X | mb) where mb denotes that the
+    variables in the Markov blanket of X take their values from event
+    e (which must assign a value to each). The Markov blanket of X is
+    X's parents, children, and children's parents.
+    """
+    Xnode = bn.variable_node(X)
+    Q = ProbDist(X)
+    for xi in bn.variable_values(X):
+        ei = extend(e, X, xi)
+        # [Equation 13.12:]
+        Q[xi] = Xnode.p(xi, e) * product(Yj.p(ei[Yj.variable], ei)
+                                         for Yj in Xnode.children)
+    # (assuming a Boolean variable here)
+    return probability(Q.normalize()[True])
+
+# _________________________________________________________________________
+# 13.4.3 Compiling approximate inference
+
+
+class complied_burglary:
+    """compiled version of burglary network"""
+
+    def Burglary(self, sample):
+        if sample['Alarm']:
+            if sample['Earthquake']:
+                return probability(0.00327)
+            else:
+                return probability(0.485)
+        else:
+            if sample['Earthquake']:
+                return probability(7.05e-05)
+            else:
+                return probability(6.01e-05)
+
+    def Earthquake(self, sample):
+        if sample['Alarm']:
+            if sample['Burglary']:
+                return probability(0.0020212)
+            else:
+                return probability(0.36755)
+        else:
+            if sample['Burglary']:
+                return probability(0.0016672)
+            else:
+                return probability(0.0014222)
+
+    def MaryCalls(self, sample):
+        if sample['Alarm']:
+            return probability(0.7)
+        else:
+            return probability(0.01)
+
+    def JongCalls(self, sample):
+        if sample['Alarm']:
+            return probability(0.9)
+        else:
+            return probability(0.05)
+
+    def Alarm(self, sample):
+        raise NotImplementedError
diff --git a/tests/test_probability4e.py b/tests/test_probability4e.py
new file mode 100644
index 000000000..1ce4d7660
--- /dev/null
+++ b/tests/test_probability4e.py
@@ -0,0 +1,342 @@
+from probability4e import *
+
+
+def tests():
+    cpt = burglary.variable_node('Alarm')
+    event = {'Burglary': True, 'Earthquake': True}
+    assert cpt.p(True, event) == 0.95
+    event = {'Burglary': False, 'Earthquake': True}
+    assert cpt.p(False, event) == 0.71
+    # #enumeration_ask('Earthquake', {}, burglary)
+
+    s = {'A': True, 'B': False, 'C': True, 'D': False}
+    assert consistent_with(s, {})
+    assert consistent_with(s, s)
+    assert not consistent_with(s, {'A': False})
+    assert not consistent_with(s, {'D': True})
+
+    random.seed(21)
+    p = rejection_sampling('Earthquake', {}, burglary, 1000)
+    assert p[True], p[False] == (0.001, 0.999)
+
+    random.seed(71)
+    p = likelihood_weighting('Earthquake', {}, burglary, 1000)
+    assert p[True], p[False] == (0.002, 0.998)
+
+# test ProbDist
+
+
+def test_probdist_basic():
+    P = ProbDist('Flip')
+    P['H'], P['T'] = 0.25, 0.75
+    assert P['H'] == 0.25
+    assert P['T'] == 0.75
+    assert P['X'] == 0.00
+
+    P = ProbDist('BiasedDie')
+    P['1'], P['2'], P['3'], P['4'], P['5'], P['6'] = 10, 15, 25, 30, 40, 80
+    P.normalize()
+    assert P['2'] == 0.075
+    assert P['4'] == 0.15
+    assert P['6'] == 0.4
+
+
+def test_probdist_frequency():
+    P = ProbDist('X', {'lo': 125, 'med': 375, 'hi': 500})
+    assert (P['lo'], P['med'], P['hi']) == (0.125, 0.375, 0.5)
+
+    P = ProbDist('Pascal-5', {'x1': 1, 'x2': 5, 'x3': 10, 'x4': 10, 'x5': 5, 'x6': 1})
+    assert (P['x1'], P['x2'], P['x3'], P['x4'], P['x5'], P['x6']) == (
+            0.03125, 0.15625, 0.3125, 0.3125, 0.15625, 0.03125)
+
+
+def test_probdist_normalize():
+    P = ProbDist('Flip')
+    P['H'], P['T'] = 35, 65
+    P = P.normalize()
+    assert (P.prob['H'], P.prob['T']) == (0.350, 0.650)
+
+    P = ProbDist('BiasedDie')
+    P['1'], P['2'], P['3'], P['4'], P['5'], P['6'] = 10, 15, 25, 30, 40, 80
+    P = P.normalize()
+    assert (P.prob['1'], P.prob['2'], P.prob['3'], P.prob['4'], P.prob['5'], P.prob['6']) == (
+                                                    0.05, 0.075, 0.125, 0.15, 0.2, 0.4)
+
+# test JoinProbDist
+
+
+def test_jointprob():
+    P = JointProbDist(['X', 'Y'])
+    P[1, 1] = 0.25
+    assert P[1, 1] == 0.25
+    P[dict(X=0, Y=1)] = 0.5
+    assert P[dict(X=0, Y=1)] == 0.5
+
+
+def test_event_values():
+    assert event_values({'A': 10, 'B': 9, 'C': 8}, ['C', 'A']) == (8, 10)
+    assert event_values((1, 2), ['C', 'A']) == (1, 2)
+
+
+def test_enumerate_joint():
+    P = JointProbDist(['X', 'Y'])
+    P[0, 0] = 0.25
+    P[0, 1] = 0.5
+    P[1, 1] = P[2, 1] = 0.125
+    assert enumerate_joint(['Y'], dict(X=0), P) == 0.75
+    assert enumerate_joint(['X'], dict(Y=2), P) == 0
+    assert enumerate_joint(['X'], dict(Y=1), P) == 0.75
+
+    Q = JointProbDist(['W', 'X', 'Y', 'Z'])
+    Q[0, 1, 1, 0] = 0.12
+    Q[1, 0, 1, 1] = 0.4
+    Q[0, 0, 1, 1] = 0.5
+    Q[0, 0, 1, 0] = 0.05
+    Q[0, 0, 0, 0] = 0.675
+    Q[1, 1, 1, 0] = 0.3
+    assert enumerate_joint(['W'], dict(X=0, Y=0, Z=1), Q) == 0
+    assert enumerate_joint(['W'], dict(X=0, Y=0, Z=0), Q) == 0.675
+    assert enumerate_joint(['W'], dict(X=0, Y=1, Z=1), Q) == 0.9
+    assert enumerate_joint(['Y'], dict(W=1, X=0, Z=1), Q) == 0.4
+    assert enumerate_joint(['Z'], dict(W=0, X=0, Y=0), Q) == 0.675
+    assert enumerate_joint(['Z'], dict(W=1, X=1, Y=1), Q) == 0.3
+
+
+def test_enumerate_joint_ask():
+    P = JointProbDist(['X', 'Y'])
+    P[0, 0] = 0.25
+    P[0, 1] = 0.5
+    P[1, 1] = P[2, 1] = 0.125
+    assert enumerate_joint_ask(
+            'X', dict(Y=1), P).show_approx() == '0: 0.667, 1: 0.167, 2: 0.167'
+
+
+def test_is_independent():
+    P = JointProbDist(['X', 'Y'])
+    P[0, 0] = P[0,1] = P[1, 1] = P[1, 0] = 0.25
+    assert enumerate_joint_ask(
+        'X', dict(Y=1), P).show_approx() == '0: 0.5, 1: 0.5'
+    assert is_independent(['X','Y'], P)
+
+# test BayesNode
+
+
+def test_bayesnode_p():
+    bn = BayesNode('X', 'Burglary', {T: 0.2, F: 0.625})
+    assert bn.p(True, {'Burglary': True, 'Earthquake': False}) == 0.2
+    assert bn.p(False, {'Burglary': False, 'Earthquake': True}) == 0.375
+    assert BayesNode('W', '', 0.75).p(False, {'Random': True}) == 0.25
+
+
+def test_bayesnode_sample():
+    X = BayesNode('X', 'Burglary', {T: 0.2, F: 0.625})
+    assert X.sample({'Burglary': False, 'Earthquake': True}) in [True, False]
+    Z = BayesNode('Z', 'P Q', {(True, True): 0.2, (True, False): 0.3,
+                               (False, True): 0.5, (False, False): 0.7})
+    assert Z.sample({'P': True, 'Q': False}) in [True, False]
+
+# test continuous variable bayesian net
+
+
+def test_gaussian_probability():
+    param = {'sigma': 0.5, 'b': 1, 'a': {'h': 0.5}}
+    event = {'h': 0.6}
+    assert gaussian_probability(param, event, 1) == 0.6664492057835993
+
+
+def test_logistic_probability():
+    param = {'mu': 0.5, 'sigma': 0.1}
+    event = {'h': 0.6}
+    assert logistic_probability(param, event, True) == 0.16857376940725355
+    assert logistic_probability(param, event, False) == 0.8314262305927465
+
+
+def test_enumeration_ask():
+    assert enumeration_ask(
+            'Burglary', dict(JohnCalls=T, MaryCalls=T),
+            burglary).show_approx() == 'False: 0.716, True: 0.284'
+    assert enumeration_ask(
+            'Burglary', dict(JohnCalls=T, MaryCalls=F),
+            burglary).show_approx() == 'False: 0.995, True: 0.00513'
+    assert enumeration_ask(
+            'Burglary', dict(JohnCalls=F, MaryCalls=T),
+            burglary).show_approx() == 'False: 0.993, True: 0.00688'
+    assert enumeration_ask(
+            'Burglary', dict(JohnCalls=T),
+            burglary).show_approx() == 'False: 0.984, True: 0.0163'
+    assert enumeration_ask(
+            'Burglary', dict(MaryCalls=T),
+            burglary).show_approx() == 'False: 0.944, True: 0.0561'
+
+
+def test_elimination_ask():
+    assert elimination_ask(
+            'Burglary', dict(JohnCalls=T, MaryCalls=T),
+            burglary).show_approx() == 'False: 0.716, True: 0.284'
+    assert elimination_ask(
+            'Burglary', dict(JohnCalls=T, MaryCalls=F),
+            burglary).show_approx() == 'False: 0.995, True: 0.00513'
+    assert elimination_ask(
+            'Burglary', dict(JohnCalls=F, MaryCalls=T),
+            burglary).show_approx() == 'False: 0.993, True: 0.00688'
+    assert elimination_ask(
+            'Burglary', dict(JohnCalls=T),
+            burglary).show_approx() == 'False: 0.984, True: 0.0163'
+    assert elimination_ask(
+            'Burglary', dict(MaryCalls=T),
+            burglary).show_approx() == 'False: 0.944, True: 0.0561'
+
+
+# test sampling
+
+
+def test_prior_sample():
+    random.seed(42)
+    all_obs = [prior_sample(burglary) for x in range(1000)]
+    john_calls_true = [observation for observation in all_obs if observation['JohnCalls'] == True]
+    mary_calls_true = [observation for observation in all_obs if observation['MaryCalls'] == True]
+    burglary_and_john = [observation for observation in john_calls_true if observation['Burglary'] == True]
+    burglary_and_mary = [observation for observation in mary_calls_true if observation['Burglary'] == True]
+    assert len(john_calls_true) / 1000 == 46 / 1000
+    assert len(mary_calls_true) / 1000 == 13 / 1000
+    assert len(burglary_and_john) / len(john_calls_true) == 1 / 46
+    assert len(burglary_and_mary) / len(mary_calls_true) == 1 / 13
+
+
+def test_prior_sample2():
+    random.seed(128)
+    all_obs = [prior_sample(sprinkler) for x in range(1000)]
+    rain_true = [observation for observation in all_obs if observation['Rain'] == True]
+    sprinkler_true = [observation for observation in all_obs if observation['Sprinkler'] == True]
+    rain_and_cloudy = [observation for observation in rain_true if observation['Cloudy'] == True]
+    sprinkler_and_cloudy = [observation for observation in sprinkler_true if observation['Cloudy'] == True]
+    assert len(rain_true) / 1000 == 0.476
+    assert len(sprinkler_true) / 1000 == 0.291
+    assert len(rain_and_cloudy) / len(rain_true) == 376 / 476
+    assert len(sprinkler_and_cloudy) / len(sprinkler_true) == 39 / 291
+
+
+def test_rejection_sampling():
+    random.seed(47)
+    assert rejection_sampling(
+            'Burglary', dict(JohnCalls=T, MaryCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.7, True: 0.3'
+    assert rejection_sampling(
+            'Burglary', dict(JohnCalls=T, MaryCalls=F),
+            burglary, 10000).show_approx() == 'False: 1, True: 0'
+    assert rejection_sampling(
+            'Burglary', dict(JohnCalls=F, MaryCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.987, True: 0.0128'
+    assert rejection_sampling(
+            'Burglary', dict(JohnCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.982, True: 0.0183'
+    assert rejection_sampling(
+            'Burglary', dict(MaryCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.965, True: 0.0348'
+
+
+def test_rejection_sampling2():
+    random.seed(42)
+    assert rejection_sampling(
+            'Cloudy', dict(Rain=T, Sprinkler=T),
+            sprinkler, 10000).show_approx() == 'False: 0.56, True: 0.44'
+    assert rejection_sampling(
+            'Cloudy', dict(Rain=T, Sprinkler=F),
+            sprinkler, 10000).show_approx() == 'False: 0.119, True: 0.881'
+    assert rejection_sampling(
+            'Cloudy', dict(Rain=F, Sprinkler=T),
+            sprinkler, 10000).show_approx() == 'False: 0.951, True: 0.049'
+    assert rejection_sampling(
+            'Cloudy', dict(Rain=T),
+            sprinkler, 10000).show_approx() == 'False: 0.205, True: 0.795'
+    assert rejection_sampling(
+            'Cloudy', dict(Sprinkler=T),
+            sprinkler, 10000).show_approx() == 'False: 0.835, True: 0.165'
+
+
+def test_likelihood_weighting():
+    random.seed(1017)
+    assert likelihood_weighting(
+            'Burglary', dict(JohnCalls=T, MaryCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.702, True: 0.298'
+    assert likelihood_weighting(
+            'Burglary', dict(JohnCalls=T, MaryCalls=F),
+            burglary, 10000).show_approx() == 'False: 0.993, True: 0.00656'
+    assert likelihood_weighting(
+            'Burglary', dict(JohnCalls=F, MaryCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.996, True: 0.00363'
+    assert likelihood_weighting(
+            'Burglary', dict(JohnCalls=F, MaryCalls=F),
+            burglary, 10000).show_approx() == 'False: 1, True: 0.000126'
+    assert likelihood_weighting(
+            'Burglary', dict(JohnCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.979, True: 0.0205'
+    assert likelihood_weighting(
+            'Burglary', dict(MaryCalls=T),
+            burglary, 10000).show_approx() == 'False: 0.94, True: 0.0601'
+
+
+def test_likelihood_weighting2():
+    random.seed(42)
+    assert likelihood_weighting(
+            'Cloudy', dict(Rain=T, Sprinkler=T),
+            sprinkler, 10000).show_approx() == 'False: 0.559, True: 0.441'
+    assert likelihood_weighting(
+            'Cloudy', dict(Rain=T, Sprinkler=F),
+            sprinkler, 10000).show_approx() == 'False: 0.12, True: 0.88'
+    assert likelihood_weighting(
+            'Cloudy', dict(Rain=F, Sprinkler=T),
+            sprinkler, 10000).show_approx() == 'False: 0.951, True: 0.0486'
+    assert likelihood_weighting(
+            'Cloudy', dict(Rain=T),
+            sprinkler, 10000).show_approx() == 'False: 0.198, True: 0.802'
+    assert likelihood_weighting(
+            'Cloudy', dict(Sprinkler=T),
+            sprinkler, 10000).show_approx() == 'False: 0.833, True: 0.167'
+
+
+def test_gibbs_ask():
+
+    g_solution = gibbs_ask('Cloudy', dict(Rain=True), sprinkler, 1000)
+    assert abs(g_solution.prob[False]-0.2) < 0.05
+    assert abs(g_solution.prob[True]-0.8) < 0.05
+
+
+# The following should probably go in .ipynb:
+
+"""
+# We can build up a probability distribution like this (p. 469):
+>>> P = ProbDist()
+>>> P['sunny'] = 0.7
+>>> P['rain'] = 0.2
+>>> P['cloudy'] = 0.08
+>>> P['snow'] = 0.02
+
+# and query it like this:  (Never mind this ELLIPSIS option
+#                           added to make the doctest portable.)
+>>> P['rain']               #doctest:+ELLIPSIS
+0.2...
+
+# A Joint Probability Distribution is dealt with like this [Figure 13.3]:
+>>> P = JointProbDist(['Toothache', 'Cavity', 'Catch'])
+>>> T, F = True, False
+>>> P[T, T, T] = 0.108; P[T, T, F] = 0.012; P[F, T, T] = 0.072; P[F, T, F] = 0.008
+>>> P[T, F, T] = 0.016; P[T, F, F] = 0.064; P[F, F, T] = 0.144; P[F, F, F] = 0.576
+
+>>> P[T, T, T]
+0.108
+
+# Ask for P(Cavity|Toothache=T)
+>>> PC = enumerate_joint_ask('Cavity', {'Toothache': T}, P)
+>>> PC.show_approx()
+'False: 0.4, True: 0.6'
+
+>>> 0.6-epsilon < PC[T] < 0.6+epsilon
+True
+
+>>> 0.4-epsilon < PC[F] < 0.4+epsilon
+True
+"""
+
+if __name__ == '__main__':
+    pytest.main()