Add unobserved_value_vars property to Model

Fixes pm.Deterministics still using rvs during pm.Sample()
Fixes find_MAP() failing when pm.Deterministics were present in model

Author: ricardoV94

pymc3/backends/base.py

@@ -61,17 +61,7 @@ def __init__(self, name, model=None, vars=None, test_point=None):
         model = modelcontext(model)
         self.model = model
         if vars is None:
-            vars = []
-            for v in model.unobserved_RVs:
-                var = getattr(v.tag, "value_var", v)
-                transform = getattr(var.tag, "transform", None)
-                if transform:
-                    # We need to create and add an un-transformed version of
-                    # each transformed variable
-                    untrans_var = transform.backward(v, var)
-                    untrans_var.name = v.name
-                    vars.append(untrans_var)
-                vars.append(var)
+            vars = model.unobserved_value_vars
 
         self.vars = vars
         self.varnames = [var.name for var in vars]

pymc3/model.py

@@ -789,6 +789,29 @@ def value_vars(self):
         """
         return [self.rvs_to_values[v] for v in self.free_RVs]
 
+    @property
+    def unobserved_value_vars(self):
+        """List of all random variables (including untransformed projections),
+        as well as deterministics used as inputs and outputs of the the model's
+        log-likelihood graph
+        """
+        vars = []
+        for rv in self.free_RVs:
+            value_var = self.rvs_to_values[rv]
+            transform = getattr(value_var.tag, "transform", None)
+            if transform is not None:
+                # We need to create and add an un-transformed version of
+                # each transformed variable
+                untrans_value_var = transform.backward(rv, value_var)
+                untrans_value_var.name = rv.name
+                vars.append(untrans_value_var)
+            vars.append(value_var)
+
+        # Remove rvs from deterministics graph
+        deterministics, _ = rvs_to_value_vars(self.deterministics, apply_transforms=True)
+
+        return vars + deterministics
+
     @property
     def basic_RVs(self):
         """List of random variables the model is defined in terms of
@@ -803,7 +826,7 @@ def basic_RVs(self):
 
     @property
     def unobserved_RVs(self):
-        """List of all random variable, including deterministic ones.
+        """List of all random variables, including deterministic ones.
 
         These are the actual random variable terms that make up the
         "sample-space" graph (i.e. you can sample these graphs by compiling them
@@ -1049,10 +1072,12 @@ def make_obs_var(
             self.add_random_variable(observed_rv_var, dims)
             self.observed_RVs.append(observed_rv_var)
 
+            # Create deterministic that combines observed and missing
             rv_var = at.zeros(data.shape)
             rv_var = at.set_subtensor(rv_var[mask.nonzero()], missing_rv_var)
             rv_var = at.set_subtensor(rv_var[antimask_idx], observed_rv_var)
             rv_var = Deterministic(name, rv_var, self, dims)
+
         elif sps.issparse(data):
             data = sparse.basic.as_sparse(data, name=name)
             rv_var.tag.observations = data

pymc3/tests/test_missing.py

@@ -80,10 +80,10 @@ def test_missing_dual_observations():
         prior_trace = sample_prior_predictive()
         assert {"beta1", "beta2", "theta", "o1", "o2"} <= set(prior_trace.keys())
         # TODO: Assert something
-        trace = sample(chains=1)
+        trace = sample(chains=1, draws=50)
 
 
-def test_internal_missing_observations():
+def test_interval_missing_observations():
     with Model() as model:
         obs1 = ma.masked_values([1, 2, -1, 4, -1], value=-1)
         obs2 = ma.masked_values([-1, -1, 6, -1, 8], value=-1)
@@ -109,8 +109,8 @@ def test_internal_missing_observations():
         assert prior_trace["theta2"].shape[-1] == obs2.shape[0]
 
         # Make sure that the observed values are newly generated samples
-        assert np.var(prior_trace["theta1_observed"]) > 0.0
-        assert np.var(prior_trace["theta2_observed"]) > 0.0
+        assert np.all(np.var(prior_trace["theta1_observed"], 0) > 0.0)
+        assert np.all(np.var(prior_trace["theta2_observed"], 0) > 0.0)
 
         # Make sure the missing parts of the combined deterministic matches the
         # sampled missing and observed variable values
@@ -121,7 +121,7 @@ def test_internal_missing_observations():
 
         assert {"theta1", "theta2"} <= set(prior_trace.keys())
 
-        trace = sample(chains=1)
+        trace = sample(chains=1, draws=50, compute_convergence_checks=False)
 
         assert np.all(0 < trace["theta1_missing"].mean(0))
         assert np.all(0 < trace["theta2_missing"].mean(0))

pymc3/tests/test_sampling.py

@@ -1098,3 +1098,12 @@ def test_sample_from_xarray_posterior(self, point_list_arg_bug_fixture):
         idat = pm.to_inference_data(trace)
         with pmodel:
             pp = pm.sample_posterior_predictive(idat.posterior, var_names=["d"])
+
+
+def test_sample_deterministic():
+    with pm.Model() as model:
+        x = pm.HalfNormal("x", 1)
+        y = pm.Deterministic("y", x + 100)
+        trace = pm.sample(chains=1, draws=50, compute_convergence_checks=False)
+
+    np.testing.assert_allclose(trace["y"], trace["x"] + 100)

pymc3/tests/test_starting.py

@@ -16,7 +16,17 @@
 
 from pytest import raises
 
-from pymc3 import Beta, Binomial, Model, Normal, Point, Uniform, find_MAP
+from pymc3 import (
+    Beta,
+    Binomial,
+    Deterministic,
+    Gamma,
+    Model,
+    Normal,
+    Point,
+    Uniform,
+    find_MAP,
+)
 from pymc3.tests.checks import close_to
 from pymc3.tests.helpers import select_by_precision
 from pymc3.tests.models import non_normal, simple_arbitrary_det, simple_model
@@ -88,6 +98,18 @@ def test_find_MAP():
     close_to(map_est2["sigma"], 1, tol)
 
 
+def test_find_MAP_issue_4488():
+    # Test for https://github.com/pymc-devs/pymc3/issues/4488
+    with Model() as m:
+        x = Gamma("x", alpha=3, beta=10, observed=np.array([1, np.nan]))
+        y = Deterministic("y", x + 1)
+        map_estimate = find_MAP()
+
+    assert not set.difference({"x", "x_missing", "x_missing_log__", "y"}, set(map_estimate.keys()))
+    assert np.isclose(map_estimate["x_missing"], 0.2)
+    np.testing.assert_array_equal(map_estimate["y"], [2.0, map_estimate["x_missing"][0] + 1])
+
+
 def test_allinmodel():
     model1 = Model()
     model2 = Model()

pymc3/tuning/scaling.py

@@ -60,7 +60,7 @@ def find_hessian(point, vars=None, model=None):
     """
     model = modelcontext(model)
     H = model.fastd2logp(vars)
-    return H(Point(point, model=model))
+    return H(Point(point, filter_model_vars=True, model=model))
 
 
 def find_hessian_diag(point, vars=None, model=None):

pymc3/tuning/starting.py

@@ -157,9 +157,7 @@ def dlogp_func(x):
 
     mx0 = RaveledVars(mx0, x0.point_map_info)
 
-    vars = get_default_varnames(
-        [v.tag.value_var for v in model.unobserved_RVs], include_transformed
-    )
+    vars = get_default_varnames(model.unobserved_value_vars, include_transformed)
     mx = {
         var.name: value
         for var, value in zip(vars, model.fastfn(vars)(DictToArrayBijection.rmap(mx0)))