Hello!
I’m running into a problem when using pm.Gamma.dist(..., beta=k / predicted_spectra) inside a custom model. Sometimes, predicted_spectra contains negative values, which makes the beta parameter invalid and causes this error:
SamplingError: Initial evaluation of model at starting point failed!
To work around it, I added a pm.Potential to enforce positivity:
pm.Potential(“positive_spectra”, pm.math.switch(predicted_spectra > 0, 0, -np.inf))
However, this doesn’t seem to prevent the error — the model still fails to initialize. Here’s a minimal version of my code:
def run_analysis(observed_spectra, time_axis, frequency_axis, betas, predictor_list):
F, T = observed_spectra.shape
N = len(predictor_list)
k = 10
coords = {"T": time_axis, "P": predictor_list, "F": frequency_axis}
with pm.Model(coords=coords) as random_walk_model:
sigma = pm.Uniform("sigma", lower=1.0, upper=10.0)
noise = pm.GaussianRandomWalk("noise", mu=0, sigma=sigma, shape=(N, T), dims=("P", "T"))
predictors = pm.Deterministic("predictors", pm.math.sigmoid(noise))
predicted_spectra = pm.math.dot(betas, predictors)
pm.Potential("positive_spectra", pm.math.switch(predicted_spectra > 0, 0, -np.inf))
log_likelihood_by_time = pm.logp(pm.Gamma.dist(alpha=k, beta=k / predicted_spectra),observed_spectra).sum(axis=0)
pm.Potential("log_likelihood_by_time", log_likelihood_by_time)
idata = pm.sample()
return idata
My question is, why doesn’t the pm.Potential("positive_spectra", ...) prevent the SamplingError from happening at model initialization?
Thank you and sorry for a long question..