I am trying to change the very simplest getting started - example of pymc3 (https://docs.pymc.io/notebooks/getting_started.html), the motivating example of linear regression into fitting a stretched exponential.
The simplest version of the model I tried is y = exp(-x**beta)
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
# Initialize random number generator
np.random.seed(1234)
# True parameter values
sigma = .1
beta = 1
# Size of dataset
size = 1000
# Predictor variable
X1 = np.random.randn(size)
# Simulate outcome variable
Y = np.exp(-X1**beta) + np.random.randn(size)*sigma
# specify the model
import pymc3 as pm
import theano.tensor as tt
print('Running on PyMC3 v{}'.format(pm.__version__))
basic_model = pm.Model()
with basic_model:
# Priors for unknown model parameters
beta = pm.HalfNormal('beta', sigma=1)
sigma = pm.HalfNormal('sigma', sigma=1)
# Expected value of outcome
mu = pm.math.exp(-X1**beta)
# Likelihood (sampling distribution) of observations
Y_obs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=Y)
with basic_model:
# draw 500 posterior samples
trace = pm.sample(500)
which yields the output
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta]
Sampling 4 chains: 0%| | 0/4000 [00:00<?, ?draws/s]/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
Bad initial energy, check any log probabilities that are inf or -inf, nan or very small:
Y_obs NaN
---------------------------------------------------------------------------
RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
File "/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py", line 160, in _start_loop
point, stats = self._compute_point()
File "/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py", line 191, in _compute_point
point, stats = self._step_method.step(self._point)
File "/opt/conda/lib/python3.7/site-packages/pymc3/step_methods/arraystep.py", line 247, in step
apoint, stats = self.astep(array)
File "/opt/conda/lib/python3.7/site-packages/pymc3/step_methods/hmc/base_hmc.py", line 144, in astep
raise SamplingError("Bad initial energy")
pymc3.exceptions.SamplingError: Bad initial energy
"""
The above exception was the direct cause of the following exception:
SamplingError Traceback (most recent call last)
SamplingError: Bad initial energy
The above exception was the direct cause of the following exception:
ParallelSamplingError Traceback (most recent call last)
<ipython-input-310-782c941fbda8> in <module>
1 with basic_model:
2 # draw 500 posterior samples
----> 3 trace = pm.sample(500)
/opt/conda/lib/python3.7/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, **kwargs)
435 _print_step_hierarchy(step)
436 try:
--> 437 trace = _mp_sample(**sample_args)
438 except pickle.PickleError:
439 _log.warning("Could not pickle model, sampling singlethreaded.")
/opt/conda/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, **kwargs)
967 try:
968 with sampler:
--> 969 for draw in sampler:
970 trace = traces[draw.chain - chain]
971 if (trace.supports_sampler_stats
/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py in __iter__(self)
391
392 while self._active:
--> 393 draw = ProcessAdapter.recv_draw(self._active)
394 proc, is_last, draw, tuning, stats, warns = draw
395 if self._progress is not None:
/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
295 else:
296 error = RuntimeError("Chain %s failed." % proc.chain)
--> 297 raise error from old_error
298 elif msg[0] == "writing_done":
299 proc._readable = True
ParallelSamplingError: Bad initial energy
INFO (theano.gof.compilelock): Waiting for existing lock by process '30255' (I am process '30252')
INFO (theano.gof.compilelock): To manually release the lock, delete /home/jovyan/.theano/compiledir_Linux-4.4--generic-x86_64-with-debian-buster-sid-x86_64-3.7.3-64/lock_dir
/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
out=out, **kwargs)
Instead of the stretched exponential, I have also tried power laws, and sine functions. It seems to me that the problem arises as soon as my model is not injective. Can this be an issue (as apparent, I am a newbie in this field)? Can I restrict sampling to only positive x values? Are there any tricks to this?
So the problem here is that
X1**beta
is only defined when X1 >= 0, or when beta is an integer. When you feed this into your observations, for most places, beta will be a float, and so many of
mu = pm.math.exp(-X1**beta)
will be nan.
I found this out with
>>> basic_model.check_test_point()
beta_log__ -0.77
sigma_log__ -0.77
Y_obs NaN
Name: Log-probability of test_point, dtype: float64
I am not sure what model you are trying to specify! There are ways to require beta to be an integer, and ways to require that X1 be positive, but I would need more details to help you describe the model.