Make draw_values draw from the joint distribution

pymc3/distributions/continuous.py

@@ -395,7 +395,8 @@ def __init__(self, mu=0, sd=None, tau=None, **kwargs):
         self.sd = tt.as_tensor_variable(sd)
         self.tau = tt.as_tensor_variable(tau)
 
-        self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu)
+        mu = tt.as_tensor_variable(mu)
+        self.mean = self.median = self.mode = self.mu = mu
         self.variance = 1. / self.tau
 
         assert_negative_support(sd, 'sd', 'Normal')

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)

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,101 @@ 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)
+        assert len(sub1.nodes) == 1
+        assert is_directed_acyclic_graph(sub1)
+
+        sub2 = get_sub_dag(dag, [self.a])
+        assert len(sub2.nodes) == 1
+        assert is_directed_acyclic_graph(sub2)
+        assert matching_dependence_dags(sub1, sub2)
+
+        sub3 = get_sub_dag(dag, [self.e])
+        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])
+        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), 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,
+                                             return_index=True)
+        assert is_directed_acyclic_graph(with_non_theano)
+        assert matching_dependence_dags(get_sub_dag(with_non_theano, self.d),
+                                        get_sub_dag(dag, [self.d]))
+        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)
+        wnt2, node2 = add_to_dependence_dag(wnt2,
+                                            params[-1],
+                                            force=True,
+                                            return_added_node=True)
+        assert matching_dependence_dags(wnt2, with_non_theano)
+        assert node2 == index[len(params) - 1]

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,8 @@ 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 +74,8 @@ 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)
@@ -84,7 +87,63 @@ def test_dep_vars(self):
         val2 = draw_values([a], point={'sd': np.array([2., 3.])})[0]
         val3 = draw_values([a], point={'sd_log__': np.array([2., 3.])})[0]
         val4 = draw_values([a], point={'sd_log__': np.array([2., 3.])})[0]
-        
+
         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
+        nr.seed(self.random_seed)
+        N = 1000
+        norm = np.random.randn(N, 3)
+        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 = 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

pymc3/tests/test_step.py

@@ -105,56 +105,46 @@ class TestStepMethods(object):  # yield test doesn't work subclassing object
         1.58740483,  1.67905741,  0.77744868,  0.15050587,  0.15050587,
         0.73979127,  0.15445515,  0.13134717,  0.85068974,  0.85068974,
         0.6974799 ,  0.16170472,  0.86405959,  0.86405959, -0.22032854]),
-       SMC: np.array([ 5.10950205e-02,  1.09811720e+00,  1.78330202e-01,  6.85938766e-01,
-        1.42354476e-01, -1.59630758e+00,  1.57176810e+00, -4.01398917e-01,
-        1.14567871e+00,  1.14954938e+00,  4.94399840e-01,  1.16253017e+00,
-        1.17432244e+00,  7.79195162e-01,  1.29017945e+00,  2.53722905e-01,
-        5.38589898e-01,  3.52121216e-01,  1.35795966e+00,  1.02086933e-01,
-        1.58845251e+00,  6.76852927e-01, -1.04716592e-02, -1.01613324e-01,
-        1.37680965e+00,  7.40036542e-01,  2.89069320e-01,  1.48153741e+00,
-        9.58156958e-01,  5.73623782e-02,  7.68850721e-01,  3.68643390e-01,
-        1.47645964e+00,  2.32596780e-01, -1.85008158e-01,  3.71335958e-01,
-        2.68600102e+00, -4.89504443e-01,  6.54265561e-02,  3.80455349e-01,
-        1.17875338e+00,  2.30233324e-01,  6.90960231e-01,  8.81668685e-01,
-       -2.19754340e-01,  1.27686862e-01,  3.28444250e-01,  1.34820635e-01,
-        5.29725257e-01,  1.43783915e+00, -1.64754264e-01,  7.41446719e-01,
-       -1.17733186e+00,  6.01215658e-02,  1.82638158e-01, -2.23232214e-02,
-       -1.79877583e-02,  8.37949150e-01,  4.41964955e-01, -8.66524743e-01,
-        4.90738093e-01,  2.42056488e-01,  4.67699626e-01,  2.91075351e-01,
-        1.49541153e+00,  8.30730845e-01,  1.03956404e+00, -5.16162910e-01,
-        2.84338859e-01,  1.72305888e+00,  9.52445566e-01,  1.48831718e+00,
-        8.03455325e-01,  1.48840970e+00,  6.98122664e-01,  3.30187139e-01,
-        7.88029712e-01,  9.31510828e-01,  1.01326878e+00,  2.26637755e-01,
-        1.70703646e-01, -8.54429841e-01,  2.97254590e-01, -2.77843274e-01,
-       -2.25544207e-01,  1.98862826e-02,  5.05953885e-01,  4.98203941e-01,
-        1.20897382e+00, -6.32958669e-05, -7.22425896e-01,  1.60930869e+00,
-       -5.02773645e-01,  2.46405678e+00,  9.16039706e-01,  1.14146060e+00,
-       -1.95781984e-01, -2.44653942e-01,  2.67851290e-01,  2.37462012e-01,
-        6.71471950e-01,  1.18319765e+00,  1.29146530e+00, -3.14177753e-01,
-       -1.31041215e-02,  1.05029405e+00,  1.31202399e+00,  7.40532839e-02,
-        9.15510041e-01,  7.71054604e-01,  9.83483263e-01,  9.03032142e-01,
-        9.14191160e-01,  9.32285366e-01,  1.13937607e+00, -4.29155928e-01,
-        3.44609229e-02, -5.46423555e-02,  1.34625982e+00, -1.28287047e-01,
-       -1.55214879e-02,  3.25294234e-01,  1.06120585e+00, -5.09891282e-01,
-        1.25789335e+00,  1.01808348e+00, -9.92590713e-01,  1.72832932e+00,
-        1.12232980e+00,  8.54801892e-01,  1.41534752e+00,  3.50798405e-01,
-        3.69381623e-01,  1.48608411e+00, -1.15506310e-02,  1.57066360e+00,
-        2.00747378e-01,  4.47219763e-01,  5.57720524e-01, -7.74295353e-02,
-        1.79192501e+00,  7.66510475e-01,  1.38852488e+00, -4.06055122e-01,
-        2.73203156e-01,  3.61014687e-01,  1.23574043e+00,  1.64565746e-01,
-       -9.89896480e-02,  9.26130265e-02,  1.06440134e+00, -1.55890408e-01,
-        4.47131846e-01, -7.59186008e-01, -1.50881256e+00, -2.13928005e-01,
-       -4.19160151e-01,  1.75815544e+00,  7.45423008e-01,  6.94781506e-01,
-        1.58596346e+00,  1.75508724e+00,  4.56070434e-01,  2.94128709e-02,
-        1.17703970e+00, -9.90230827e-02,  8.42796845e-01,  1.79154944e+00,
-        5.92779197e-01,  2.73562285e-01,  1.61597907e+00,  1.23514403e+00,
-        4.86261080e-01, -3.10434934e-01,  5.57873722e-01,  6.50365217e-01,
-       -3.41009850e-01,  9.26851109e-01,  8.28936486e-01,  9.16180689e-02,
-        1.30226405e+00,  3.73945789e-01,  6.04560122e-02,  6.00698708e-01,
-        9.68764731e-02,  1.41904148e+00,  6.94182961e-03,  3.17504138e-01,
-        5.90956041e-01, -5.78113887e-01,  5.26615565e-01, -4.19715252e-01,
-        8.92891364e-01,  1.30207363e-01,  4.19899637e-01,  7.10275704e-01,
-        9.27418179e-02,  1.85758044e+00,  4.76988907e-01, -1.36341398e-01]),
+       SMC: np.array([ 0.40152748, -0.1440789 ,  1.87105436,  1.65027354,  0.78140894,
+       -0.33437271,  0.55987446,  1.05976848,  0.52126327,  0.5295624 ,
+       -0.7120724 ,  0.39250673,  0.92590897,  0.776836  ,  0.30528805,
+        1.32178809,  1.30972392,  0.77107019,  1.11885364,  0.59633151,
+        0.63584096, -0.29117982,  0.97372731,  1.06270256,  0.87424729,
+        0.49249202, -0.55942483, -0.17608982,  0.47118016,  1.0026767 ,
+        1.42476886,  1.16505966,  0.71572226,  1.14267914, -0.27628211,
+        0.66712824,  0.58322462,  0.28193361,  0.30175522, -0.11615552,
+       -0.02127047,  0.01085484,  1.21229396,  0.50109798,  0.2046552 ,
+        0.95648093,  0.26673391, -0.703456  ,  1.23223409, -0.87686456,
+        1.45480993,  1.04172093,  1.73512969,  1.00835375,  0.56551883,
+        0.43457948,  1.85267864,  0.51961398,  0.20641743,  0.70484816,
+        1.04491792, -0.70236338,  1.47248532,  0.57438209, -0.15590465,
+        0.51528505,  1.49158593,  0.02418851, -0.04563402,  1.50712686,
+        1.01211014, -0.1058956 ,  1.91153929,  1.09281243,  0.78028316,
+        0.08148316,  0.3989925 ,  0.30230531,  1.59469562, -0.53948736,
+       -0.35653048,  0.44440402,  1.02983002,  0.05184227,  0.78152799,
+        0.99204159,  0.44148902, -0.12657838,  0.97114256,  0.67963455,
+        1.33757129,  0.71977859,  0.09706076, -0.13609892, -0.39969385,
+        0.04687582,  0.053386  ,  0.33382962, -0.36082645,  0.86597207,
+        0.09824643, -0.85212079,  0.54518473, -0.26622955,  0.71836765,
+        0.81359943,  1.39550066,  0.25118273,  1.03965837, -0.65995684,
+       -0.25522586,  2.12497766,  0.69534904,  0.74613619, -0.10312994,
+        1.3244944 , -0.036056  ,  0.90976629,  0.49647046,  0.80779428,
+        0.18921903, -0.18365952,  0.56968353, -0.8232526 , -0.88612154,
+       -0.47326386,  0.18939692,  0.2298177 ,  0.65693251,  1.08908496,
+        1.04748985,  0.53615771, -0.4611776 ,  1.12076823, -0.79971572,
+        1.78908277,  1.32673932,  1.43691077,  0.2564599 ,  0.08480867,
+        0.26340606, -0.86864626,  1.05716355,  0.18611255,  0.44701292,
+       -0.06966819,  0.3325726 ,  0.94594745, -0.0904025 ,  0.14349182,
+        0.83638941,  0.57657934,  0.9549692 , -0.18496471,  0.87838048,
+        0.66938294,  0.54401984,  0.47804147,  0.32545637, -0.82626784,
+        0.93390148,  0.39170683, -0.22244643,  0.36576256,  0.62426937,
+       -0.16594267,  1.55050592,  0.60508809, -1.28925325,  1.1470063 ,
+        0.71030941,  1.20896922,  1.23267962,  0.67278808,  0.5846423 ,
+       -0.09343583, -0.28323718,  0.87891542,  0.54779014,  0.17131075,
+        1.02287448,  0.61819842,  1.28724788,  0.641085  ,  1.48324063,
+       -1.68770188,  0.03750369,  0.47352403,  0.22929128,  0.637757  ,
+        0.61735636,  0.17260147,  1.10929764, -0.33766643,  0.27064342,
+       -0.54594464, -1.23229206, -0.18328842, -0.78636148,  1.38189874]),
     }
 
     def setup_class(self):
@@ -202,7 +192,6 @@ def check_trace(self, step_method):
                 trace = sample(0, tune=n_steps,
                                discard_tuned_samples=False,
                                step=step_method(), random_seed=1, chains=1)
-
         assert_array_almost_equal(
             trace['x'],
             self.master_samples[step_method],
@@ -428,7 +417,7 @@ def test_bad_init(self):
         with Model():
             HalfNormal('a', sd=1, testval=-1, transform=None)
             with pytest.raises(ValueError) as error:
-                sample(init=None)
+                sample(init=None, cores=1)
             error.match('Bad initial')
 
     def test_linalg(self, caplog):