Computing out-of-sample deviance

This is an interesting use case. I think in principle it should work, for example:

y = theano.shared(np.asarray([6.,8.,8.]))
with pm.Model() as m:
    p=pm.Beta('p', 1.,1.)
    obs = pm.Binomial('y', p=p, n=10, observed=y)
    tr = pm.sample()

In [8]: pm.waic(tr, m)
Out[8]: WAIC_r(WAIC=10.390547057421303, WAIC_se=0.929771615409605, p_WAIC=0.5304594854066094, var_warn=0)

y.set_value(np.ones(4)*8)

In [12]: pm.waic(tr, m)
Out[12]: WAIC_r(WAIC=12.335752053177222, WAIC_se=0.0, p_WAIC=0.43684314201765617, var_warn=0)

So I would suggest checking your model to see if there are:
1, Random variables with explicit shape
2, Deterministic node wrapped in pm.Deterministic. If so remove the wrap:
ymu = pm.Deterministic('mu', X.dot(beta)) --> ymu = X.dot(beta)