I try to estimate a set of subsurface model parameters (say, H1, H2, H3, R1, R2, R3) from earthquake ground motion data. I can predict the data from the model using a numerical method
predictRF(H1, H2, H3, R1, R2, R3), that takes about 0.1 seconds per call. However, the Baysian inference takes overly long. Trying different samplers I run into
NotImplementedErrors, which I would like to understand and ideally get rid of.
The following setup runs for 50 hours:
with open(infile, 'rb') as f: rf_obs = pickle.load(f) with pm.Model() as model: R1 = pm.TruncatedNormal('Ruc', mu=1.75, sigma=0.1, lower=1.42) R2 = pm.TruncatedNormal('Rlvl', mu=1.75, sigma=0.1, lower=1.42) R3 = pm.TruncatedNormal('Rloc', mu=1.75, sigma=0.1, lower=1.42) H1 = pm.TruncatedNormal('Huc', mu=35, sigma=10, lower=0) H2 = pm.TruncatedNormal('Hlvl', mu=5, sigma=1, lower=0) H3 = pm.TruncatedNormal('Haml', mu=5, sigma=5, lower=0) rf_pred = predictRF(predictRF(H1, H2, H3, R1, R2, R3)) L1 = pm.Laplace('L1', mu=rf_pred, b=1, observed=rf_obs) inference = pm.sample(2000, tune=1500, cores=12, return_inferencedata=True)
I need to increase the number of samples to arrive at a stable parameter estimate, but for any initialization method, I get
Initializing NUTS failed. Falling back to elementwise auto-assignment.
inference suggests that the 50 hour runtime is because most steps are discarded. Setting more informative priors and/or a start point close to a parameter estimate I got from simulated annealing doesn’t help. So I wanted to try different samplers:
step=pm.Metropolis() inference = pm.sample(2000, tune=1500, cores=12, step=step return_inferencedata=True)ut things
works, but still takes 20 hours. When I do
step = pm.HamiltonianMC()
I end up with
File "(...)/theano/graph/op.py", line 299, in grad raise NotImplementedError() NotImplementedError
but I wonder what is not implemented here? Interestingly enough, I get the same error with
step = pm.NUTS()
even though I assumed this would dub the default behavior, that worked (see first example).
Could you help me with getting the HamiltonianMC and/or NUTS sampler working? Which parameters do you think might help to get a better acceptance rate and speed up things?
Thanks in adavance!