diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index 5783609765..08bfcf8a0f 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -1,16 +1,16 @@
-import collections
 import numbers
-
 import numpy as np
 import theano.tensor as tt
 from theano import function
 import theano
+from networkx import DiGraph, topological_sort
 from ..memoize import memoize
 from ..model import (
-    Model, get_named_nodes_and_relations, FreeRV,
-    ObservedRV, MultiObservedRV
+    Model, modelcontext, FreeRV, ObservedRV,
+    not_shared_or_constant_variable, get_sub_dag
 )
 from ..vartypes import string_types
+from ..util import WrapAsHashable
 
 __all__ = ['DensityDist', 'Distribution', 'Continuous', 'Discrete',
            'NoDistribution', 'TensorType', 'draw_values', 'generate_samples']
@@ -211,10 +211,10 @@ def random(self, *args, **kwargs):
             return self.rand(*args, **kwargs)
         else:
             raise ValueError("Distribution was not passed any random method "
-                            "Define a custom random method and pass it as kwarg random")
+                             "Define a custom random method and pass it as kwarg random")
 
 
-def draw_values(params, point=None, size=None):
+def draw_values(params, point=None, size=None, model=None):
     """
     Draw (fix) parameter values. Handles a number of cases:
 
@@ -232,99 +232,110 @@ def draw_values(params, point=None, size=None):
             b) are *RVs with a random method
 
     """
-    # Distribution parameters may be nodes which have named node-inputs
-    # specified in the point. Need to find the node-inputs, their
-    # parents and children to replace them.
-    leaf_nodes = {}
-    named_nodes_parents = {}
-    named_nodes_children = {}
-    for param in params:
-        if hasattr(param, 'name'):
-            # Get the named nodes under the `param` node
-            nn, nnp, nnc = get_named_nodes_and_relations(param)
-            leaf_nodes.update(nn)
-            # Update the discovered parental relationships
-            for k in nnp.keys():
-                if k not in named_nodes_parents.keys():
-                    named_nodes_parents[k] = nnp[k]
-                else:
-                    named_nodes_parents[k].update(nnp[k])
-            # Update the discovered child relationships
-            for k in nnc.keys():
-                if k not in named_nodes_children.keys():
-                    named_nodes_children[k] = nnc[k]
-                else:
-                    named_nodes_children[k].update(nnc[k])
-
-    # Init givens and the stack of nodes to try to `_draw_value` from
-    givens = {}
-    stored = set()  # Some nodes
-    stack = list(leaf_nodes.values())  # A queue would be more appropriate
-    while stack:
-        next_ = stack.pop(0)
-        if next_ in stored:
-            # If the node already has a givens value, skip it
-            continue
-        elif isinstance(next_, (tt.TensorConstant,
-                                tt.sharedvar.SharedVariable)):
-            # If the node is a theano.tensor.TensorConstant or a
-            # theano.tensor.sharedvar.SharedVariable, its value will be
-            # available automatically in _compile_theano_function so
-            # we can skip it. Furthermore, if this node was treated as a
-            # TensorVariable that should be compiled by theano in
-            # _compile_theano_function, it would raise a `TypeError:
-            # ('Constants not allowed in param list', ...)` for
-            # TensorConstant, and a `TypeError: Cannot use a shared
-            # variable (...) as explicit input` for SharedVariable.
-            stored.add(next_.name)
-            continue
+    # Get the dependence graph, either by copying a subgraph from the model's
+    # dependence_dag or by creating it from scratch
+    try:
+        model = modelcontext(model)
+        graph = model.dependence_dag
+    except Exception:
+        graph = DiGraph()
+    dependence_dag, index = get_sub_dag(graph,
+                                        params)
+    # Store list of nodes in reversed topological sort order, i.e. the nodes
+    # without ancestors at the end. This way, default pop() will get the
+    # correct node, and append() will add a node for the next iteration
+    queue = list(reversed(list(topological_sort(dependence_dag))))
+
+    # Init drawn values and updatable point and givens
+    drawn = {}
+    givens = []
+    if point is None:
+        point = {}
+    else:
+        point = point.copy()
+    nodes_missing_inputs = {}
+
+    while queue:
+        try:
+            # Pop nodes from queue because we may use stack new nodes
+            # onto it in the middle of the loop
+            node = queue.pop()
+        except Exception:
+            break
+
+        # We may have flaged this node to be computed by compiling the
+        # theano function, because deterministic relations take precedence
+        # over conditional relations
+        if isinstance(node, tuple):
+            node, is_determined = node
         else:
-            # If the node does not have a givens value, try to draw it.
-            # The named node's children givens values must also be taken
-            # into account.
-            children = named_nodes_children[next_]
-            temp_givens = [givens[k] for k in givens if k in children]
-            try:
-                # This may fail for autotransformed RVs, which don't
-                # have the random method
-                givens[next_.name] = (next_, _draw_value(next_,
-                                                         point=point,
-                                                         givens=temp_givens,
-                                                         size=size))
-                stored.add(next_.name)
-            except theano.gof.fg.MissingInputError:
-                # The node failed, so we must add the node's parents to
-                # the stack of nodes to try to draw from. We exclude the
-                # nodes in the `params` list.
-                stack.extend([node for node in named_nodes_parents[next_]
-                              if node is not None and
-                              node.name not in stored and
-                              node not in params])
-
-    # the below makes sure the graph is evaluated in order
-    # test_distributions_random::TestDrawValues::test_draw_order fails without it
-    params = dict(enumerate(params))  # some nodes are not hashable
-    evaluated = {}
-    to_eval = set()
-    missing_inputs = set(params)
-    while to_eval or missing_inputs:
-        if to_eval == missing_inputs:
-            raise ValueError('Cannot resolve inputs for {}'.format([str(params[j]) for j in to_eval]))
-        to_eval = set(missing_inputs)
-        missing_inputs = set()
-        for param_idx in to_eval:
-            param = params[param_idx]
-            if hasattr(param, 'name') and param.name in givens:
-                evaluated[param_idx] = givens[param.name][1]
+            is_determined = False
+
+        # Node's value had already been determined so we jump onto the next
+        if node in drawn:
+            continue
+        try:
+            if is_determined:
+                node_value = _compute_value(node,
+                                            givens=givens,
+                                            size=size)
             else:
-                try:  # might evaluate in a bad order,
-                    evaluated[param_idx] = _draw_value(param, point=point, givens=givens.values(), size=size)
-                    if isinstance(param, collections.Hashable) and named_nodes_parents.get(param):
-                        givens[param.name] = (param, evaluated[param_idx])
-                except theano.gof.fg.MissingInputError:
-                    missing_inputs.add(param_idx)
+                node_value = _draw_value(node,
+                                         point=point,
+                                         givens=givens,
+                                         size=size)
+            drawn[node] = node_value
+        except theano.gof.fg.MissingInputError as e:
+            # Expected to fail for auto-transformed RVs whos values were
+            # not provided in point
+            nodes_missing_inputs[node] = e
+            continue
 
-    return [evaluated[j] for j in params] # set the order back
+        # If the node is a theano Variable, which is not TensorConstant,
+        # nor SharedVariable, we must add its value to point (only if it
+        # has conditional children) and givens (only if it has
+        # deterministic children)
+        if not_shared_or_constant_variable(node):
+            edges = dependence_dag[node]
+            if any(not d['deterministic'] for v, d in edges.items()):
+                point[node.name] = node_value
+            if any(d['deterministic'] for v, d in edges.items()):
+                givens.append((node, node_value))
+
+        # If the node has deterministic children, check if the children's
+        # values can be computed. This must be done because deterministic
+        # relations must take precedence over conditional relations
+        # amongst variables.
+        deterministic_children = (child for child, d in
+                                  dependence_dag[node].items()
+                                  if d.get('deterministic', 0))
+        for child in deterministic_children:
+            # Check if all of the child's deterministic parents have their
+            # values set, allowing us to compute the child's value.
+            if not any((p not in drawn for p, c in
+                        dependence_dag.pred[child].items()
+                        if c.get('deterministic', 0))):
+                # Append child to the queue, to
+                # compute its value ahead of its schedule
+                queue.append((child, True))
+
+    # Now that the DAG has been transversed, drawing values, we can place them
+    # in a list following the indexing given by the input list `params`
+    output = []
+    for ind in range(len(params)):
+        node = index[ind]
+        value = drawn.get(node, None)
+        if value is None:
+            # We failed to draw the params[i] value. This could be due to an
+            # ignored MissingInputError, in which case we reraise it, or it
+            # could be some other form of unexpected RuntimeError.
+            if node in nodes_missing_inputs:
+                raise nodes_missing_inputs[node]
+            else:
+                raise RuntimeError('Failed to draw value for parameter {}'.
+                                   format(params[ind]))
+        output.append(value)
+    return output
 
 
 @memoize
@@ -357,7 +368,7 @@ def _draw_value(param, point=None, givens=None, size=None):
 
     Parameters
     ----------
-    param : number, array like, theano variable or pymc3 random variable
+    param : WrapAsHashable, theano variable or pymc3 random variable
         The value or distribution. Constants or shared variables
         will be converted to an array and returned. Theano variables
         are evaluated. If `param` is a pymc3 random variables, draw
@@ -365,23 +376,27 @@ def _draw_value(param, point=None, givens=None, size=None):
         in `point`.
     point : dict, optional
         A dictionary from pymc3 variable names to their values.
-    givens : dict, optional
-        A dictionary from theano variables to their values. These values
+    givens : list, optional
+        A list of tuples of theano variables and their values. These values
         are used to evaluate `param` if it is a theano variable.
     size : int, optional
         Number of samples
     """
-    if isinstance(param, (numbers.Number, np.ndarray)):
-        return param
-    elif isinstance(param, tt.TensorConstant):
-        return param.value
-    elif isinstance(param, tt.sharedvar.SharedVariable):
-        return param.get_value()
-    elif isinstance(param, (tt.TensorVariable, MultiObservedRV)):
+    if isinstance(param, numbers.Number):
+        output = param
+    elif isinstance(param, np.ndarray):
+        output = param
+    elif isinstance(param, theano.tensor.TensorConstant):
+        output = param.value
+    elif isinstance(param, theano.tensor.sharedvar.SharedVariable):
+        output = param.get_value()
+    elif isinstance(param, WrapAsHashable):
+        output = param.get_value()
+    elif isinstance(param, tt.TensorVariable):
         if point and hasattr(param, 'model') and param.name in point:
-            return point[param.name]
+            output = point[param.name]
         elif hasattr(param, 'random') and param.random is not None:
-            return param.random(point=point, size=size)
+            output = param.random(point=point, size=size)
         elif (hasattr(param, 'distribution') and
                 hasattr(param.distribution, 'random') and
                 param.distribution.random is not None):
@@ -400,22 +415,59 @@ def _draw_value(param, point=None, givens=None, size=None):
                     # reset shape to account for shape changes
                     # with theano.shared inputs
                     dist_tmp.shape = np.array([])
-                    val = dist_tmp.random(point=point, size=None)
-                    dist_tmp.shape = val.shape
-                return dist_tmp.random(point=point, size=size)
+                    val = np.atleast_1d(dist_tmp.random(point=point,
+                                                        size=None))
+                    # Sometimes point may change the size of val but not the
+                    # distribution's shape
+                    if point and size is not None:
+                        temp_size = np.atleast_1d(size)
+                        if all(val.shape[:len(temp_size)] == temp_size):
+                            dist_tmp.shape = val.shape[len(temp_size):]
+                        else:
+                            dist_tmp.shape = val.shape
+                output = dist_tmp.random(point=point, size=size)
             else:
-                return param.distribution.random(point=point, size=size)
+                output = param.distribution.random(point=point, size=size)
         else:
-            if givens:
-                variables, values = list(zip(*givens))
-            else:
-                variables = values = []
-            func = _compile_theano_function(param, variables)
-            if size and values and not all(var.dshape == val.shape for var, val in zip(variables, values)):
-                return np.array([func(*v) for v in zip(*values)])
-            else:
-                return func(*values)
-    raise ValueError('Unexpected type in draw_value: %s' % type(param))
+            output = _compute_value(param, givens=givens, size=size)
+    else:
+        raise ValueError('Unexpected type in draw_value: %s' % type(param))
+    # Maybe at this point we should control the output shape depending on size?
+    return output
+
+
+def _compute_value(param, givens=None, size=None):
+    """Compute the value of a theano variable from a givens input list
+
+    Parameters
+    ----------
+    param : Theano variable
+        The variable will be compiled into a function and then evaluated
+        using the givens input list
+    givens : list, optional
+        A list of tuples of theano variables and their values. These values
+        are used to evaluate `param` if it is a theano variable.
+    size : int, optional
+        Number of samples
+    """
+    if givens:
+        variables, values = list(zip(*givens))
+    else:
+        variables = values = []
+    func = _compile_theano_function(param, variables)
+    if size is not None:
+        size = np.atleast_1d(size)
+    dshaped_variables = all((hasattr(var, 'dshape') for var in variables))
+    if (values and dshaped_variables and
+        not all(var.dshape == getattr(val, 'shape', tuple())
+                for var, val in zip(variables, values))):
+        output = np.array([func(*v) for v in zip(*values)])
+    elif (size is not None and any((val.ndim > var.ndim)
+          for var, val in zip(variables, values))):
+        output = np.array([func(*v) for v in zip(*values)])
+    else:
+        output = func(*values)
+    return output
 
 
 def to_tuple(shape):
@@ -424,6 +476,7 @@ def to_tuple(shape):
         return tuple()
     return tuple(np.atleast_1d(shape))
 
+
 def _is_one_d(dist_shape):
     if hasattr(dist_shape, 'dshape') and dist_shape.dshape in ((), (0,), (1,)):
         return True
@@ -433,6 +486,7 @@ def _is_one_d(dist_shape):
         return True
     return False
 
+
 def generate_samples(generator, *args, **kwargs):
     """Generate samples from the distribution of a random variable.
 
diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py
index 5c45d34f8b..37ffeb770a 100755
--- a/pymc3/distributions/multivariate.py
+++ b/pymc3/distributions/multivariate.py
@@ -338,15 +338,17 @@ def __init__(self, nu, Sigma=None, mu=None, cov=None, tau=None, chol=None,
         self.mean = self.median = self.mode = self.mu = self.mu
 
     def random(self, point=None, size=None):
-        nu, mu = draw_values([self.nu, self.mu], point=point, size=size)
         if self._cov_type == 'cov':
-            cov, = draw_values([self.cov], point=point, size=size)
+            nu, mu, cov = draw_values([self.nu, self.mu, self.cov],
+                                      point=point, size=size)
             dist = MvNormal.dist(mu=np.zeros_like(mu), cov=cov)
         elif self._cov_type == 'tau':
-            tau, = draw_values([self.tau], point=point, size=size)
+            nu, mu, tau = draw_values([self.nu, self.mu, self.tau],
+                                      point=point, size=size)
             dist = MvNormal.dist(mu=np.zeros_like(mu), tau=tau)
         else:
-            chol, = draw_values([self.chol_cov], point=point, size=size)
+            nu, mu, chol = draw_values([self.nu, self.mu, self.chol],
+                                       point=point, size=size)
             dist = MvNormal.dist(mu=np.zeros_like(mu), chol=chol)
 
         samples = dist.random(point, size)
diff --git a/pymc3/model.py b/pymc3/model.py
index 2449d70321..2c0a16b1e1 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -2,9 +2,16 @@
 import functools
 import itertools
 import threading
+import re
 import six
+from copy import copy
+try:
+    from collections.abc import Hashable
+except ImportError:
+    from collections import Hashable
 
 import numpy as np
+import networkx
 from pandas import Series
 import scipy.sparse as sps
 import theano.sparse as sparse
@@ -18,7 +25,7 @@
 from .theanof import gradient, hessian, inputvars, generator
 from .vartypes import typefilter, discrete_types, continuous_types, isgenerator
 from .blocking import DictToArrayBijection, ArrayOrdering
-from .util import get_transformed_name
+from .util import get_transformed_name, WrapAsHashable
 
 __all__ = [
     'Model', 'Factor', 'compilef', 'fn', 'fastfn', 'modelcontext',
@@ -79,72 +86,6 @@ def incorporate_methods(source, destination, methods, default=None,
         else:
             setattr(destination, method, None)
 
-def get_named_nodes_and_relations(graph):
-    """Get the named nodes in a theano graph (i.e., nodes whose name
-    attribute is not None) along with their relationships (i.e., the
-    node's named parents, and named children, while skipping unnamed
-    intermediate nodes)
-
-    Parameters
-    ----------
-    graph - a theano node
-
-    Returns:
-    leaf_nodes: A dictionary of name:node pairs, of the named nodes that
-        are also leafs of the graph
-    node_parents: A dictionary of node:set([parents]) pairs. Each key is
-        a theano named node, and the corresponding value is the set of
-        theano named nodes that are parents of the node. These parental
-        relations skip unnamed intermediate nodes.
-    node_children: A dictionary of node:set([children]) pairs. Each key
-        is a theano named node, and the corresponding value is the set
-        of theano named nodes that are children of the node. These child
-        relations skip unnamed intermediate nodes.
-
-    """
-    if graph.name is not None:
-        node_parents = {graph: set()}
-        node_children = {graph: set()}
-    else:
-        node_parents = {}
-        node_children = {}
-    return _get_named_nodes_and_relations(graph, None, {}, node_parents, node_children)
-
-def _get_named_nodes_and_relations(graph, parent, leaf_nodes,
-                                        node_parents, node_children):
-    if getattr(graph, 'owner', None) is None:  # Leaf node
-        if graph.name is not None:  # Named leaf node
-            leaf_nodes.update({graph.name: graph})
-            if parent is not None:  # Is None for the root node
-                try:
-                    node_parents[graph].add(parent)
-                except KeyError:
-                    node_parents[graph] = set([parent])
-                node_children[parent].add(graph)
-            # Flag that the leaf node has no children
-            node_children[graph] = set()
-    else:  # Intermediate node
-        if graph.name is not None:  # Intermediate named node
-            if parent is not None:  # Is only None for the root node
-                try:
-                    node_parents[graph].add(parent)
-                except KeyError:
-                    node_parents[graph] = set([parent])
-                node_children[parent].add(graph)
-            # The current node will be set as the parent of the next
-            # nodes only if it is a named node
-            parent = graph
-            # Init the nodes children to an empty set
-            node_children[graph] = set()
-        for i in graph.owner.inputs:
-            temp_nodes, temp_inter, temp_tree = \
-                _get_named_nodes_and_relations(i, parent, leaf_nodes,
-                                               node_parents, node_children)
-            leaf_nodes.update(temp_nodes)
-            node_parents.update(temp_inter)
-            node_children.update(temp_tree)
-    return leaf_nodes, node_parents, node_children
-
 
 class Context(object):
     """Functionality for objects that put themselves in a context using
@@ -635,6 +576,7 @@ def __new__(cls, *args, **kwargs):
 
     def __init__(self, name='', model=None, theano_config=None):
         self.name = name
+        self.dependence_dag = networkx.DiGraph()
         if self.parent is not None:
             self.named_vars = treedict(parent=self.parent.named_vars)
             self.free_RVs = treelist(parent=self.parent.free_RVs)
@@ -650,6 +592,11 @@ def __init__(self, name='', model=None, theano_config=None):
             self.potentials = treelist()
             self.missing_values = treelist()
 
+    def __setstate__(self, state):
+        self.__dict__.update(state)
+        if 'dependence_dag' not in state:
+            self.dependence_dag = build_dependence_dag_from_model(self)
+
     @property
     def model(self):
         return self
@@ -817,7 +764,8 @@ def Var(self, name, dist, data=None, total_size=None):
                               ' and added transformed {orig_name} to model.'.format(
                                 transform=dist.transform.name,
                                 name=name,
-                                orig_name=get_transformed_name(name, dist.transform)))
+                                orig_name=get_transformed_name(name, dist.transform))
+                              )
                 self.deterministics.append(var)
                 self.add_random_variable(var)
                 return var
@@ -845,7 +793,7 @@ def Var(self, name, dist, data=None, total_size=None):
         self.add_random_variable(var)
         return var
 
-    def add_random_variable(self, var):
+    def add_random_variable(self, var, accept_cons_shared=False):
         """Add a random variable to the named variables of the model."""
         if self.named_vars.tree_contains(var.name):
             raise ValueError(
@@ -853,6 +801,11 @@ def add_random_variable(self, var):
         self.named_vars[var.name] = var
         if not hasattr(self, self.name_of(var.name)):
             setattr(self, self.name_of(var.name), var)
+        # The model should automatically construct a DependenceDAG instance
+        # that encodes the relations between its variables
+        add_to_dependence_dag(self.dependence_dag,
+                              var,
+                              accept_cons_shared=accept_cons_shared)
 
     @property
     def prefix(self):
@@ -1000,7 +953,7 @@ def flatten(self, vars=None, order=None, inputvar=None):
                     inputvar.tag.test_value = flatten_list(vars).tag.test_value
                 else:
                     inputvar.tag.test_value = np.asarray([], inputvar.dtype)
-        replacements = {self.named_vars[name]: inputvar[slc].reshape(shape).astype(dtype)
+        replacements = {self.named_vars[name]: (inputvar[slc].reshape(shape).astype(dtype))
                         for name, slc, shape, dtype in order.vmap}
         view = {vm.var: vm for vm in order.vmap}
         flat_view = FlatView(inputvar, replacements, view)
@@ -1024,8 +977,8 @@ def check_test_point(self, test_point=None, round_vals=2):
         if test_point is None:
             test_point = self.test_point
 
-        return Series({RV.name:np.round(RV.logp(self.test_point), round_vals) for RV in self.basic_RVs},
-            name='Log-probability of test_point')
+        return Series({RV.name: np.round(RV.logp(self.test_point), round_vals) for RV in self.basic_RVs},
+                      name='Log-probability of test_point')
 
     def _repr_latex_(self, name=None, dist=None):
         tex_vars = []
@@ -1098,7 +1051,8 @@ def Point(*args, **kwargs):
 
 
 class FastPointFunc(object):
-    """Wraps so a function so it takes a dict of arguments instead of arguments."""
+    """Wraps so a function so it takes a dict of arguments instead of arguments.
+    """
 
     def __init__(self, f):
         self.f = f
@@ -1399,7 +1353,7 @@ def _walk_up_rv(rv):
 
 def _latex_repr_rv(rv):
     """Make latex string for a Deterministic variable"""
-    return r'$\text{%s} \sim \text{Deterministic}(%s)$' % (rv.name, r',~'.join(_walk_up_rv(rv)))
+    return (r'$\text{%s} \sim \text{Deterministic}(%s)$' % (rv.name, r',~'.join(_walk_up_rv(rv))))
 
 
 def Deterministic(name, var, model=None):
@@ -1417,7 +1371,7 @@ def Deterministic(name, var, model=None):
     model = modelcontext(model)
     var = var.copy(model.name_for(name))
     model.deterministics.append(var)
-    model.add_random_variable(var)
+    model.add_random_variable(var, accept_cons_shared=True)
     var._repr_latex_ = functools.partial(_latex_repr_rv, var)
     var.__latex__ = var._repr_latex_
     return var
@@ -1438,7 +1392,7 @@ def Potential(name, var, model=None):
     model = modelcontext(model)
     var.name = model.name_for(name)
     model.potentials.append(var)
-    model.add_random_variable(var)
+    model.add_random_variable(var, accept_cons_shared=True)
     return var
 
 
@@ -1518,3 +1472,314 @@ def all_continuous(vars):
         return False
     else:
         return True
+
+
+class ConstantNodeException(Exception):
+    pass
+
+
+def not_shared_or_constant_variable(x):
+    return (isinstance(x, theano.Variable) and
+            not(isinstance(x, theano.Constant) or
+                isinstance(x, theano.tensor.sharedvar.SharedVariable))
+            ) or (isinstance(x, (FreeRV, MultiObservedRV, TransformedRV)))
+
+
+def is_autonamed_node(node):
+    _theano_autonamed_ops = [(r'.+\.T$',
+                              theano.tensor.elemwise.DimShuffle,
+                              None),
+                             (r'max$',
+                              theano.tensor.basic.MaxAndArgmax,
+                              None),
+                             (r'argmax$',
+                              theano.tensor.basic.MaxAndArgmax,
+                              None),
+                             (r'mean$',
+                              theano.tensor.elemwise.Elemwise,
+                              'Elemwise{true_div,no_inplace}'),
+                             (r'var$',
+                              theano.tensor.elemwise.Elemwise,
+                              'Elemwise{true_div,no_inplace}'),
+                             (r'std$',
+                              theano.tensor.elemwise.Elemwise,
+                              'Elemwise{sqrt,no_inplace}'),
+                             ]
+    for name_pattern, owner_type, owner_name in _theano_autonamed_ops:
+        if node.name and re.match(name_pattern, node.name):
+            owner = getattr(node, 'owner', None)
+            if owner is not None and isinstance(owner.op, owner_type):
+                if (owner_name is None or
+                        getattr(owner.op, 'name', None) == owner_name):
+                    return True
+    return False
+
+
+def get_first_level_conditionals(root):
+    """Performs a breadth first search on the supplied root node's logpt or
+    transformed logpt graph searching for named input nodes, which are
+    different from the supplied root. Each explored branch will stop
+    either when it reaches the leaf node or when it finds its first named node.
+
+    Parameters
+    ----------
+    root : theano.Variable
+        The node from which to get the transformed.logpt or logpt and perform
+        the search. If root does not have either of these attributes, the
+        function returns None.
+
+    Returns
+    -------
+    conditional_on : set, with named nodes that are not theano.Constant nor
+    SharedVariable. The input `root` is conditionally dependent on these nodes
+    and is one step away from them in the Bayesian network that specifies the
+    relationships, hence the name `get_first_level_conditionals`.
+    """
+    transformed = getattr(root, 'transformed', None)
+    try:
+        cond = transformed.logpt
+    except AttributeError:
+        cond = getattr(root, 'logpt', None)
+    if cond is None:
+        return None
+    conditional_on = set()
+    queue = copy(getattr(cond.owner, 'inputs', []))
+    while queue:
+        parent = queue.pop(0)
+        if (parent is not None and getattr(parent, 'name', None) is not None
+                and not_shared_or_constant_variable(parent)):
+            # We don't include as a conditional relation either logpt depending
+            # on root or on transformed because they are both deterministic
+            # relations.
+            # Also, some theano ops place default names to their associated
+            # output nodes. We want to ignore them, so we test that the name
+            # pattern matches with the ones we don't want, that the op is an
+            # instance of the conflicting op, and, in the case of mean, var,
+            # and std, where the op is built from other ops, that the op.name
+            # matches the default one.
+            if (parent == root or parent == transformed or
+                    is_autonamed_node(parent)):
+                continue
+            conditional_on.add(parent)
+        else:
+            parent_owner = getattr(parent, 'owner', None)
+            queue.extend(getattr(parent_owner, 'inputs', []))
+    if len(conditional_on) == 0:
+        return None
+    return conditional_on
+
+
+def add_to_dependence_dag(dag, node, force=False,
+                          accept_cons_shared=False):
+    """Add a node and its conditional and deterministic parents along with
+    their relations, recursively into the `networkx.DiGraph` instance.
+
+    Parameters
+    ----------
+    dag : networkx.DiGraph instance
+        The Digraph instance on which to add the node. The addition
+        is done inplace.
+    node : The variable to add to the DAG
+        By default `theano.Variable`'s, which could be a pymc random
+        variable, Deterministic or Potential, are allowed. TensorConstants
+        and SharedVariables are not allowed by default, but this behavior
+        can be changed with either the `force` or `accept_cons_shared`
+        inputs.
+        Other unhashable types are only accepted if `force=True`, and they
+        are wrapped by `WrapAsHashable` instances.
+    force : bool (optional)
+        If True, any type of node, except None, is allowed to be added.
+    accept_cons_shared : bool (optional)
+        If True, `theano` `TensorConstant`s and `theano` `SharedVariable`s
+        are allowed to be added to the DAG. These are treated separately
+        a priori because `_draw_value` handles these cases differently.
+
+    Returns
+    -------
+    The node that was actually added into the `dag`. Usually this node is the
+    input node, but depending on the node's type, it could be a
+    `WrapAsHashable` instance.
+    """
+    if node is None:
+        raise TypeError('None is not allowed to be added as a node in '
+                        'variable dependency graph.')
+    if not isinstance(dag, networkx.DiGraph):
+        raise TypeError('Input `dag` must be an instance of networkx.DiGraph. '
+                        'Got {}, which has type {}, instead.'.
+                        format(dag, type(dag)))
+    if not isinstance(node, (theano.Variable, MultiObservedRV)):
+        if not force:
+            raise TypeError(
+                "By default, it is not allowed to add nodes that "
+                "are not `theano.Variable`'s nor pymc3 random variates "
+                "to a variable dependence graph."
+                "Got node `{}` of type `{}`. "
+                "However, this kind of node could be added by "
+                "passing `force=True` to `add`. It would be "
+                "wrapped by a `WrapAsHashable` instead.".
+                format(node, type(node)))
+        node = WrapAsHashable(node)
+    elif not (not_shared_or_constant_variable(node) or
+              hasattr(node, 'distribution')):
+        if not (force or accept_cons_shared):
+            raise ConstantNodeException(
+                'Supplied node, of type `{}`, does not have a '
+                '`distribution` attribute or is an instance of a `theano` '
+                '`Constant` or `SharedVariable`. This node could be '
+                'accepted by passing either `force=True` or '
+                '`accept_cons_shared=True` to `add`.'.
+                format(type(node)))
+    if not isinstance(node, Hashable):
+        node = WrapAsHashable(node)
+    if node in dag:
+        # Should we raise a warning with we attempt to add a node that is
+        # already in the DAG??
+        return node
+
+    # Add node into the nodes set and then initiate all node relations and
+    # values to their defaults
+    dag.add_node(node)
+
+    # Try to get the conditional parents of node and add them
+    cond = get_first_level_conditionals(node)
+    if cond is not None:
+        for conditional_parent in walk_down_ownership(cond):
+            if conditional_parent not in dag:
+                try:
+                    add_to_dependence_dag(dag, conditional_parent)
+                except ConstantNodeException:
+                    continue
+            dag.add_edge(conditional_parent,
+                         node,
+                         deterministic=0)
+
+    # Try to get the deterministic parents of node and add them
+    if not_shared_or_constant_variable(node):
+        for deterministic_parent in walk_down_ownership([node], ignore=True):
+            if deterministic_parent not in dag:
+                try:
+                    add_to_dependence_dag(dag, deterministic_parent)
+                except ConstantNodeException:
+                    continue
+            dag.add_edge(deterministic_parent,
+                         node,
+                         deterministic=1)
+    if not networkx.is_directed_acyclic_graph(dag):
+        raise RuntimeError('The dependence graph is no longer a directed '
+                           'acyclic graph (DAG). The addition of node `{}`, '
+                           'of type `{}`, and the edges from its predecessors '
+                           'introduced a loop inside the variable dependence '
+                           'graph. Consider raising an issue with developers.'.
+                           format(node, type(node)))
+
+    return node
+
+
+def walk_down_ownership(node_list, ignore=False):
+    """This function goes through an iterable of nodes provided in
+    `node_list`, yielding the named nodes in the ownership graph, whoes name is
+    not None.
+    With the optional input `ignore`, a node without a name can be yielded.
+    """
+    for node in node_list:
+        if (not ignore and hasattr(node, 'name') and node.name is not None and
+                not is_autonamed_node(node)):
+            yield node
+        elif not_shared_or_constant_variable(node):
+            owner = getattr(node, 'owner', None)
+            if owner is not None:
+                for parent in walk_down_ownership(owner.inputs):
+                    yield parent
+
+
+def get_sub_dag(dag, input_nodes, force=True):
+    """Get a new DiGraph instance which is like a right outer join
+    of `dag` with a list of input nodes provided in `input_nodes`.
+    What this means is that it will look for the `input_nodes` inside
+    `dag`, the nodes which are contained in `dag`, along with all their
+    predecessors will be used to get a subgraph of `dag`. Then, the
+    remaining nodes will be added to the `DiGraph` with
+    `add_to_dependence_dag`. Finally, this instance is returned.
+    In summary, it copies the shared part of `dag`, given the nodes in
+    `input_nodes`, and then adds onto that.
+
+    Parameters
+    ----------
+    dag : DiGraph instance
+        Dependence DAG DiGraph on which to perform the right outer join
+        operation.
+    input_nodes : list or scalar
+        If it is a scalar `input_nodes` will be converted to a list as
+        `[input_nodes]`. `input_nodes` is a list of nodes that will be
+        used to create a new `DependenceDAG` instance. The part of the DAG
+        that is shared with `self`, will be copied, and the rest will be
+        added.
+    force : bool (optional)
+        If True, the nodes that must be added, will be added with the
+        force flag set to True. [Default is True]
+
+    Returns
+    -------
+    A tuple. The first element is a new `DiGraph` instance that results from
+    the right outer join operation.
+    The second element in the output tuple is the dictionary of indices to
+    nodes `{index: node}`. Each key will be the position of the node provided
+    in `input_nodes` and the value will be the node that was added to the
+    returned `DiGraph` instance, which could either be a
+    `WrapAsHashable` instance or `input_nodes[index]` itself.
+    """
+    if not isinstance(dag, networkx.DiGraph):
+        raise TypeError('Input `dag` must be an instance of networkx.DiGraph. '
+                        'Got {}, which has type {}, instead.'.
+                        format(dag, type(dag)))
+    if not isinstance(input_nodes, list):
+        input_nodes = [input_nodes]
+    index = {}
+    copied_bunch = set()
+    nodes_to_add = []
+    for i, node in enumerate(input_nodes):
+        if node in dag:
+            index[i] = node
+            copied_bunch.add(node)
+            copied_bunch.update(networkx.ancestors(dag, node))
+        else:
+            nodes_to_add.append((i, node))
+
+    subgraph = dag.subgraph(copied_bunch).copy()
+    for node_index, node in nodes_to_add:
+        added_node = add_to_dependence_dag(subgraph,
+                                           node,
+                                           force=force)
+        if node_index is not None:
+            index[node_index] = added_node
+    return subgraph, index
+
+
+def matching_dependence_dags(a, b):
+    """A helper function intended to be used during debugging. Is true is
+    two DiGraph instances that represent a dependence DAG contain the same
+    nodes, the same edges and the `deterministic` attribute of all edges are
+    the same.
+    """
+    if (not isinstance(a, networkx.DiGraph) or
+            not isinstance(a, networkx.DiGraph)):
+        return False
+    if len(a.nodes) != len(b.nodes):
+        return False
+    for n in b.nodes:
+        if n not in a:
+            return False
+        for edge in b.edges(n, data=True):
+            if not a.has_edge(edge[0], edge[1]):
+                return False
+            data = a[edge[0]][edge[1]]
+            if data['deterministic'] != edge[2]['deterministic']:
+                return False
+    return True
+
+
+def build_dependence_dag_from_model(model):
+    dag = networkx.DiGraph()
+    for node in model.free_RVs + model.deterministics + model.observed_RVs:
+        add_to_dependence_dag(dag, node, accept_cons_shared=True)
+    return dag
diff --git a/pymc3/model_graph.py b/pymc3/model_graph.py
index e68a28cb17..4e0f81554a 100644
--- a/pymc3/model_graph.py
+++ b/pymc3/model_graph.py
@@ -4,6 +4,8 @@
 
 from .util import get_default_varnames
 import pymc3 as pm
+from .model import build_dependence_dag_from_model
+import networkx as nx
 
 
 def powerset(iterable):
@@ -184,3 +186,175 @@ def model_to_graphviz(model=None):
     """
     model = pm.modelcontext(model)
     return ModelGraph(model).make_graph()
+
+
+class OtherModelGraph(object):
+    def __init__(self, model):
+        self.model = model
+        try:
+            graph = model.dependence_dag
+        except AttributeError:
+            graph = build_dependence_dag_from_model(model)
+        self.set_graph(graph)
+
+    def set_graph(self, graph):
+        self.graph = graph
+        self.node_names = {}
+        unnamed_count = 0
+        for n in self.graph.nodes():
+            try:
+                name = n.name
+            except AttributeError:
+                name = 'Unnamed {}'.format(unnamed_count)
+                unnamed_count += 1
+            self.node_names[n] = name
+
+    def draw(self, pos=None, draw_nodes=False, ax=None,
+             edge_kwargs=None,
+             label_kwargs=None,
+             node_kwargs=None):
+        graph = self.graph
+        if edge_kwargs is None:
+            edge_kwargs = {}
+        if node_kwargs is None:
+            node_kwargs = {}
+        if label_kwargs is None:
+            label_kwargs = {}
+        label_kwargs.setdefault('bbox', {'boxstyle': 'round',
+                                         'facecolor': 'lightgray'})
+        if pos is None:
+            try:
+                pos = nx.drawing.nx_agraph.graphviz_layout(graph, prog='dot')
+            except Exception:
+                pos = nx.shell_layout(graph)
+        d = nx.get_edge_attributes(graph, 'deterministic')
+        edgelist = list(d.keys())
+        edge_color = [float(v) for v in d.values()]
+        labels = {n: n.name for n in graph}
+
+        if draw_nodes:
+            nx.draw_networkx_edges(graph, pos=pos, ax=ax, **node_kwargs)
+        nx.draw_networkx_edges(graph, pos=pos, ax=ax, edgelist=edgelist,
+                               edge_color=edge_color, **edge_kwargs)
+        nx.draw_networkx_labels(graph, pos=pos, ax=ax, labels=labels,
+                                **label_kwargs)
+
+    def get_plates(self, graph=None, ignore_transformed=True):
+        """ Groups nodes by the shape of the underlying distribution, and if
+        the nodes form a disconnected component of the graph.
+
+        Parameters
+        ----------
+        graph: networkx.DiGraph (optional)
+            The graph object from which to get the plates. If None, self.graph
+            will be used.
+        ignore_transformed: bool (optional)
+            If True, the transformed variables will be ignored while getting
+            the plates.
+
+        Returns
+        -------
+        list of tuples: (shape, set(nodes_in_plate))
+        """
+        if graph is None:
+            graph = self.graph
+        if ignore_transformed:
+            transforms = set([n.transformed for n in graph
+                              if hasattr(n, 'transformed')])
+            nbunch = [n for n in graph if n not in transforms]
+            graph = nx.subgraph(graph, nbunch)
+        shape_plates = {}
+        for node in graph:
+            if hasattr(node, 'observations'):
+                shape = node.observations.shape
+            elif hasattr(node, 'dshape'):
+                shape = node.dshape
+            else:
+                try:
+                    shape = node.tag.test_value.shape
+                except AttributeError:
+                    shape = tuple()
+            if shape == (1,):
+                shape = tuple()
+            if shape not in shape_plates:
+                shape_plates[shape] = set()
+            shape_plates[shape].add(node)
+        plates = []
+        for shape, nodes in shape_plates.items():
+            # We want to find the disconnected components that have a common
+            # shape. These will be the plates
+            subgraph = nx.subgraph(graph, nodes).to_undirected()
+            for G in nx.connected_component_subgraphs(subgraph, copy=False):
+                plates.append((shape, set(G.nodes())))
+        return plates
+
+    def make_graph(self, ignore_transformed=True, edge_cmap=None):
+        """Make graphviz Digraph of PyMC3 model
+
+        Returns
+        -------
+        graphviz.Digraph
+        """
+        try:
+            import graphviz
+        except ImportError:
+            raise ImportError('This function requires the python library graphviz, along with binaries. '
+                              'The easiest way to install all of this is by running\n\n'
+                              '\tconda install -c conda-forge python-graphviz')
+
+        G = self.graph
+        if ignore_transformed:
+            transforms = set([n.transformed for n in G
+                              if hasattr(n, 'transformed')])
+            nbunch = [n for n in G if n not in transforms]
+            G = nx.subgraph(G, nbunch)
+        graph = graphviz.Digraph(self.model.name)
+        nclusters = 0
+        for shape, nodes in self.get_plates(graph=G):
+            label = ' x '.join(map('{:,d}'.format, shape))
+            if label:
+                # must be preceded by 'cluster' to get a box around it
+                with graph.subgraph(name='cluster {}'.format(nclusters)) as sub:
+                    nclusters += 1
+                    for node in nodes:
+                        self._make_node(node, sub)
+                    # plate label goes bottom right
+                    sub.attr(label=label, labeljust='r', labelloc='b', style='rounded')
+            else:
+                for node in nodes:
+                    self._make_node(node, graph)
+
+        for from_node, to_node, ats in G.edges(data=True):
+            if edge_cmap is None:
+                edge_color = '#000000'
+            else:
+                from matplotlib import colors
+                val = float(ats['deterministic'])
+                edge_color = colors.to_hex(edge_cmap(val), keep_alpha=True)
+            graph.edge(self.node_names[from_node],
+                       self.node_names[to_node],
+                       color=edge_color)
+        return graph
+
+    def _make_node(self, node, graph):
+        """Attaches the given variable to a graphviz Digraph"""
+        # styling for node
+        attrs = {}
+        if isinstance(node, pm.model.ObservedRV):
+            attrs['style'] = 'filled'
+
+        # Get name for node
+        if hasattr(node, 'distribution'):
+            distribution = node.distribution.__class__.__name__
+        else:
+            distribution = 'Deterministic'
+            attrs['shape'] = 'box'
+        var_name = self.node_names[node]
+
+        graph.node(var_name,
+                '{var_name} ~ {distribution}'.format(var_name=var_name, distribution=distribution),
+                **attrs)
+
+
+def crude_draw(model, *args, **kwargs):
+    OtherModelGraph(model).draw(*args, **kwargs)
diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index af4e36a7f2..9b5defd761 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -1128,7 +1128,10 @@ def sample_posterior_predictive(trace, samples=None, model=None, vars=None, size
 
     # draw once to inspect the shape
     var_values = list(zip(varnames,
-                          draw_values(vars, point=model.test_point, size=size)))
+                          draw_values(vars,
+                                      point=model.test_point,
+                                      size=size,
+                                      model=model)))
     ppc_trace = defaultdict(list)
     for varname, value in var_values:
         ppc_trace[varname] = np.zeros((samples,) + value.shape, value.dtype)
@@ -1141,7 +1144,7 @@ def sample_posterior_predictive(trace, samples=None, model=None, vars=None, size
             else:
                 param = trace[idx % len_trace]
 
-            values = draw_values(vars, point=param, size=size)
+            values = draw_values(vars, point=param, size=size, model=model)
             for k, v in zip(vars, values):
                 ppc_trace[k.name][slc] = v
 
@@ -1279,6 +1282,7 @@ def sample_posterior_predictive_w(traces, samples=None, models=None, weights=Non
             var = variables[idx]
             # TODO sample_posterior_predictive_w is currently only work for model with
             # one observed.
+            # TODO supply the proper model to draw_values
             ppc[var.name].append(draw_values([var],
                                              point=param,
                                              size=size[idx]
@@ -1331,7 +1335,8 @@ def sample_prior_predictive(samples=500, model=None, vars=None, random_seed=None
         np.random.seed(random_seed)
     names = get_default_varnames(model.named_vars, include_transformed=False)
     # draw_values fails with auto-transformed variables. transform them later!
-    values = draw_values([model[name] for name in names], size=samples)
+    values = draw_values([model[name] for name in names], size=samples,
+                         model=model)
 
     data = {k: v for k, v in zip(names, values)}
 
diff --git a/pymc3/tests/test_model_graph.py b/pymc3/tests/test_model_graph.py
index c261410cc9..72e6b5c2c1 100644
--- a/pymc3/tests/test_model_graph.py
+++ b/pymc3/tests/test_model_graph.py
@@ -1,6 +1,7 @@
 import numpy as np
 import pymc3 as pm
-from pymc3.model_graph import ModelGraph, model_to_graphviz
+from pymc3.model_graph import ModelGraph, model_to_graphviz, OtherModelGraph
+from theano import tensor as tt
 
 from .helpers import SeededTest
 
@@ -77,3 +78,90 @@ def test_graphviz(self):
         for key in self.compute_graph:
             assert key in g.source
 
+
+def radon_model2():
+    """Similar in shape to the Radon model"""
+    n_homes = 919
+    counties = 85
+    uranium = np.random.normal(-.1, 0.4, size=n_homes)
+    xbar = np.random.normal(1, 0.1, size=n_homes)
+    floor_measure = np.random.randint(0, 2, size=n_homes)
+    log_radon = np.random.normal(1, 1, size=n_homes)
+    radon_dispensers = np.random.normal(1, 0.2, size=n_homes)
+
+    d, r = divmod(919, 85)
+    county = np.hstack((
+        np.tile(np.arange(counties, dtype=int), d),
+        np.arange(r)
+        ))
+    with pm.Model() as model:
+        sigma_a = pm.HalfCauchy('sigma_a', 5)
+        gamma = pm.Normal('gamma', mu=0., sd=1e5, shape=3)
+        mu_a = pm.Deterministic('mu_a', gamma[0] + gamma[1]*uranium + gamma[2]*xbar)
+        eps_a = pm.Normal('eps_a', mu=0, sd=sigma_a, shape=counties)
+        a = pm.Deterministic('a', mu_a + eps_a[county])
+        b = pm.Normal('b', mu=0., sd=1e15)
+        sigma_y = pm.Uniform('sigma_y', lower=0, upper=100)
+        p_rad = pm.Beta('p_rad', alpha=1, beta=1)
+        disp_on = pm.Bernoulli('disp_on', p=p_rad, shape=n_homes)
+        rad_cloud = pm.Deterministic('rad_cloud', tt.mean(disp_on * radon_dispensers))
+        y_hat = a + b * floor_measure + rad_cloud
+        y_like = pm.Normal('y_like', mu=y_hat, sd=sigma_y, observed=log_radon)
+
+    plates_without_trans = [((3,), {gamma}),
+                            ((), {b}),
+                            ((), {rad_cloud}),
+                            ((), {sigma_a}),
+                            ((), {p_rad}),
+                            ((), {sigma_y}),
+                            ((919,), {a, mu_a, y_like}),
+                            ((919,), {disp_on}),
+                            ((85,), {eps_a})]
+    plates_with_trans = [((), {b}),
+                         ((), {rad_cloud}),
+                         ((), {sigma_a, model['sigma_a_log__']}),
+                         ((), {sigma_y, model['sigma_y_interval__']}),
+                         ((), {p_rad, model['p_rad_logodds__']}),
+                         ((3,), {gamma}),
+                         ((919,), {a, mu_a, y_like}),
+                         ((919,), {disp_on}),
+                         ((85,), {eps_a})]
+    return model, plates_without_trans, plates_with_trans
+
+
+class TestSimpleModel2(SeededTest):
+    @classmethod
+    def setup_class(cls):
+        cls.model, cls.plates_without_trans, cls.plates_with_trans= radon_model2()
+        cls.model_graph = OtherModelGraph(cls.model)
+
+    def test_plates(self):
+        assert all([plate in self.plates_without_trans
+                    for plate in self.model_graph.get_plates(ignore_transformed=True)])
+        assert all([plate in self.plates_with_trans
+                    for plate in self.model_graph.get_plates(ignore_transformed=False)])
+
+    def test_graphviz(self):
+        # just make sure everything runs without error
+
+        transforms = set([n.transformed for n in self.model_graph.graph
+                          if hasattr(n, 'transformed')])
+        g = self.model_graph.make_graph(ignore_transformed=True)
+        for node, key in self.model_graph.node_names.items():
+            if node not in transforms:
+                assert key in g.source
+            else:
+                assert key not in g.source
+        g = OtherModelGraph(self.model).make_graph(ignore_transformed=True)
+        for node, key in self.model_graph.node_names.items():
+            if node not in transforms:
+                assert key in g.source
+            else:
+                assert key not in g.source
+
+        g = self.model_graph.make_graph(ignore_transformed=False)
+        for node, key in self.model_graph.node_names.items():
+            assert key in g.source
+        g = OtherModelGraph(self.model).make_graph(ignore_transformed=False)
+        for node, key in self.model_graph.node_names.items():
+            assert key in g.source
diff --git a/pymc3/tests/test_model_helpers.py b/pymc3/tests/test_model_helpers.py
index 2b191bd144..b75b47bb31 100644
--- a/pymc3/tests/test_model_helpers.py
+++ b/pymc3/tests/test_model_helpers.py
@@ -2,7 +2,13 @@
 import numpy.ma as ma
 import numpy.testing as npt
 import pandas as pd
+from networkx import DiGraph, is_directed_acyclic_graph, topological_sort
 import pymc3 as pm
+from pymc3.model import (build_dependence_dag_from_model,
+                         matching_dependence_dags,
+                         add_to_dependence_dag,
+                         get_sub_dag)
+from pymc3.util import WrapAsHashable
 import scipy.sparse as sps
 
 import theano
@@ -138,3 +144,100 @@ def test_as_tensor(self):
         assert masked_output.missing_values is not None
 
         return None
+
+
+class TestDependenceDAG(object):
+    def setup_method(self):
+        self.obs = np.random.randn(1000,) + 2
+        with pm.Model() as model:
+            self.a = pm.Normal('a', mu=0, sd=100)
+            self.b = pm.Normal('b', mu=self.a, sd=1e-8)
+            self.c = pm.Normal('c', mu=self.a, sd=1e-8)
+            self.d = pm.Deterministic('d', self.b + self.c)
+            self.e = pm.Normal('e', mu=self.d, sd=1, observed=self.obs)
+        self.model = model
+        self.expected_full_dag = DiGraph()
+        add_to_dependence_dag(self.expected_full_dag, self.e)
+
+    def test_built_DependenceDAG(self):
+        assert matching_dependence_dags(self.expected_full_dag,
+                                        self.model.dependence_dag)
+        assert matching_dependence_dags(
+            build_dependence_dag_from_model(self.model),
+            self.model.dependence_dag)
+
+    def test_get_sub_dag(self):
+        dag = self.model.dependence_dag
+        sub1 = get_sub_dag(dag, self.a)[0]
+        assert len(sub1.nodes) == 1
+        assert is_directed_acyclic_graph(sub1)
+
+        sub2 = get_sub_dag(dag, [self.a])[0]
+        assert len(sub2.nodes) == 1
+        assert is_directed_acyclic_graph(sub2)
+        assert matching_dependence_dags(sub1, sub2)
+
+        sub3 = get_sub_dag(dag, [self.e])[0]
+        assert len(sub3.nodes) == 5
+        assert is_directed_acyclic_graph(sub3)
+        assert matching_dependence_dags(sub3, dag)
+
+        hard_expr = (theano.tensor.exp(self.b + self.e * self.e) *
+                     self.e * self.b + self.a)
+        hard = get_sub_dag(dag, [self.e, hard_expr])[0]
+        assert len(hard.nodes) == 6
+        assert is_directed_acyclic_graph(hard)
+        sorted_nodes = list(topological_sort(hard))
+        expected = [(self.a,),
+                    (self.b, self.c),
+                    (self.b, self.c),
+                    (self.d,),
+                    (self.e,),
+                    (hard_expr,)]
+        assert all((n in e for n, e in zip(sorted_nodes, expected)))
+        assert matching_dependence_dags(get_sub_dag(hard, self.e)[0], dag)
+
+        params = [self.d,
+                  0.,
+                  np.zeros((10, 2), dtype=np.float32),
+                  theano.tensor.constant(5.354),
+                  theano.shared(np.array([2, 6, 8])),
+                  ]
+        with_non_theano, index = get_sub_dag(dag,
+                                             params)
+        assert is_directed_acyclic_graph(with_non_theano)
+        assert matching_dependence_dags(get_sub_dag(with_non_theano,
+                                                    self.d)[0],
+                                        get_sub_dag(dag, [self.d])[0])
+        assert index[0] == params[0]
+        assert all([isinstance(index[i], WrapAsHashable)
+                    for i in range(1, 3)])
+        for i in range(len(params)):
+            supplied = params[i]
+            obj = index[i]
+            if not isinstance(obj, WrapAsHashable):
+                continue
+            assert obj in with_non_theano
+            if obj.node_is_hashable:
+                assert obj.node == supplied
+            else:
+                assert id(obj.node) == id(params[i])
+            obj_value = obj.get_value()
+            if isinstance(supplied, theano.tensor.sharedvar.SharedVariable):
+                expected_value = supplied.get_value()
+            elif isinstance(supplied, theano.tensor.TensorConstant):
+                expected_value = supplied.value
+            else:
+                expected_value = supplied
+            if isinstance(obj_value, np.ndarray):
+                assert np.all(obj_value == expected_value)
+            else:
+                assert obj_value == expected_value
+
+        wnt2 = with_non_theano.copy()
+        assert matching_dependence_dags(with_non_theano, wnt2)
+        node2 = add_to_dependence_dag(wnt2,
+                                      params[-1],
+                                      force=True)
+        assert matching_dependence_dags(wnt2, with_non_theano)
+        assert node2 == index[len(params) - 1]
diff --git a/pymc3/tests/test_random.py b/pymc3/tests/test_random.py
index f61fe3b4a7..2aa74f6fec 100644
--- a/pymc3/tests/test_random.py
+++ b/pymc3/tests/test_random.py
@@ -1,10 +1,11 @@
 import pymc3 as pm
 import numpy as np
+import numpy.random as nr
 import numpy.testing as npt
 import pytest
 import theano.tensor as tt
 import theano
-
+from .helpers import SeededTest
 from pymc3.distributions.distribution import _draw_value, draw_values
 
 
@@ -60,7 +61,7 @@ def test_vals(self):
 
     def test_simple_model(self):
         with pm.Model():
-            mu = 2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5))
+            mu = (2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5)))
             a = pm.Normal('a', mu=mu, sd=5, shape=2)
 
         val1 = draw_values([a])
@@ -72,7 +73,7 @@ def test_simple_model(self):
 
     def test_dep_vars(self):
         with pm.Model():
-            mu = 2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5))
+            mu = (2 * tt.constant(np.array([5., 6.])) + theano.shared(np.array(5)))
             sd = pm.HalfNormal('sd', shape=2)
             tau = 1 / sd ** 2
             a = pm.Normal('a', mu=mu, tau=tau, shape=2)
@@ -88,3 +89,59 @@ def test_dep_vars(self):
         assert all([np.all(val1 != val2), np.all(val1 != val3),
                     np.all(val1 != val4), np.all(val2 != val3),
                     np.all(val2 != val4), np.all(val3 != val4)])
+
+
+class TestJointDistributionDrawValues(SeededTest):
+    def test_joint_distribution(self):
+        with pm.Model() as model:
+            a = pm.Normal('a', mu=0, sd=100)
+            b = pm.Normal('b', mu=a, sd=1e-8)
+            c = pm.Normal('c', mu=a, sd=1e-8)
+            d = pm.Deterministic('d', b + c)
+
+        # Expected RVs
+        N = 1000
+        norm = np.random.randn(3, N)
+        eA = norm[0] * 100
+        eB = eA + norm[1] * 1e-8
+        eC = eA + norm[2] * 1e-8
+        eD = eB + eC
+
+        # Drawn RVs
+        nr.seed(self.random_seed)
+#        A, B, C, D = list(zip(*[draw_values([a, b, c, d]) for i in range(N)]))
+        A, B, C, D = draw_values([a, b, c, d], size=N)
+        A = np.array(A).flatten()
+        B = np.array(B).flatten()
+        C = np.array(C).flatten()
+        D = np.array(D).flatten()
+
+        # Assert that the drawn samples match the expected values
+        assert np.allclose(eA, A)
+        assert np.allclose(eB, B)
+        assert np.allclose(eC, C)
+        assert np.allclose(eD, D)
+
+        # Assert that A, B and C have the expected difference
+        assert np.all(np.abs(A - B) < 1e-6)
+        assert np.all(np.abs(A - C) < 1e-6)
+        assert np.all(np.abs(B - C) < 1e-6)
+
+        # Marginal draws
+        mA = np.array([draw_values([a]) for i in range(N)]).flatten()
+        mB = np.array([draw_values([b]) for i in range(N)]).flatten()
+        mC = np.array([draw_values([c]) for i in range(N)]).flatten()
+        # Also test the with model context of draw_values
+        with model:
+            mD = np.array([draw_values([d]) for i in range(N)]).flatten()
+
+        # Assert that the marginal distributions have different sample values
+        assert not np.all(np.abs(B - mB) < 1e-2)
+        assert not np.all(np.abs(C - mC) < 1e-2)
+        assert not np.all(np.abs(D - mD) < 1e-2)
+
+        # Assert that the marginal distributions do not have high cross
+        # correlation
+        assert np.abs(np.corrcoef(mA, mB)[0, 1]) < 0.1
+        assert np.abs(np.corrcoef(mA, mC)[0, 1]) < 0.1
+        assert np.abs(np.corrcoef(mB, mC)[0, 1]) < 0.1
diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py
index 4b1cfe17b8..28b1814dcb 100644
--- a/pymc3/tests/test_step.py
+++ b/pymc3/tests/test_step.py
@@ -467,206 +467,206 @@ class TestStepMethods(object):  # yield test doesn't work subclassing object
         ),
         SMC: np.array(
             [
-                0.02957155,
-                -0.10448987,
-                1.83330321,
-                0.45125892,
-                0.08847866,
-                1.42446262,
-                -0.50373428,
-                1.7926401,
-                -0.43595123,
-                1.29051484,
-                1.48825784,
-                0.58940796,
-                0.39452458,
-                0.6458593,
-                1.41254985,
-                0.33640003,
-                -0.82261885,
-                0.19660577,
-                -1.31652677,
-                -0.43394062,
-                0.47844794,
-                0.32073809,
-                0.43811662,
-                0.15053655,
-                0.36644253,
-                0.67791414,
-                -0.02508511,
-                0.7298907,
-                1.98694909,
-                1.7186239,
-                0.36670138,
-                1.31585297,
-                1.0010273,
-                1.03616273,
-                0.14999443,
-                0.27743457,
-                -0.01086402,
-                1.49293981,
-                1.56440459,
-                -0.05856325,
-                1.21029567,
-                0.3526936,
-                0.72190671,
-                0.06857266,
-                2.5848243,
-                0.03416341,
-                0.11731779,
-                0.71426943,
-                -0.01581274,
-                0.31147419,
-                0.59718523,
-                0.25670663,
-                0.71827237,
-                0.19043905,
-                1.32676133,
-                -0.02546965,
-                1.56694401,
-                -0.27818479,
-                -0.32227988,
-                -0.58671167,
-                0.18231422,
-                -0.51437682,
-                0.60519212,
-                2.01615164,
-                1.25931896,
-                1.0957265,
-                2.35609054,
-                -1.03617366,
-                -0.38976571,
-                1.2889846,
-                0.39140397,
-                -0.88146889,
-                0.11897506,
-                0.68374636,
-                0.35969604,
-                1.55791486,
-                0.60414685,
-                0.76464213,
-                0.5734933,
-                0.55888177,
-                1.14165163,
-                0.29172266,
-                1.0656623,
-                1.49059082,
-                -0.78259634,
-                1.39116132,
-                -1.30272835,
-                1.47197121,
-                -0.44181174,
-                0.61235063,
-                0.06672622,
-                0.15037854,
-                1.3167083,
-                1.24502001,
-                0.42421749,
-                -0.31524216,
-                0.0315209,
-                0.52553104,
-                -0.12301504,
-                -1.01849236,
-                1.66776635,
-                0.98537117,
-                -0.52051714,
-                1.15968179,
-                -0.41569573,
-                -0.30220252,
-                -0.24221614,
-                1.01834671,
-                0.79537263,
-                0.47201206,
-                0.53230357,
-                0.95621545,
-                0.67045256,
-                0.56275657,
-                0.66782958,
-                1.04360486,
-                1.48513602,
-                1.3954355,
-                0.72727048,
-                0.65556693,
-                -0.80579428,
-                -0.6951795,
-                -0.41933972,
-                0.79093036,
-                -0.13152736,
-                0.19672923,
-                0.61027041,
-                0.34030084,
-                1.78314475,
-                0.99044934,
-                0.51198086,
-                -0.2394725,
-                -0.50457759,
-                0.57011011,
-                1.47181142,
-                0.73780165,
-                1.67648809,
-                0.48628167,
-                1.12112297,
-                0.49702466,
-                0.67702156,
-                1.09022274,
-                1.06579346,
-                0.80085527,
-                0.11645159,
-                -1.58625901,
-                1.60302073,
-                0.85477186,
-                -0.06355492,
-                0.94015912,
-                0.81545671,
-                1.36875334,
-                -0.09772053,
-                -0.39652235,
-                0.11642491,
-                -0.13081192,
-                0.38424794,
-                0.07590825,
-                0.04093988,
-                1.41189211,
-                1.21406156,
-                0.73694675,
-                0.0660595,
-                0.25210641,
-                0.06578272,
-                0.26790674,
-                0.77504874,
-                -0.7038388,
-                -0.32709336,
-                1.25433616,
-                -0.01871021,
-                -0.064151,
-                0.7459666,
-                -0.28951973,
-                -0.8878012,
-                -0.43049727,
-                -0.30480856,
-                2.2053675,
-                0.1190498,
-                0.21880221,
-                0.77411644,
-                1.22875531,
-                0.63578569,
-                0.96799128,
-                0.50535637,
-                -0.68404005,
-                0.66134393,
-                0.60923846,
-                -0.39533348,
-                0.70696449,
-                0.14433537,
-                0.06392199,
-                0.53761533,
-                0.58593561,
-                0.45146962,
-                1.16348982,
-                0.9917836,
-                -0.35454457,
-                0.52929041,
-                0.83350277,
+                0.53154128,
+                1.69014751,
+                0.38863621,
+                1.36550742,
+                0.54937705,
+                0.85452502,
+                1.34000193,
+                0.04963276,
+                0.73207585,
+                -0.45504452,
+                0.99087969,
+                0.53800418,
+                1.69481083,
+                0.19015456,
+                0.54587521,
+                0.51198155,
+                0.17514563,
+                -0.62023686,
+                0.73211032,
+                1.0269751,
+                0.82004582,
+                1.07714066,
+                1.27243655,
+                0.8603388,
+                -0.96709536,
+                -0.4963755,
+                0.47436677,
+                0.34392296,
+                0.08501226,
+                0.95779747,
+                1.21125461,
+                -0.04609757,
+                0.29714065,
+                0.89447118,
+                0.00472546,
+                0.50365803,
+                1.73127064,
+                1.04164544,
+                -0.22236077,
+                1.33631993,
+                0.96357527,
+                1.06122196,
+                0.12798557,
+                0.4665272,
+                -0.1162582,
+                1.62002463,
+                1.44557222,
+                -0.49218985,
+                1.2175313,
+                0.25761981,
+                0.82879531,
+                0.16321047,
+                1.34260731,
+                -0.05709803,
+                0.18903618,
+                0.76825821,
+                0.08211472,
+                0.53817434,
+                0.53379232,
+                -0.47094362,
+                1.14433914,
+                0.03770296,
+                1.30737805,
+                0.39671022,
+                1.22440728,
+                0.09600212,
+                -0.49796137,
+                -0.44963869,
+                0.95208986,
+                -0.04308567,
+                0.45937807,
+                2.59887219,
+                0.36326674,
+                1.16659465,
+                2.26888158,
+                -0.64081701,
+                0.13085995,
+                1.5847621,
+                0.29342994,
+                -0.7802778,
+                0.62631291,
+                0.56155063,
+                0.63017152,
+                1.88801376,
+                0.32864795,
+                0.19722366,
+                0.62506725,
+                -0.04154236,
+                0.74923865,
+                0.64958051,
+                1.05205509,
+                1.12818507,
+                0.35463018,
+                1.49394093,
+                -1.32280176,
+                0.48922758,
+                -0.67185953,
+                0.01282045,
+                -0.00832875,
+                0.60746178,
+                1.04869967,
+                0.43197615,
+                0.14665959,
+                -0.08117829,
+                0.43216574,
+                0.87241428,
+                -0.07985268,
+                -0.93380876,
+                1.73662159,
+                0.23926283,
+                -0.69068641,
+                1.17829179,
+                -0.16332134,
+                -0.5112194,
+                -0.43442261,
+                0.34948852,
+                1.11002685,
+                0.42445302,
+                0.68379355,
+                -0.12877628,
+                0.59561974,
+                0.67230016,
+                1.67895815,
+                1.51053172,
+                1.14415702,
+                1.00682996,
+                1.09882483,
+                0.28820149,
+                -0.75250142,
+                -0.66453929,
+                -0.0991988,
+                0.2907921,
+                0.04077525,
+                -0.44036405,
+                0.44894708,
+                0.68646345,
+                0.03986746,
+                0.50061203,
+                1.18904715,
+                0.36231722,
+                -0.16347099,
+                0.35343108,
+                1.15870795,
+                0.5973369,
+                1.50731862,
+                0.69983246,
+                1.50854043,
+                0.97489667,
+                0.25267479,
+                0.26369507,
+                1.59775053,
+                1.56383915,
+                0.1721522,
+                -0.96935772,
+                1.47191825,
+                0.79858327,
+                0.69071774,
+                -0.17667758,
+                0.61438524,
+                0.99424152,
+                -0.23558854,
+                -0.27873225,
+                0.16615446,
+                0.02589567,
+                -0.38007309,
+                0.24960815,
+                1.17127086,
+                1.96577002,
+                0.83224965,
+                0.9386008,
+                -0.16018964,
+                0.25239747,
+                -0.09936852,
+                -0.20376629,
+                1.39291948,
+                -0.2083358,
+                0.51435573,
+                1.38304537,
+                0.23272616,
+                -0.15257289,
+                0.77293718,
+                0.33558962,
+                -0.99534345,
+                -0.03472376,
+                0.07169895,
+                1.62726823,
+                0.08074445,
+                0.38765492,
+                0.7844363,
+                0.89340893,
+                0.28605836,
+                0.83632054,
+                0.54210362,
+                -0.55168005,
+                0.91756515,
+                0.16982656,
+                -0.36404392,
+                1.0011945,
+                -0.2659181,
+                0.31691263,
             ]
         ),
     }
diff --git a/pymc3/tests/test_util.py b/pymc3/tests/test_util.py
index 52e567cc75..46031933e0 100644
--- a/pymc3/tests/test_util.py
+++ b/pymc3/tests/test_util.py
@@ -1,10 +1,16 @@
 import pytest
-
+import numbers
+try:
+    from collections.abc import Hashable
+except ImportError:
+    from collections import Hashable
 import pymc3 as pm
 import numpy as np
 from numpy.testing import assert_almost_equal
 from .helpers import SeededTest
 from pymc3.distributions.transforms import Transform
+from pymc3.util import WrapAsHashable
+import theano
 
 
 class TestTransformName(object):
@@ -84,3 +90,54 @@ def test_soft_update_parent(self):
                             test_point['interv_interval__'])
 
 
+class TestWrapAsHashable(object):
+    def setup_method(self):
+        class test_obj(object):
+            def get_value(self):
+                return 'test_obj_value'
+        self.objects = ['string',
+                        123,
+                        -123.543,
+                        np.arange(10),
+                        {'key': 'value'},
+                        [0, 1, 2],
+                        theano.tensor.constant(10),
+                        theano.shared(np.arange(5)),
+                        test_obj()]
+        self.wraps = [WrapAsHashable(obj) for obj in self.objects]
+
+    def test_hashable(self):
+        for obj, wrap in zip(self.objects, self.wraps):
+            assert isinstance(wrap, Hashable)
+            if isinstance(obj, Hashable):
+                assert hash(obj) == hash(wrap)
+            else:
+                assert hash((type(obj), id(obj))) == hash(wrap)
+
+    def test_eq(self):
+        for obj, wrap in zip(self.objects, self.wraps):
+            wrap2 = WrapAsHashable(obj)
+            assert wrap2 == wrap
+            assert wrap == wrap2
+            assert wrap != obj
+
+    def test_get_value(self):
+        for obj, wrap in zip(self.objects, self.wraps):
+            if isinstance(obj, (numbers.Number,
+                                np.ndarray,
+                                theano.tensor.TensorConstant,
+                                theano.tensor.sharedvar.SharedVariable)):
+                wrap_value = wrap.get_value()
+                if isinstance(obj, numbers.Number):
+                    assert obj == wrap_value
+                elif isinstance(obj, np.ndarray):
+                    assert np.all(obj == wrap_value)
+                elif isinstance(obj, theano.tensor.TensorConstant):
+                    assert np.all(obj.value == wrap_value)
+                elif isinstance(obj, theano.tensor.sharedvar.SharedVariable):
+                    assert np.all(obj.get_value() == wrap_value)
+            elif hasattr(obj, 'get_value'):
+                assert np.all(obj.get_value() == wrap.get_value())
+            else:
+                with pytest.raises(TypeError):
+                    wrap.get_value()
diff --git a/pymc3/util.py b/pymc3/util.py
index 080a2121ed..0e482f4098 100644
--- a/pymc3/util.py
+++ b/pymc3/util.py
@@ -1,6 +1,13 @@
 import re
+import numbers
+try:
+    from collections.abc import Hashable
+except ImportError:
+    from collections import Hashable
 import functools
+import numpy as np
 from numpy import asscalar
+import theano
 
 LATEX_ESCAPE_RE = re.compile(r'(%|_|\$|#|&)', re.MULTILINE)
 
@@ -166,3 +173,58 @@ def enhanced(*args, **kwargs):
             newwrapper = functools.partial(wrapper, *args, **kwargs)
             return newwrapper
     return enhanced
+
+
+class WrapAsHashable(object):
+    """Generic class that can provide a `__hash__` and `__eq__` method to any
+    object, which will be referred to as node. It also provides a multipurpose
+    `get_value` function to get the value of the wrapped node.
+
+    If a node is hashable, then __hash__ and __eq__ will both return the same
+    as the node's __hash__ and __eq__. If the node is unhashable, the __eq__
+    will compare the node's id with the other object's id, and the __hash__
+    will be made from a hash of (type(node), id(node)).
+
+    """
+    def __init__(self, node):
+        self.node_is_hashable = isinstance(node, Hashable)
+        self.node = node
+
+    def __hash__(self):
+        if self.node_is_hashable:
+            return self.node.__hash__()
+        else:
+            return hash((type(self.node), id(self.node)))
+
+    def __eq__(self, other):
+        if isinstance(other, WrapAsHashable):
+            if self.node_is_hashable and other.node_is_hashable:
+                return self.node == other.node
+            elif not self.node_is_hashable and not other.node_is_hashable:
+                return id(self.node) == id(other.node)
+            else:
+                return False
+        else:
+            return False
+
+    def __ne__(self, other):
+        return not self == other
+
+    def get_value(self):
+        node = self.node
+        if isinstance(node, numbers.Number):
+            return node
+        elif isinstance(node, np.ndarray):
+            return node
+        elif isinstance(node, theano.tensor.TensorConstant):
+            return node.value
+        elif isinstance(node, theano.tensor.sharedvar.SharedVariable):
+            return node.get_value()
+        else:
+            try:
+                return node.get_value()
+            except AttributeError:
+                raise TypeError("WrapAsHashable instance's node, {}, does not "
+                                "define a `get_value` method and is not of "
+                                "known types. type(node) = {}.".
+                                format(node, type(node)))
diff --git a/requirements-dev.txt b/requirements-dev.txt
index d8cfdee554..ac10e22a40 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -22,3 +22,4 @@ scipy>=0.12.0
 seaborn>=0.8.1
 sphinx-autobuild==0.7.1
 sphinx>=1.5.5
+networkx>=1.11.0
diff --git a/requirements.txt b/requirements.txt
index 3684a701e5..f5370f913d 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -8,3 +8,4 @@ tqdm>=4.8.4
 six>=1.10.0
 h5py>=2.7.0
 enum34>=1.1.6; python_version < '3.4'
+networkx>=1.11.0
diff --git a/scripts/create_testenv.sh b/scripts/create_testenv.sh
index fe1c4cc249..49e272deca 100755
--- a/scripts/create_testenv.sh
+++ b/scripts/create_testenv.sh
@@ -34,7 +34,7 @@ then
     source activate ${ENVNAME}
 fi
 conda install --yes numpy scipy mkl-service
-conda install --yes -c conda-forge python-graphviz
+conda install --yes -c conda-forge python-graphviz networkx
 
 pip install --upgrade pip