Add auto_deterministics list to Model

Ensures that when missing variables are present in the model, the automatic deterministic (x_observed + x_missing) only appears in predictive sampling and not normal sampling.
Fixes `x` missing from prior_predictive when missing values were present (only `x_missing` was present)

Author: ricardoV94

Diff for: pymc3/model.py

@@ -621,6 +621,7 @@ def __init__(self, name="", model=None, aesara_config=None, coords=None, check_b
             self.rvs_to_values = treedict(parent=self.parent.rvs_to_values)
             self.free_RVs = treelist(parent=self.parent.free_RVs)
             self.observed_RVs = treelist(parent=self.parent.observed_RVs)
+            self.auto_deterministics = treelist(parent=self.parent.auto_deterministics)
             self.deterministics = treelist(parent=self.parent.deterministics)
             self.potentials = treelist(parent=self.parent.potentials)
         else:
@@ -629,6 +630,7 @@ def __init__(self, name="", model=None, aesara_config=None, coords=None, check_b
             self.rvs_to_values = treedict()
             self.free_RVs = treelist()
             self.observed_RVs = treelist()
+            self.auto_deterministics = treelist()
             self.deterministics = treelist()
             self.potentials = treelist()
 
@@ -1076,7 +1078,7 @@ def make_obs_var(
             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)
+            rv_var = Deterministic(name, rv_var, self, dims, auto=True)
 
         elif sps.issparse(data):
             data = sparse.basic.as_sparse(data, name=name)
@@ -1594,21 +1596,28 @@ def __call__(self, *args, **kwargs):
 compilef = fastfn
 
 
-def Deterministic(name, var, model=None, dims=None):
+def Deterministic(name, var, model=None, dims=None, auto=False):
     """Create a named deterministic variable
 
     Parameters
     ----------
     name: str
     var: Aesara variables
+    auto: bool
+        Add automatically created deterministics (e.g., when imputing missing values)
+        to a separate model.auto_deterministics list for filtering during sampling.
+
 
     Returns
     -------
     var: var, with name attribute
     """
     model = modelcontext(model)
     var = var.copy(model.name_for(name))
-    model.deterministics.append(var)
+    if auto:
+        model.auto_deterministics.append(var)
+    else:
+        model.deterministics.append(var)
     model.add_random_variable(var, dims)
 
     return var

Diff for: pymc3/sampling.py

@@ -1679,7 +1679,7 @@ def sample_posterior_predictive(
     if var_names is not None:
         vars_ = [model[x] for x in var_names]
     else:
-        vars_ = model.observed_RVs
+        vars_ = model.observed_RVs + model.auto_deterministics
 
     if random_seed is not None:
         # np.random.seed(random_seed)
@@ -1955,7 +1955,7 @@ def sample_prior_predictive(
         )
 
     if var_names is None:
-        prior_pred_vars = model.observed_RVs
+        prior_pred_vars = model.observed_RVs + model.auto_deterministics
         prior_vars = (
             get_default_varnames(model.unobserved_RVs, include_transformed=True) + model.potentials
         )

Diff for: pymc3/tests/test_idata_conversion.py

@@ -309,7 +309,7 @@ def test_missing_data_model(self):
         assert "y_missing" in model.named_vars
 
         test_dict = {
-            "posterior": ["x", "y", "y_missing"],
+            "posterior": ["x", "y_missing"],
             "observed_data": ["y_observed"],
             "log_likelihood": ["y_observed"],
         }

Diff for: pymc3/tests/test_missing.py

@@ -23,7 +23,7 @@
 from pymc3.distributions.transforms import Interval
 from pymc3.exceptions import ImputationWarning
 from pymc3.model import Model
-from pymc3.sampling import sample, sample_prior_predictive
+from pymc3.sampling import sample, sample_posterior_predictive, sample_prior_predictive
 
 
 @pytest.mark.parametrize(
@@ -125,6 +125,16 @@ def test_interval_missing_observations():
 
         assert np.all(0 < trace["theta1_missing"].mean(0))
         assert np.all(0 < trace["theta2_missing"].mean(0))
+        assert "theta1" not in trace.varnames
+        assert "theta2" not in trace.varnames
+
+        # Make sure that the observed values are newly generated samples and that
+        # the observed and deterministic matche
+        pp_trace = sample_posterior_predictive(trace)
+        assert np.all(np.var(pp_trace["theta1"], 0) > 0.0)
+        assert np.all(np.var(pp_trace["theta2"], 0) > 0.0)
+        assert np.mean(pp_trace["theta1"][:, ~obs1.mask] - pp_trace["theta1_observed"]) == 0.0
+        assert np.mean(pp_trace["theta2"][:, ~obs2.mask] - pp_trace["theta2_observed"]) == 0.0
 
 
 def test_double_counting():
@@ -139,3 +149,17 @@ def test_double_counting():
 
     logp_val = m2.logp({"x_missing_log__": np.array([0])})
     assert logp_val == -4.0
+
+
+def test_missing_logp():
+    with Model() as m:
+        theta1 = Normal("theta1", 0, 5, observed=[0, 1, 2, 3, 4])
+        theta2 = Normal("theta2", mu=theta1, observed=[0, 1, 2, 3, 4])
+    m_logp = m.logp()
+
+    with Model() as m_missing:
+        theta1 = Normal("theta1", 0, 5, observed=np.array([0, 1, np.nan, 3, np.nan]))
+        theta2 = Normal("theta2", mu=theta1, observed=np.array([np.nan, np.nan, 2, np.nan, 4]))
+    m_missing_logp = m_missing.logp({"theta1_missing": [2, 4], "theta2_missing": [0, 1, 3]})
+
+    assert m_logp == m_missing_logp

Diff for: pymc3/tests/test_starting.py

@@ -105,7 +105,7 @@ def test_find_MAP_issue_4488():
         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 not set.difference({"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])