That definitely sounds like a substantial issue. I don’t have a CustomDist, and I am not explicitly using anything called Scan. Maybe seeing the model code will give you some ideas.
params = {}
coords = {"health_authority": data["HA"]}
with pm.Model(coords=coords) as model:
scanners = pm.Geometric("scanners", p=0.9, observed=data["count_scanners"])
params["β[fte|pop,scanners]"] = pm.Exponential("β[fte|pop,scanners]", lam=1)
fte = pm.Exponential(
"fte",
lam=1 / (params["β[fte|pop,scanners]"] * norm_lha_pop_size * scanners),
observed=data["current_fte"],
)
params["β[vol|pop]"] = pm.Exponential(
"β[vol|pop]", lam=1, dims="health_authority"
)
params["γ[vol|pop]"] = pm.Exponential(
"γ[vol|pop]", lam=1
)
params["β[vol]"] = pm.Exponential("β[vol]", lam=1)
volume = pm.Poisson(
"volume",
mu=np.abs((params["γ[vol|pop]"] + params["β[vol|pop]"]) * norm_lha_pop_size + params["β[vol]"]),
observed=data["annual_exam_volume"],
)
Only thing that comes to mind off the top of my head, although I don’t know how relevant it is, is that one of the health authorities in coords has only a single observation. Maybe the slope will be poorly-identified.