Custom Likelihood Op Fails v3.8 and v3.9 #4057
Comments
|
This worked in v3.7, and I've just downgraded to v3.7 and confirmed it still works (with chains>1) in v3.7: pymc3 version: 3.7
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, mu]
Sampling 4 chains: 100%|██████████████████████████████████████████████████████████████████████████████████| 4000/4000 [00:00<00:00, 4000.38draws/s]I've also tried it on Ubuntu 16.04 (i.e., not a WSL environment), and confirmed that the issues with v3.9.3 and v3.8 (chains > 1) exist there as well. |
|
Thanks for reporting, can confirm the bug. Seems it is relate to io with Arviz - @OriolAbril could you take a look? |
|
This is related to #4002 and arviz-devs/arviz#1279. It is currently an incompatibility between ArviZ and PyMC3. And we have to decide how to best fix it from either side or from both at the same time. TL;DR from the other issue, the root of the problem is using the In the meantime I'd recommend using a Potential instead of a DensityDist, or going all the way to define a custom Distribution. |
|
Excellent - thanks for the recommendation. I can confirm that the following works. Are there any downsides to using Potential instead of DensityDist? import numpy as np
import pymc3 as pm
import theano.tensor as tt
import theano.tests.unittest_tools as utt
class Loglike(tt.Op):
itypes = [tt.dvector]
otypes = [tt.dscalar]
def __init__(self, data):
self.data = data
self.loglike_grad = LoglikeGrad(self.data)
def perform(self, node, inputs, outputs):
theta, = inputs
mu, sigma = theta
logp = -len(self.data)*np.log(np.sqrt(2.0*np.pi)*sigma)
logp += -np.sum((self.data - mu)**2.0)/(2.0*sigma**2.0)
outputs[0][0] = np.array(logp)
def grad(self, inputs, grad_outputs):
theta, = inputs
grads = self.loglike_grad(theta)
return [grad_outputs[0]*grads]
class LoglikeGrad(tt.Op):
itypes = [tt.dvector]
otypes = [tt.dvector]
def __init__(self, data):
self.data = data
def perform(self, node, inputs, outputs):
theta, = inputs
mu, sigma = theta
dmu = np.sum(self.data - mu)/sigma**2.0
dsigma = np.sum((self.data - mu)**2.0)/sigma**3.0 - len(self.data)/sigma
outputs[0][0] = np.array([dmu, dsigma])
def main():
mu_true = 5.0
sigma_true = 1.0
data = np.random.normal(loc=mu_true, scale=sigma_true, size=10000)
loglike = Loglike(data)
utt.verify_grad(loglike, [np.array([3.0, 2.0])])
# verify_grad passes with no errors
with pm.Model() as model:
mu = pm.Normal('mu', mu=4.0, sigma=2.0, testval=4.0)
sigma = pm.HalfNormal('sigma', sigma=5.0, testval=2.0)
theta = tt.as_tensor_variable([mu, sigma])
like = pm.Potential('like', loglike(theta))
with model:
trace = pm.sample()
print(pm.summary(trace).to_string())
if __name__ == "__main__":
print("pymc3 version: ", pm.__version__)
main() |
|
Any updates on this? |
|
ArviZ 0.11.0 added the with ...:
trace = pm.sample(..., compute_convergence_checks=False)
idata = az.from_pymc3(idata, density_dist_obs=False)
# use idata for plotting and stats like ess, summary... *independently* of using az.function or pm.function
# both use the same arviz function under the hood.As this concerns old pymc3 versions, stick with 0.11.0. The .1 and .2 releases were necessary to fix some issues with the latest cmdstanpy release at the time and while fixing the issues with latest cmdstanpy and working seamlessly with latest pymc3 which is where we test, they break arviz compatibility with old pymc3 versions. Things will be completely different in pymc3 4+, the converter to inferencedata will be part of the pymc3 codebase and DensityDist may no longer exist in the same form, therefore I am closing this issue. |
Description of your problem
Implementing a custom likelihood theano.tensor.Op fails in v3.8 with chains > 1 and always in v3.9. See the following example:
This runs fine with pymc v3.8 as long as I run with only 1 chain. If I run on v3.8 with chains > 1, or with any number of chains on v3.9.3, I get the following error after sampling finishes:
Here is the output on v3.8 with chains=1
Here is the full 3.8 traceback with chains = 2:
And here is the full v3.9.3 traceback with chains=anything:
Versions and main components
The text was updated successfully, but these errors were encountered: