Hi,
I’m new to Bayesian calibration as well as PyMC3. I’m trying to apply [Using a “black box” likelihood function — PyMC3 3.10.0 documentation] Using a “black box” likelihood function — PyMC3 3.10.0 documentation) to calibrate parameter m (a list with 2 elements). This is what I did:
Defining likelihood-based deviance:
def poisson_deviance(m, obs_inc):
df = pd.read_csv('dataset.csv')
year_lc = find_year_lc(df, m) #df is dataframe containing all the dependent variable (10+) required for the models
sim_m = find_lc_incidence(year_lc, 'Male')[0]
sim_f = find_lc_incidence(year_lc, 'Female')[0]
obs_inc_m = obs_inc[0]
obs_inc_f = obs_inc[1]
poisson_dev = []
for i, j in zip(obs_inc_m, sim_m):
poisson_dev.append(2 * (i * (np.log(i / j)) - (i - j)))
for m, n in zip(obs_inc_f, sim_f):
poisson_dev.append(2 * (m * (np.log(m / n)) - (m - n)))
return np.sum(poisson_dev)
Theano Op
class LogLike(tt.Op):
itypes = [tt.dvector] # expects a vector of parameter values when called
otypes = [tt.dscalar] # outputs a single scalar value (the log likelihood)
def __init__(self, loglike, data):
# add inputs as class attributes
self.likelihood = loglike
self.data = data
def perform(self, node, inputs, outputs):
# the method that is used when calling the Op
(theta,) = inputs # this will contain my variables
# call the log-likelihood function
logl = self.likelihood(theta, self.data)
outputs[0][0] = np.array(logl) # output the log-likelihood
Lastly, PyMC3 model
# creating Op
logl = LogLike(poisson_deviance, lc_obs_inc) #lc_obs_inc is observed data (array with shape (2,5))
# PyMC3 model
bayes = pm.Model()
with bayes:
# uniform priors on m and c
m = pm.Uniform("m", lower=1e-9, upper=1e-6, shape=2, testval = [7.90e-08, 1.47e-07])
# convert m and c to a tensor vector
theta = tt.as_tensor_variable([m[0], m[1]])
# use a DensityDist (use a lamdba function to "call" the Op)
pm.DensityDist("likelihood", lambda v: logl(v), observed = {'v': theta})
step = pm.Slice()
trace = pm.sample(500, tune=100, step = step, discard_tuned_samples=True)
# plot the traces
_ = pm.traceplot(trace)
And I got this error:
RemoteTraceback:
"""
Traceback (most recent call last):
File "c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\parallel_sampling.py", line 116, in _unpickle_step_method
self._step_method = pickle.loads(self._step_method)
AttributeError: Can't get attribute 'LogLike' on <module '__main__' (built-in)>
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\parallel_sampling.py", line 135, in run
self._unpickle_step_method()
File "c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\parallel_sampling.py", line 118, in _unpickle_step_method
raise ValueError(unpickle_error)
ValueError: The model could not be unpickled. This is required for sampling with more than one core and multiprocessing context spawn or forkserver.
"""
The above exception was the direct cause of the following exception:
ValueError Traceback (most recent call last)
ValueError: The model could not be unpickled. This is required for sampling with more than one core and multiprocessing context spawn or forkserver.
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
<ipython-input-28-8f565fb50b4c> in <module>
12
13 step = pm.Slice()
---> 14 trace = pm.sample(500, tune=100, step = step, discard_tuned_samples=True)
15
16 # plot the traces
c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
556 _print_step_hierarchy(step)
557 try:
--> 558 trace = _mp_sample(**sample_args, **parallel_args)
559 except pickle.PickleError:
560 _log.warning("Could not pickle model, sampling singlethreaded.")
c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
1474 try:
1475 with sampler:
-> 1476 for draw in sampler:
1477 trace = traces[draw.chain - chain]
1478 if trace.supports_sampler_stats and draw.stats is not None:
c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\parallel_sampling.py in __iter__(self)
477
478 while self._active:
--> 479 draw = ProcessAdapter.recv_draw(self._active)
480 proc, is_last, draw, tuning, stats, warns = draw
481 self._total_draws += 1
c:\users\ncc\appdata\local\programs\python\python38\lib\site-packages\pymc3\parallel_sampling.py in recv_draw(processes, timeout)
357 else:
358 error = RuntimeError("Chain %s failed." % proc.chain)
--> 359 raise error from old_error
360 elif msg[0] == "writing_done":
361 proc._readable = True
**RuntimeError:** Chain 0 failed.
I did not specify the gradient since the function is totally a black box and has no analytical form. I have no idea how to reparameterize it. Any idea to solve this problem?
My PyMC3 version is 3.11.1, Theano version 1.0.5, Python 3.8.