PYMC3 ValueError: Mass matrix contains zeros on the diagonal. Some derivatives might always be zero

Question

I am getting this error while using PYMC3 for calculating posteriors:

with pm.Model() as model:
    p = pm.Gamma('p', alpha=1, beta=3, shape=regions.shape[0])
    q = pm.Gamma('q', alpha=1, beta=3, shape=regions.shape[0])
    m = pm.Lognormal('m', mu=np.log(total_M), sd=.25, shape=regions.shape[0])
    t = pm.Uniform('t', lower=0, upper=100, observed=sales.t)
    cid = pm.Categorical('cid', p=np.repeat(1./sales.shape[0], sales.shape[0]), observed=sales.region )
    sigma = pm.Gamma('sigma', alpha=1, beta=3)
    mu = m[cid]*(((p[cid]+q[cid])**2)/p[cid])*((np.exp(-(p[cid]+q[cid])*t))/((1+(q[cid]/p[cid])*np.exp(-(p[cid]+q[cid])*t))**2))
    Y_obs = pm.Normal('Ft', mu= mu, sd=sigma, observed= sales.sales)
    trace = pm.sample(100000,init = 'adapt_diag', progressbar = True, tune  = 1000)

I have tried changing mu = mu to mu = np.log(mu), it resolved the error but gave me bad results as compared to my other mates.

Can you provide a reproducible example (i.e., examples for `regions` and `sales`)? Does it work if you just call `pm.sample()` at the end? The automatic intialization (`jitter+adapt_diag`) may do a better job, and 100,000 is a lot of samples to take using `NUTS`, unless you have a great reason! — colcarroll, Mar 26 '18 at 02:33
@colcarroll: I have already tried with pm.sample(). It gives same error. Regions are country names and sales are the sales values in numbers. This is Bass model of diffusion for predicting sales. — Ravender Singh Dahiya, Mar 28 '18 at 11:27

score 0 · Answer 1 · answered Mar 30 '18 at 06:22

I have tried with changing the init methods with various options and running them in parallel. It has resolved the issues but it has been running since last 48 hours. Can someone help with any suggestion and feedback on the same?

Below given are my codes with changes in the last few lines:

 1. 

with pm.Model() as model:
    p = pm.Gamma('p', alpha=1, beta=3, shape=regions.shape[0])
    q = pm.Gamma('q', alpha=1, beta=3, shape=regions.shape[0])
    m = pm.Lognormal('m', mu=np.log(total_M), sd=.375, shape=regions.shape[0])
    t = pm.Uniform('t', lower=0, upper=100, observed=sales.t)
   # p1 = pm.Deterministic('p1', np.repeat(1./sales.shape[0],sales.shape[0]))
   # cid = pm.Categorical('cid', p=p1, observed=sales.region )
    cid = pm.Categorical('cid', p=np.repeat(1./sales.shape[0], sales.shape[0]), observed=sales.region )
    sigma = pm.Gamma('sigma', alpha=1, beta=3)
    mu = m[cid]*(((p[cid]+q[cid])**2)/p[
        cid])*((np.exp(-(p[cid]+q[cid])*t))/((1+(q[cid]/p[cid])*np.exp(-(p[cid]+q[cid])*t))**2))
    Y_obs = pm.Normal('Ft', mu=mu, sd=sigma, observed=sales.sales)
    **trace = pm.sample(init = 'advi+adapt_diag', tune  = 1000)**

 2. 

with pm.Model() as model:
    p = pm.Gamma('p', alpha=1, beta=3, shape=regions.shape[0])
    q = pm.Gamma('q', alpha=1, beta=3, shape=regions.shape[0])
    m = pm.Lognormal('m', mu=np.log(total_M), sd=.375, shape=regions.shape[0])
    t = pm.Uniform('t', lower=0, upper=100, observed=sales.t)
   # p1 = pm.Deterministic('p1', np.repeat(1./sales.shape[0],sales.shape[0]))
   # cid = pm.Categorical('cid', p=p1, observed=sales.region )
    cid = pm.Categorical('cid', p=np.repeat(1./sales.shape[0], sales.shape[0]), observed=sales.region )
    sigma = pm.Gamma('sigma', alpha=1, beta=3)
    mu = m[cid]*(((p[cid]+q[cid])**2)/p[
        cid])*((np.exp(-(p[cid]+q[cid])*t))/((1+(q[cid]/p[cid])*np.exp(-(p[cid]+q[cid])*t))**2))
    Y_obs = pm.Normal('Ft', mu=mu, sd=sigma, observed=sales.sales)
    **trace = pm.sample(200000, init = 'advi', tune  = 1000)**

 3. 

with pm.Model() as model:
    p = pm.Gamma('p', alpha=1, beta=3, shape=regions.shape[0])
    q = pm.Gamma('q', alpha=1, beta=3, shape=regions.shape[0])
    m = pm.Lognormal('m', mu=np.log(total_M), sd=.375, shape=regions.shape[0])
    t = pm.Uniform('t', lower=0, upper=100, observed=sales.t)
   # p1 = pm.Deterministic('p1', np.repeat(1./sales.shape[0],sales.shape[0]))
   # cid = pm.Categorical('cid', p=p1, observed=sales.region )
    cid = pm.Categorical('cid', p=np.repeat(1./sales.shape[0], sales.shape[0]), observed=sales.region )
    sigma = pm.Gamma('sigma', alpha=1, beta=3)
    mu = m[cid]*(((p[cid]+q[cid])**2)/p[
        cid])*((np.exp(-(p[cid]+q[cid])*t))/((1+(q[cid]/p[cid])*np.exp(-(p[cid]+q[cid])*t))**2))
    Y_obs = pm.Normal('Ft', mu=mu, sd=sigma, observed=sales.sales)
    **start = pm.find_MAP() 
    step = pm.NUTS(scaling=start)
    trace = pm.sample(5000, step, start = start, progressbar=True)**

Also dealing with a similar issue, I can only imagine that the model should be better specified, and maybe this forum post off of pymc's website can give you some ideas, https://discourse.pymc.io/t/pymc3-slows-rapidly-with-increasing-numbers-of-parameters/134/2 — Coolio2654, Jan 22 '19 at 00:10

PYMC3 ValueError: Mass matrix contains zeros on the diagonal. Some derivatives might always be zero

1 Answers1