Fix shape issues in model comparison

Author: aseyboldt

pymc3/model.py

@@ -915,7 +915,7 @@ def as_tensor(data, name, model, distribution):
 
     if hasattr(data, 'mask'):
         from .distributions import NoDistribution
-        testval = distribution.default()
+        testval = np.broadcast_to(distribution.default(), data.shape)[data.mask]
         fakedist = NoDistribution.dist(shape=data.mask.sum(), dtype=dtype,
                                        testval=testval, parent_dist=distribution)
         missing_values = FreeRV(name=name + '_missing', distribution=fakedist,

pymc3/stats.py

@@ -8,6 +8,7 @@
 from collections import namedtuple
 from .model import modelcontext
 from .util import get_default_varnames
+from pymc3.theanof import floatX
 
 from scipy.misc import logsumexp
 from scipy.stats.distributions import pareto
@@ -123,16 +124,35 @@ def dic(trace, model=None):
     return 2 * mean_deviance - deviance_at_mean
 
 
-def log_post_trace(trace, model):
+def _log_post_trace(trace, model):
     """Calculate the elementwise log-posterior for the sampled trace.
 
     Parameters
     ----------
     trace : result of MCMC run
     model : PyMC Model
         Optional model. Default None, taken from context.
+
+    Returns
+    -------
+    logp : array of shape (n_samples, n_observations)
+        The contribution of the observations to the logp of the whole model.
     """
-    return np.vstack([obs.logp_elemwise(pt) for obs in model.observed_RVs] for pt in trace)
+    def logp_vals_point(pt):
+        if len(model.observed_RVs) == 0:
+            return floatX(np.array([], dtype='d'))
+
+        logp_vals = []
+        for var in model.observed_RVs:
+            logp = var.logp_elemwise(pt)
+            if var.missing_values:
+                logp = logp[~var.observations.mask]
+            logp_vals.append(logp.ravel())
+
+        return np.concatenate(logp_vals)
+
+    logp = (logp_vals_point(pt) for pt in trace)
+    return np.stack(logp)
 
 
 def waic(trace, model=None, pointwise=False):
@@ -160,7 +180,9 @@ def waic(trace, model=None, pointwise=False):
     """
     model = modelcontext(model)
 
-    log_py = log_post_trace(trace, model)
+    log_py = _log_post_trace(trace, model)
+    if log_py.size == 0:
+        raise ValueError('The model does not contain observed values.')
 
     lppd_i = logsumexp(log_py, axis=0, b=1.0 / log_py.shape[0])
 
@@ -210,7 +232,9 @@ def loo(trace, model=None, pointwise=False):
     """
     model = modelcontext(model)
 
-    log_py = log_post_trace(trace, model)
+    log_py = _log_post_trace(trace, model)
+    if log_py.size == 0:
+        raise ValueError('The model does not contain observed values.')
 
     # Importance ratios
     r = np.exp(-log_py)

pymc3/tests/test_model_helpers.py

@@ -97,7 +97,7 @@ def test_as_tensor(self):
         # Create a fake model and fake distribution to be used for the test
         fake_model = pm.Model()
         with fake_model:
-            fake_distribution = pm.Normal('fake_dist', mu=0, sd=1)
+            fake_distribution = pm.Normal.dist(mu=0, sd=1)
             # Create the testval attribute simply for the sake of model testing
             fake_distribution.testval = None
 

pymc3/tests/test_stats.py

@@ -8,11 +8,45 @@
 from ..backends import ndarray
 from ..stats import df_summary, autocorr, hpd, mc_error, quantiles, make_indices
 from ..theanof import floatX_array
+import pymc3.stats as pmstats
 from numpy.random import random, normal
 from numpy.testing import assert_equal, assert_almost_equal, assert_array_almost_equal
 from scipy import stats as st
 
 
+def test_log_post_trace():
+    with pm.Model() as model:
+        pm.Normal('y')
+        trace = pm.sample()
+
+    logp = pmstats._log_post_trace(trace, model)
+    assert logp.shape == (len(trace), 0)
+
+    with pm.Model() as model:
+        pm.Normal('a')
+        pm.Normal('y', observed=np.zeros((2, 3)))
+        trace = pm.sample()
+
+    logp = pmstats._log_post_trace(trace, model)
+    assert logp.shape == (len(trace), 6)
+    npt.assert_allclose(logp, -0.5 * np.log(2 * np.pi), atol=1e-7)
+
+    with pm.Model() as model:
+        pm.Normal('a')
+        pm.Normal('y', observed=np.zeros((2, 3)))
+        data = pd.DataFrame(np.zeros((3, 4)))
+        data.values[1, 1] = np.nan
+        pm.Normal('y2', observed=data)
+        data = data.copy()
+        data.values[:] = np.nan
+        pm.Normal('y3', observed=data)
+        trace = pm.sample()
+
+    logp = pmstats._log_post_trace(trace, model)
+    assert logp.shape == (len(trace), 17)
+    npt.assert_allclose(logp, -0.5 * np.log(2 * np.pi), atol=1e-7)
+
+
 class TestStats(SeededTest):
     @classmethod
     def setup_class(cls):