I have a model that attempts to estimate parameters of a curve that has a linear component followed by an exponential curve. That is, for time-t < t_2 the curve is linear, and for t \geq t_2 the curve is exponential. The conditional equation for the curve is
I’m attempting to estimate all five parameters B, C, D, t_2, H_{ult}. As you can see, for times-t<t_2 parameter D does not exist and i believe this is why i am getting the “Mass matrix contains zeros on the diagonal” error message. Below is the full output:
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [sigma, Hult, t2, D, C, B]
40%|█████████████████████████████████████████████████████████████▋ | 403/1000 [00:12<00:17, 33.28it/s]
Traceback (most recent call last):
File "Heat.py", line 47, in <module>
trace = pm.sample()
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/sampling.py", line 469, in sample
trace = _sample_many(**sample_args)
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/sampling.py", line 515, in _sample_many
step=step, random_seed=random_seed[i], **kwargs)
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/sampling.py", line 559, in _sample
for it, strace in enumerate(sampling):
File "/home/nick/.local/lib/python3.6/site-packages/tqdm/_tqdm.py", line 931, in __iter__
for obj in iterable:
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/sampling.py", line 655, in _iter_sample
point, states = step.step(point)
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py", line 247, in step
apoint, stats = self.astep(array)
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/step_methods/hmc/base_hmc.py", line 115, in astep
self.potential.raise_ok(self._logp_dlogp_func._ordering.vmap)
File "/home/nick/.local/lib/python3.6/site-packages/pymc3/step_methods/hmc/quadpotential.py", line 201, in raise_ok
raise ValueError('\n'.join(errmsg))
ValueError: Mass matrix contains zeros on the diagonal.
The derivative of RV `t2`.ravel()[0] is zero.
The derivative of RV `D`.ravel()[0] is zero.
The code used to obtain this error is:
import numpy as npy
import pymc3 as pm
from theano import tensor as tt
from theano import shared, scan
from theano.ifelse import ifelse
Hd = npy.array([629.52234607, 18969.57367046, 238960.91467992, 328497.85275376])
rows = npy.array([2, 9, 37, 81])
tend = 21
dt = 0.25
tSeq = npy.arange(0, tend+dt, dt)
H0 = tt.as_tensor_variable(npy.asarray(0., dtype='float64'))
def HeatCurve(B, C, D, t2, Hult, tn):
Hn = ifelse(tt.lt(tn,t2),
0.5*Hult*(1-tt.exp(-B*tn**C)),
0.5*Hult*(1-tt.exp(-B*tn**C))+0.5*Hult*(tn-t2)/(tn-t2+D))
return Hn
def CalHeat(tn, Hn_, B, C, D, t2, Hult):
Hn = HeatCurve(B, C, D, t2, Hult, tn)
return Hn
""" Build model """
with pm.Model() as TempModel:
BoundedNormal = pm.Bound(pm.Normal, lower=1e-6, upper=1e6)
B = BoundedNormal('B', mu=0.01172, sd=0.001172, testval=0.01172)
C = BoundedNormal('C', mu=1.6, sd=0.16, testval=1.6)
D = pm.Normal('D', mu=6.2, sd=0.62, testval=6.2)
t2 = pm.Normal('t2', mu=3.5, sd=0.35, testval=3.5)
Hult = BoundedNormal('Hult', mu=338e3, sd=338e2, testval=338e3)
(Hn, _) = scan(
fn = CalHeat,
sequences = shared(tSeq),
outputs_info = H0,
non_sequences = [B, C, D, t2, Hult])
Hn_trim = Hn[rows]
Hm = pm.Deterministic('Hm', Hn_trim)
sigma = pm.HalfNormal('sigma', sd=1)#, shape=len(rows))
H_obs = pm.Normal('H_obs', mu=Hm, sd=sigma, observed=Hd)
trace = pm.sample()
When i do a gradient check with
dHdB, dHdC, dHdD, dHdt2, dHdHult = tt.grad(Hn, [B, C, D, t2, Hult])
the dHdD and dHdt2 are indeed zero for t<t_2 as expected, therefore the error seems to make sense. Does anyone have any suggestions for how i might avoid this error?