diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py
index 76d05f9d20..94593b9c07 100644
--- a/pymc3/distributions/continuous.py
+++ b/pymc3/distributions/continuous.py
@@ -298,7 +298,7 @@ def __init__(self, sd=None, tau=None, *args, **kwargs):
         assert_negative_support(sd, 'sd', 'HalfNormal')
 
     def random(self, point=None, size=None, repeat=None):
-        sd = draw_values([self.sd], point=point)
+        sd = draw_values([self.sd], point=point)[0]
         return generate_samples(stats.halfnorm.rvs, loc=0., scale=sd,
                                 dist_shape=self.shape,
                                 size=size)
@@ -578,7 +578,7 @@ def __init__(self, lam, *args, **kwargs):
         assert_negative_support(lam, 'lam', 'Exponential')
 
     def random(self, point=None, size=None, repeat=None):
-        lam = draw_values([self.lam], point=point)
+        lam = draw_values([self.lam], point=point)[0]
         return generate_samples(np.random.exponential, scale=1. / lam,
                                 dist_shape=self.shape,
                                 size=size)
@@ -962,7 +962,7 @@ def _random(self, beta, size=None):
         return beta * np.abs(np.tan(np.pi * (u - 0.5)))
 
     def random(self, point=None, size=None, repeat=None):
-        beta = draw_values([self.beta], point=point)
+        beta = draw_values([self.beta], point=point)[0]
         return generate_samples(self._random, beta,
                                 dist_shape=self.shape,
                                 size=size)
diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py
index a15168b07a..0042999dbe 100644
--- a/pymc3/distributions/discrete.py
+++ b/pymc3/distributions/discrete.py
@@ -170,7 +170,7 @@ def __init__(self, p, *args, **kwargs):
         self.mode = tt.cast(tround(p), 'int8')
 
     def random(self, point=None, size=None, repeat=None):
-        p = draw_values([self.p], point=point)
+        p = draw_values([self.p], point=point)[0]
         return generate_samples(stats.bernoulli.rvs, p,
                                 dist_shape=self.shape,
                                 size=size)
@@ -286,7 +286,7 @@ def __init__(self, mu, *args, **kwargs):
         self.mode = tt.floor(mu).astype('int32')
 
     def random(self, point=None, size=None, repeat=None):
-        mu = draw_values([self.mu], point=point)
+        mu = draw_values([self.mu], point=point)[0]
         return generate_samples(stats.poisson.rvs, mu,
                                 dist_shape=self.shape,
                                 size=size)
@@ -398,7 +398,7 @@ def __init__(self, p, *args, **kwargs):
         self.mode = 1
 
     def random(self, point=None, size=None, repeat=None):
-        p = draw_values([self.p], point=point)
+        p = draw_values([self.p], point=point)[0]
         return generate_samples(np.random.geometric, p,
                                 dist_shape=self.shape,
                                 size=size)
@@ -559,7 +559,7 @@ def __init__(self, c, *args, **kwargs):
         self.mean = self.median = self.mode = self.c = c = tt.as_tensor_variable(c)
 
     def random(self, point=None, size=None, repeat=None):
-        c = draw_values([self.c], point=point)
+        c = draw_values([self.c], point=point)[0]
         dtype = np.array(c).dtype
 
         def _random(c, dtype=dtype, size=None):
diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index 178249dd87..d926d27e25 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -1,3 +1,4 @@
+import numbers
 import numpy as np
 import theano.tensor as tt
 from theano import function
@@ -48,7 +49,8 @@ def dist(cls, *args, **kwargs):
         dist.__init__(*args, **kwargs)
         return dist
 
-    def __init__(self, shape, dtype, testval=None, defaults=(), transform=None, broadcastable=None):
+    def __init__(self, shape, dtype, testval=None, defaults=(),
+                 transform=None, broadcastable=None):
         self.shape = np.atleast_1d(shape)
         if False in (np.floor(self.shape) == self.shape):
             raise TypeError("Expected int elements in shape")
@@ -70,8 +72,10 @@ def get_test_val(self, val, defaults):
             return self.getattr_value(val)
 
         if val is None:
-            raise AttributeError(str(self) + " has no finite default value to use, checked: " +
-                                 str(defaults) + " pass testval argument or adjust so value is finite.")
+            raise AttributeError("%s has no finite default value to use, "
+                                 "checked: %s. Pass testval argument or "
+                                 "adjust so value is finite."
+                                 % (self, str(defaults)))
 
     def getattr_value(self, val):
         if isinstance(val, string_types):
@@ -97,7 +101,8 @@ def TensorType(dtype, shape, broadcastable=None):
 
 class NoDistribution(Distribution):
 
-    def __init__(self, shape, dtype, testval=None, defaults=(), transform=None, parent_dist=None, *args, **kwargs):
+    def __init__(self, shape, dtype, testval=None, defaults=(),
+                 transform=None, parent_dist=None, *args, **kwargs):
         super(NoDistribution, self).__init__(shape=shape, dtype=dtype,
                                              testval=testval, defaults=defaults,
                                              *args, **kwargs)
@@ -137,7 +142,8 @@ def __init__(self, shape=(), dtype=None, defaults=('mode',),
 class Continuous(Distribution):
     """Base class for continuous distributions"""
 
-    def __init__(self, shape=(), dtype=None, defaults=('median', 'mean', 'mode'), *args, **kwargs):
+    def __init__(self, shape=(), dtype=None, defaults=('median', 'mean', 'mode'),
+                 *args, **kwargs):
         if dtype is None:
             dtype = theano.config.floatX
         super(Continuous, self).__init__(
@@ -190,18 +196,13 @@ def draw_values(params, point=None):
             if param.name in named_nodes:
                 named_nodes.pop(param.name)
             for name, node in named_nodes.items():
-                if not isinstance(node, (tt.sharedvar.TensorSharedVariable,
+                if not isinstance(node, (tt.sharedvar.SharedVariable,
                                          tt.TensorConstant)):
-                    givens[name] = (node, draw_value(node, point=point))
-    values = [None for _ in params]
-    for i, param in enumerate(params):
-        # "Homogonise" output
-        values[i] = np.atleast_1d(draw_value(
-            param, point=point, givens=givens.values()))
-    if len(values) == 1:
-        return values[0]
-    else:
-        return values
+                    givens[name] = (node, _draw_value(node, point=point))
+    values = []
+    for param in params:
+        values.append(_draw_value(param, point=point, givens=givens.values()))
+    return values
 
 
 @memoize
@@ -229,43 +230,45 @@ def _compile_theano_function(param, vars, givens=None):
                     allow_input_downcast=True)
 
 
-def draw_value(param, point=None, givens=()):
-    if hasattr(param, 'name'):
-        if hasattr(param, 'model'):
-            if point is not None and param.name in point:
-                value = point[param.name]
-            elif hasattr(param, 'random') and param.random is not None:
-                value = param.random(point=point, size=None)
-            else:
-                value = param.tag.test_value
-        else:
-            input_pairs = ([g[0] for g in givens],
-                           [g[1] for g in givens])
+def _draw_value(param, point=None, givens=None):
+    """Draw a random value from a distribution or return a constant.
 
-            value = _compile_theano_function(param,
-                                             input_pairs[0])(*input_pairs[1])
+    Parameters
+    ----------
+    param : number, array like, 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
+        a new value from it and return that, unless a value is specified
+        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
+        are used to evaluate `param` if it is a theano variable.
+    """
+    if isinstance(param, numbers.Number):
+        return param
+    elif isinstance(param, 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):
+        if point and hasattr(param, 'model') and param.name in point:
+            return point[param.name]
+        elif hasattr(param, 'random') and param.random is not None:
+            return param.random(point=point, size=None)
+        else:
+            if givens:
+                variables, values = list(zip(*givens))
+            else:
+                variables = values = []
+            func = _compile_theano_function(param, variables)
+            return func(*values)
     else:
-        value = param
-
-    # Sanitising values may be necessary.
-    if hasattr(value, 'value'):
-        value = value.value
-    elif hasattr(value, 'get_value'):
-        value = value.get_value()
-
-    if hasattr(param, 'dtype'):
-        value = np.atleast_1d(value).astype(param.dtype)
-    if hasattr(param, 'shape'):
-        try:
-            shape = param.shape.tag.test_value
-        except AttributeError:
-            try:
-                shape = param.shape.eval()
-            except AttributeError:
-                shape = param.shape
-        if len(shape) == 0 and len(value) == 1:
-            value = value[0]
-    return value
+        raise ValueError('Unexpected type in draw_value: %s' % type(param))
 
 
 def broadcast_shapes(*args):
@@ -490,8 +493,8 @@ class Bound(object):
         # or within the model context
         NegativeNormal = pm.Bound(pm.Normal, upper=0.0)
         par2 = NegativeNormal('par2', mu=0.0, sd=1.0, testval=1.0)
-        
-        # or you can define it implicitly within the model context 
+
+        # or you can define it implicitly within the model context
         par3 = pm.Bound(pm.Normal, lower=-1.0, upper=1.0)(
                 'par3', mu=0.0, sd=1.0, testval=1.0)
     """
diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py
index be117dedcd..bd6f347f4c 100755
--- a/pymc3/distributions/multivariate.py
+++ b/pymc3/distributions/multivariate.py
@@ -410,7 +410,7 @@ def __init__(self, a, transform=transforms.stick_breaking,
                               np.nan)
 
     def random(self, point=None, size=None):
-        a = draw_values([self.a], point=point)
+        a = draw_values([self.a], point=point)[0]
 
         def _random(a, size=None):
             return stats.dirichlet.rvs(a, None if size == a.shape else size)
diff --git a/pymc3/step_methods/elliptical_slice.py b/pymc3/step_methods/elliptical_slice.py
index 5e6928fa58..968c4de2ae 100644
--- a/pymc3/step_methods/elliptical_slice.py
+++ b/pymc3/step_methods/elliptical_slice.py
@@ -88,7 +88,7 @@ def astep(self, q0, logp):
 
         # Draw from the normal prior by multiplying the Cholesky decomposition
         # of the covariance with draws from a standard normal
-        chol = draw_values([self.prior_chol])
+        chol = draw_values([self.prior_chol])[0]
         nu = np.dot(chol, nr.randn(chol.shape[0]))
         y = logp(q0) - nr.standard_exponential()
 
diff --git a/pymc3/tests/test_random.py b/pymc3/tests/test_random.py
new file mode 100644
index 0000000000..9e04d7c59e
--- /dev/null
+++ b/pymc3/tests/test_random.py
@@ -0,0 +1,92 @@
+import pymc3 as pm
+import numpy as np
+import numpy.testing as npt
+import pytest
+import theano.tensor as tt
+import theano
+
+from pymc3.distributions.distribution import _draw_value, draw_values
+
+
+def test_draw_value():
+    npt.assert_equal(_draw_value(np.array([5, 6])), [5, 6])
+    npt.assert_equal(_draw_value(np.array(5.)), 5)
+
+    npt.assert_equal(_draw_value(tt.constant([5., 6.])), [5, 6])
+    assert _draw_value(tt.constant(5)) == 5
+    npt.assert_equal(_draw_value(2 * tt.constant([5., 6.])), [10, 12])
+
+    val = theano.shared(np.array([5., 6.]))
+    npt.assert_equal(_draw_value(val), [5, 6])
+    npt.assert_equal(_draw_value(2 * val), [10, 12])
+
+    a = tt.scalar('a')
+    a.tag.test_value = 6
+    npt.assert_equal(_draw_value(2 * a, givens=[(a, 1)]), 2)
+
+    assert _draw_value(5) == 5
+    assert _draw_value(5.) == 5
+    assert isinstance(_draw_value(5.), type(5.))
+    assert isinstance(_draw_value(5), type(5))
+
+    with pm.Model():
+        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_value(a)
+    val2 = _draw_value(a)
+    assert np.all(val1 != val2)
+
+    with pytest.raises(ValueError) as err:
+        _draw_value([])
+    err.match('Unexpected type')
+
+
+class TestDrawValues(object):
+    def test_empty(self):
+        assert draw_values([]) == []
+
+    def test_vals(self):
+        npt.assert_equal(draw_values([np.array([5, 6])])[0], [5, 6])
+        npt.assert_equal(draw_values([np.array(5.)])[0], 5)
+
+        npt.assert_equal(draw_values([tt.constant([5., 6.])])[0], [5, 6])
+        assert draw_values([tt.constant(5)])[0] == 5
+        npt.assert_equal(draw_values([2 * tt.constant([5., 6.])])[0], [10, 12])
+
+        val = theano.shared(np.array([5., 6.]))
+        npt.assert_equal(draw_values([val])[0], [5, 6])
+        npt.assert_equal(draw_values([2 * val])[0], [10, 12])
+
+    def test_simple_model(self):
+        with pm.Model():
+            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])
+        val2 = draw_values([a])
+        assert np.all(val1[0] != val2[0])
+
+        point = {'a': np.array([3., 4.])}
+        npt.assert_equal(draw_values([a], point=point), [point['a']])
+
+    def test_dep_vars(self):
+        with pm.Model():
+            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)
+
+        point = {'a': np.array([1., 2.])}
+        npt.assert_equal(draw_values([a], point=point), [point['a']])
+
+        with pytest.raises(theano.gof.MissingInputError):
+            draw_values([a])
+
+        # We need the untransformed vars
+        with pytest.raises(theano.gof.MissingInputError):
+            draw_values([a], point={'sd': np.array([2., 3.])})
+
+        val1 = draw_values([a], point={'sd_log__': np.array([2., 3.])})[0]
+        val2 = draw_values([a], point={'sd_log__': np.array([2., 3.])})[0]
+        assert np.all(val1 != val2)
diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py
index 6f52e1b9f8..ab5dbd9868 100644
--- a/pymc3/tests/test_sampling.py
+++ b/pymc3/tests/test_sampling.py
@@ -13,14 +13,9 @@
 import theano
 from .models import simple_init
 from .helpers import SeededTest
+from scipy import stats
 
-# Test if multiprocessing is available
-import multiprocessing
 import pytest
-try:
-    multiprocessing.Pool(2)
-except:
-    pass
 
 
 @pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
@@ -63,7 +58,8 @@ def test_sample(self):
         with self.model:
             for njobs in test_njobs:
                 for steps in [1, 10, 300]:
-                    pm.sample(steps, tune=0, step=self.step, njobs=njobs, random_seed=self.random_seed)
+                    pm.sample(steps, tune=0, step=self.step, njobs=njobs,
+                              random_seed=self.random_seed)
 
     def test_sample_init(self):
         with self.model:
@@ -72,7 +68,6 @@ def test_sample_init(self):
                           n_init=1000, draws=50,
                           random_seed=self.random_seed)
 
-
     def test_sample_args(self):
         with self.model:
             with pytest.raises(TypeError) as excinfo:
@@ -120,7 +115,7 @@ def test_sample_tune_len(self):
 class TestSoftUpdate(SeededTest):
     def setup_method(self):
         super(TestSoftUpdate, self).setup_method()
-        
+
     def test_soft_update_all_present(self):
         start = {'a': 1, 'b': 2}
         test_point = {'a': 3, 'b': 4}
@@ -187,7 +182,6 @@ def test_constant_named(self):
             assert np.isclose(res, 0.)
 
 
-
 class TestChooseBackend(object):
     def test_choose_backend_none(self):
         with mock.patch('pymc3.sampling.NDArray') as nd:
@@ -209,3 +203,53 @@ def test_choose_backend_shortcut(self):
                                       'name': None}}
         pm.sampling._choose_backend('test_backend', 'chain', shortcuts=shortcuts)
         assert backend.called
+
+
+class TestSamplePPC(object):
+    def test_normal_scalar(self):
+        with pm.Model() as model:
+            a = pm.Normal('a', mu=0, sd=1)
+            trace = pm.sample()
+
+        with model:
+            ppc = pm.sample_ppc(trace, samples=1000, vars=[])
+            assert len(ppc) == 0
+            ppc = pm.sample_ppc(trace, samples=1000, vars=[a])
+            assert 'a' in ppc
+            assert ppc['a'].shape == (1000,)
+        _, pval = stats.kstest(ppc['a'], stats.norm().cdf)
+        assert pval > 0.001
+
+        with model:
+            ppc = pm.sample_ppc(trace, samples=10, size=5, vars=[a])
+            assert ppc['a'].shape == (10, 5)
+
+    def test_normal_vector(self):
+        with pm.Model() as model:
+            a = pm.Normal('a', mu=0, sd=1, shape=2)
+            trace = pm.sample()
+
+        with model:
+            ppc = pm.sample_ppc(trace, samples=10, vars=[])
+            assert len(ppc) == 0
+            ppc = pm.sample_ppc(trace, samples=10, vars=[a])
+            assert 'a' in ppc
+            assert ppc['a'].shape == (10, 2)
+
+            ppc = pm.sample_ppc(trace, samples=10, vars=[a], size=4)
+            assert 'a' in ppc
+            assert ppc['a'].shape == (10, 4, 2)
+
+    def test_sum_normal(self):
+        with pm.Model() as model:
+            a = pm.Normal('a', sd=0.2)
+            b = pm.Normal('b', mu=a)
+            trace = pm.sample()
+
+        with model:
+            ppc = pm.sample_ppc(trace, samples=1000, vars=[b])
+            assert len(ppc) == 1
+            assert ppc['b'].shape == (1000,)
+            scale = np.sqrt(1 + 0.2 ** 2)
+            _, pval = stats.kstest(ppc['b'], stats.norm(scale=scale).cdf)
+            assert pval > 0.001