diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 7d39b1f4c7..79ed4fff16 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -10,6 +10,16 @@
 - Add log CDF functions to continuous distributions: `Beta`, `Cauchy`, `ExGaussian`, `Exponential`, `Flat`, `Gumbel`, `HalfCauchy`, `HalfFlat`, `HalfNormal`, `Laplace`, `Logistic`, `Lognormal`, `Normal`, `Pareto`, `StudentT`, `Triangular`, `Uniform`, `Wald`, `Weibull`.
 - Behavior of `sample_posterior_predictive` is now to produce posterior predictive samples, in order, from all values of the `trace`. Previously, by default it would produce 1 chain worth of samples, using a random selection from the `trace` (#3212)
 - Show diagnostics for initial energy errors in HMC and NUTS.
+- PR #3273 has added the `distributions.distribution._DrawValuesContext` context
+  manager. This is used to store the values already drawn in nested `random`
+  and `draw_values` calls, enabling `draw_values` to draw samples from the
+  joint probability distribution of RVs and not the marginals. Custom
+  distributions that must call `draw_values` several times in their `random`
+  method, or that invoke many calls to other distribution's `random` methods
+  (e.g. mixtures) must do all of these calls under the same `_DrawValuesContext`
+  context manager instance. If they do not, the conditional relations between
+  the distribution's parameters could be broken, and `random` could return
+  values drawn from an incorrect distribution.
 
 ### Maintenance
 
@@ -17,6 +27,7 @@
 - Fixed Triangular distribution `c` attribute handling in `random` and updated sample codes for consistency (#3225)
 - Refactor SMC and properly compute marginal likelihood (#3124)
 - Removed use of deprecated `ymin` keyword in matplotlib's `Axes.set_ylim` (#3279)
+- Fix for #3210. Now `distribution.draw_values(params)`, will draw the `params` values from their joint probability distribution and not from combinations of their marginals (Refer to PR #3273).
 
 ### Deprecations
 
diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index 5783609765..0dcb6bdb48 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -1,4 +1,4 @@
-import collections
+import six
 import numbers
 
 import numpy as np
@@ -8,7 +8,7 @@
 from ..memoize import memoize
 from ..model import (
     Model, get_named_nodes_and_relations, FreeRV,
-    ObservedRV, MultiObservedRV
+    ObservedRV, MultiObservedRV, Context, InitContextMeta
 )
 from ..vartypes import string_types
 
@@ -214,6 +214,48 @@ def random(self, *args, **kwargs):
                             "Define a custom random method and pass it as kwarg random")
 
 
+class _DrawValuesContext(six.with_metaclass(InitContextMeta, Context)):
+    """ A context manager class used while drawing values with draw_values
+    """
+
+    def __new__(cls, *args, **kwargs):
+        # resolves the parent instance
+        instance = super(_DrawValuesContext, cls).__new__(cls)
+        if cls.get_contexts():
+            potencial_parent = cls.get_contexts()[-1]
+            # We have to make sure that the context is a _DrawValuesContext
+            # and not a Model
+            if isinstance(potencial_parent, cls):
+                instance._parent = potencial_parent
+            else:
+                instance._parent = None
+        else:
+            instance._parent = None
+        return instance
+
+    def __init__(self):
+        if self.parent is not None:
+            # All _DrawValuesContext instances that are in the context of
+            # another _DrawValuesContext will share the reference to the
+            # drawn_vars dictionary. This means that separate branches
+            # in the nested _DrawValuesContext context tree will see the
+            # same drawn values
+            self.drawn_vars = self.parent.drawn_vars
+        else:
+            self.drawn_vars = dict()
+
+    @property
+    def parent(self):
+        return self._parent
+
+
+def is_fast_drawable(var):
+    return isinstance(var, (numbers.Number,
+                            np.ndarray,
+                            tt.TensorConstant,
+                            tt.sharedvar.SharedVariable))
+
+
 def draw_values(params, point=None, size=None):
     """
     Draw (fix) parameter values. Handles a number of cases:
@@ -232,97 +274,134 @@ 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
-        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]
+    # Get fast drawable values (i.e. things in point or numbers, arrays,
+    # constants or shares, or things that were already drawn in related
+    # contexts)
+    if point is None:
+        point = {}
+    with _DrawValuesContext() as context:
+        params = dict(enumerate(params))
+        drawn = context.drawn_vars
+        evaluated = {}
+        symbolic_params = []
+        for i, p in params.items():
+            # If the param is fast drawable, then draw the value immediately
+            if is_fast_drawable(p):
+                v = _draw_value(p, point=point, size=size)
+                evaluated[i] = v
+                continue
+
+            name = getattr(p, 'name', None)
+            if p in drawn:
+                # param was drawn in related contexts
+                v = drawn[p]
+                evaluated[i] = v
+            elif name is not None and name in point:
+                # param.name is in point
+                v = point[name]
+                evaluated[i] = drawn[p] = v
             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])
+                # param still needs to be drawn
+                symbolic_params.append((i, p))
+
+        if not symbolic_params:
+            # We only need to enforce the correct order if there are symbolic
+            # params that could be drawn in variable order
+            return [evaluated[i] for i in params]
+
+        # 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 symbolic_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 = {p.name: (p, v) for p, v in drawn.items()
+                  if getattr(p, 'name', None) is not None}
+        stack = list(leaf_nodes.values())  # A queue would be more appropriate
+        while stack:
+            next_ = stack.pop(0)
+            if next_ in drawn:
+                # 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.
+                continue
+            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
+                    value = _draw_value(next_,
+                                        point=point,
+                                        givens=temp_givens,
+                                        size=size)
+                    givens[next_.name] = (next_, value)
+                    drawn[next_] = value
                 except theano.gof.fg.MissingInputError:
-                    missing_inputs.add(param_idx)
+                    # 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 drawn 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
+        # The remaining params that must be drawn are all hashable
+        to_eval = set()
+        missing_inputs = set([j for j, p in symbolic_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 param in drawn:
+                    evaluated[param_idx] = drawn[param]
+                else:
+                    try:  # might evaluate in a bad order,
+                        value = _draw_value(param,
+                                            point=point,
+                                            givens=givens.values(),
+                                            size=size)
+                        evaluated[param_idx] = drawn[param] = value
+                        givens[param.name] = (param, value)
+                    except theano.gof.fg.MissingInputError:
+                        missing_inputs.add(param_idx)
 
     return [evaluated[j] for j in params] # set the order back
 
@@ -400,8 +479,16 @@ 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
+                    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
                 return dist_tmp.random(point=point, size=size)
             else:
                 return param.distribution.random(point=point, size=size)
@@ -411,10 +498,20 @@ def _draw_value(param, point=None, givens=None, size=None):
             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)])
+            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:
-                return func(*values)
+                output = func(*values)
+            return output
     raise ValueError('Unexpected type in draw_value: %s' % type(param))
 
 
@@ -499,6 +596,20 @@ def generate_samples(generator, *args, **kwargs):
             samples = generator(size=broadcast_shape, *args, **kwargs)
         elif dist_shape == broadcast_shape:
             samples = generator(size=size_tup + dist_shape, *args, **kwargs)
+        elif len(dist_shape) == 0 and size_tup and broadcast_shape[:len(size_tup)] == size_tup:
+            # Input's dist_shape is scalar, but it has size repetitions.
+            # So now the size matches but we have to manually broadcast to
+            # the right dist_shape
+            samples = [generator(*args, **kwargs)]
+            if samples[0].shape == broadcast_shape:
+                samples = samples[0]
+            else:
+                suffix = broadcast_shape[len(size_tup):] + dist_shape
+                samples.extend([generator(*args, **kwargs).
+                                reshape(broadcast_shape)[..., np.newaxis]
+                                for _ in range(np.prod(suffix,
+                                                       dtype=int) - 1)])
+                samples = np.hstack(samples).reshape(size_tup + suffix)
         else:
             samples = None
     # Args have been broadcast correctly, can just ask for the right shape out
@@ -515,9 +626,11 @@ def generate_samples(generator, *args, **kwargs):
     if samples is None:
         raise TypeError('''Attempted to generate values with incompatible shapes:
             size: {size}
+            size_tup: {size_tup}
+            broadcast_shape[:len(size_tup)] == size_tup: {test}
             dist_shape: {dist_shape}
             broadcast_shape: {broadcast_shape}
-        '''.format(size=size, dist_shape=dist_shape, broadcast_shape=broadcast_shape))
+        '''.format(size=size, size_tup=size_tup, dist_shape=dist_shape, broadcast_shape=broadcast_shape, test=broadcast_shape[:len(size_tup)] == size_tup))
 
     # reshape samples here
     if samples.shape[0] == 1 and size == 1:
diff --git a/pymc3/distributions/mixture.py b/pymc3/distributions/mixture.py
index 38f34d6c0a..08a3379be7 100644
--- a/pymc3/distributions/mixture.py
+++ b/pymc3/distributions/mixture.py
@@ -4,7 +4,8 @@
 from pymc3.util import get_variable_name
 from ..math import logsumexp
 from .dist_math import bound, random_choice
-from .distribution import Discrete, Distribution, draw_values, generate_samples
+from .distribution import (Discrete, Distribution, draw_values,
+                           generate_samples, _DrawValuesContext)
 from .continuous import get_tau_sd, Normal
 
 
@@ -147,8 +148,9 @@ def logp(self, value):
                      broadcast_conditions=False)
 
     def random(self, point=None, size=None):
-        w = draw_values([self.w], point=point)[0]
-        comp_tmp = self._comp_samples(point=point, size=None)
+        with _DrawValuesContext() as draw_context:
+            w = draw_values([self.w], point=point)[0]
+            comp_tmp = self._comp_samples(point=point, size=None)
         if np.asarray(self.shape).size == 0:
             distshape = np.asarray(np.broadcast(w, comp_tmp).shape)[..., :-1]
         else:
@@ -163,7 +165,8 @@ def random(self, point=None, size=None):
                                      dist_shape=distshape,
                                      size=size).squeeze()
         if (size is None) or (distshape.size == 0):
-            comp_samples = self._comp_samples(point=point, size=size)
+            with draw_context:
+                comp_samples = self._comp_samples(point=point, size=size)
             if comp_samples.ndim > 1:
                 samples = np.squeeze(comp_samples[np.arange(w_samples.size), ..., w_samples])
             else:
@@ -172,13 +175,14 @@ def random(self, point=None, size=None):
             if w_samples.ndim == 1:
                 w_samples = np.reshape(np.tile(w_samples, size), (size,) + w_samples.shape)
             samples = np.zeros((size,)+tuple(distshape))
-            for i in range(size):
-                w_tmp = w_samples[i, :]
-                comp_tmp = self._comp_samples(point=point, size=None)
-                if comp_tmp.ndim > 1:
-                    samples[i, :] = np.squeeze(comp_tmp[np.arange(w_tmp.size), ..., w_tmp])
-                else:
-                    samples[i, :] = np.squeeze(comp_tmp[w_tmp])
+            with draw_context:
+                for i in range(size):
+                    w_tmp = w_samples[i, :]
+                    comp_tmp = self._comp_samples(point=point, size=None)
+                    if comp_tmp.ndim > 1:
+                        samples[i, :] = np.squeeze(comp_tmp[np.arange(w_tmp.size), ..., w_tmp])
+                    else:
+                        samples[i, :] = np.squeeze(comp_tmp[w_tmp])
 
         return samples
 
diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py
index b807648839..d233055282 100755
--- a/pymc3/distributions/multivariate.py
+++ b/pymc3/distributions/multivariate.py
@@ -16,7 +16,8 @@
 from pymc3.theanof import floatX
 from . import transforms
 from pymc3.util import get_variable_name
-from .distribution import Continuous, Discrete, draw_values, generate_samples
+from .distribution import (Continuous, Discrete, draw_values, generate_samples,
+                           _DrawValuesContext)
 from ..model import Deterministic
 from .continuous import ChiSquared, Normal
 from .special import gammaln, multigammaln
@@ -338,18 +339,19 @@ 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)
-            dist = MvNormal.dist(mu=np.zeros_like(mu), cov=cov)
-        elif self._cov_type == 'tau':
-            tau, = draw_values([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)
-            dist = MvNormal.dist(mu=np.zeros_like(mu), chol=chol)
+        with _DrawValuesContext():
+            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)
+                dist = MvNormal.dist(mu=np.zeros_like(mu), cov=cov)
+            elif self._cov_type == 'tau':
+                tau, = draw_values([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)
+                dist = MvNormal.dist(mu=np.zeros_like(mu), chol=chol)
 
-        samples = dist.random(point, size)
+            samples = dist.random(point, size)
 
         chi2 = np.random.chisquare
         return (np.sqrt(nu) * samples.T / chi2(nu, size)).T + mu
diff --git a/pymc3/model.py b/pymc3/model.py
index 9d705de25a..c9ec85fac6 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -1386,6 +1386,16 @@ def __init__(self, name, data, distribution, total_size=None, model=None):
         self.distribution = distribution
         self.scaling = _get_scaling(total_size, self.logp_elemwiset.shape, self.logp_elemwiset.ndim)
 
+    # Make hashable by id for draw_values
+    def __hash__(self):
+        return id(self)
+
+    def __eq__(self, other):
+        return self.id == other.id
+
+    def __ne__(self, other):
+        return not self == other
+
 
 def _walk_up_rv(rv):
     """Walk up theano graph to get inputs for deterministic RV."""
diff --git a/pymc3/tests/test_random.py b/pymc3/tests/test_random.py
index f61fe3b4a7..1eed3dcd9b 100644
--- a/pymc3/tests/test_random.py
+++ b/pymc3/tests/test_random.py
@@ -1,11 +1,13 @@
 import pymc3 as pm
 import numpy as np
+from numpy import random as nr
 import numpy.testing as npt
 import pytest
 import theano.tensor as tt
 import theano
 
 from pymc3.distributions.distribution import _draw_value, draw_values
+from .helpers import SeededTest
 
 
 def test_draw_value():
@@ -88,3 +90,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,
             ]
         ),
     }