python-3.x - PYMC3 ValueError : Mass matrix contains zeros on the diagonal. 某些导数可能始终为零

Question

使用 PYMC3 计算后验时出现此错误:

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) 

我尝试将 mu = mu 更改为 mu = np.log(mu)，它解决了错误，但与我的其他伙伴相比，结果很糟糕。

Answer 1

我尝试使用各种选项更改 init 方法并并行运行它们。
它已解决问题，但自过去 48 小时以来一直在运行。
有人可以提供任何建议和反馈吗？

下面给出的是我在最后几行更改的代码:

 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)**