Bambi, Count data and thresholds

Дата публикации: 24-09-2026 09:28:55

tcapretto:
bmb.Prior("HalfNormal", sigma=1)
Thank you all for the help. Apologies for the delay, been on holiday.
For anyone looking in future this is the full example with negative binomial likelihood, editing from the poisson version above. In testing with some dummy data seems to produce sensible figures.
import numpy as np
import pandas as pd
import pymc as pm
import xarray as xr
from scipy.stats import nbinom
import bambi as bmb
class NegativeBinomialWithThreshold(pm.NegativeBinomial):
"""A NegativeBinomial distribution that records ``response > threshold`` as a deterministic."""
##Modified from fantastic forum post here: https://discourse.pymc.io/t/bambi-count-data-and-thresholds/17926
rv_op = pm.NegativeBinomial.rv_op
threshold_name = "TestThreshold"
def __new__(cls, name, mu, threshold, *args, **kwargs):
dims = kwargs.get("dims")
response = super().__new__(cls, name, mu, threshold, *args, **kwargs)
pm.Deterministic(cls.threshold_name, response > threshold, dims=dims)
return response
@classmethod
def dist(cls, mu, threshold, *args, **kwargs):
"""Return the unchanged NegativeBinomial distribution used for prediction and log-likelihood."""
del threshold
return super().dist( mu, *args, **kwargs)
likelihood = bmb.Likelihood(
name="NegativeBinomialWithThreshold",
params=["mu", "threshold","alpha"],
parent="mu",
dist=NegativeBinomialWithThreshold,
)
custom_negativebinomial = bmb.Family(
name="custom_negativebinomial",
likelihood=likelihood,
link={"mu": "log", "threshold": "identity"},
)
custom_negativebinomial.set_default_priors({"threshold": 5, "alpha": bmb.Prior("HalfNormal", sigma=1)})
def probability_above_threshold(model, idata):
"""Return an ``interpret`` transform for ``P(Y > threshold | mu)``.
This works with either a constant threshold or posterior draws from a threshold prior.
"""
if "threshold" in idata.posterior:
threshold = idata.posterior["threshold"]
else:
threshold = model.parameters["threshold"].prior
alpha=idata.posterior["alpha"]
def transform(mu):
# For a count Y, Y > t is equivalent to Y > floor(t).
n = 1.0 / alpha
p = 1.0 / (1.0 + alpha * mu)
return xr.apply_ufunc(nbinom.sf, np.floor(threshold), n,p)
return transform
def run_example(draws=500, tune=500):
"""Fit the custom family and obtain predictions with ``bmb.interpret``.
Returns
-------
dict
The data, fitted model, posterior draws, prediction summaries, and plots.
"""
rng = np.random.default_rng(1234)
x = rng.normal(size=150)
data = pd.DataFrame({"x": x, "y": rng.negative_binomial(len(x),0.5+x/1500)})
# Alternatives for the auxiliary parameter are:
#
# - Omit "threshold" to use the family's default constant, 5.
# - Pass another constant with priors={"threshold": 7}.
# - Pass a prior, as below, to propagate uncertainty about the threshold.
#Due to edits alpha prior is currently required. Needs investigating.
threshold_prior = bmb.Prior("Normal", mu=5, sigma=0.75)
alpha_prior=bmb.Prior("Exponential", lam=1)
model = bmb.Model(
"y ~ x",
data=data,
family=custom_negativebinomial,
priors={"threshold": threshold_prior, "alpha":alpha_prior},
)
idata = model.fit(
draws=draws,
tune=tune,
random_seed=1234,
)
x_grid = np.linspace(data["x"].min(), data["x"].max(), 30)
conditional = {"x": x_grid}
# Expected Poisson count, E[Y | x] = mu.
mean_predictions = bmb.interpret.predictions(
model,
idata,
conditional=conditional,
)
# Posterior predictive counts.
count_predictions = bmb.interpret.predictions(
model,
idata,
conditional=conditional,
target="y",
)

# Transform the posterior expected count into P(Y > threshold | x), using the Poisson survival
# function. This propagates uncertainty in both mu and threshold.
probability_transform = probability_above_threshold(model, idata)
over_threshold_predictions = bmb.interpret.predictions(
model,
idata,
conditional=conditional,
transforms={"y": probability_transform},
)
mean_plot = bmb.interpret.plot_predictions(
model,
idata,
conditional=conditional,
)
over_threshold_plot = bmb.interpret.plot_predictions(
model,
idata,
conditional=conditional,
transforms={"y": probability_transform},
)
return {
"data": data,
"model": model,
"idata": idata,
"mean_predictions": mean_predictions,
"count_predictions": count_predictions,
"over_threshold_predictions": over_threshold_predictions,
"mean_plot": mean_plot,
"over_threshold_plot": over_threshold_plot,
}
if __name__ == "__main__":
example = run_example()
print(example["mean_predictions"].summary)
print(example["over_threshold_predictions"].summary)
example["mean_plot"].show()
example["over_threshold_plot"].show()

Схожие новости

#	Наименование новости	Тональность	Информативность	Дата публикации
1	🚀 Release pymc-extras v0.15.0	0	19.63	11-09-2026
2	Proposal / feedback: topology-aware posterior predictive and simulator summaries	0	11.43	04-10-2026
3	🚀 Release pymc-extras v0.15.1	0	19.63	16-09-2026
4	🚀 Release v6.3.2	0	18.52	08-09-2026
5	Small-N Bayesian Poisson regression, sign flip on leave-one-out — sanity check on approach?	0	9.38	30-08-2026
6	Setting and justifying priors for a discrete "what went wrong" model when I have no labeled data	0	11.32	03-09-2026
7	Genomic epidemiology of the ongoing 2026 Bundibugyo Virus Disease outbreak in the Democratic Republic of the Congo	0	8.06	16-07-2026
8	New contributor looking for guidance: from issue fixes to sustained PyMC contributions	0	10.39	12-09-2026
9	Sampling PyMC models in JupyterLite with a WebAssembly backend for PyTensor	0	6.38	28-09-2026
10	Hiring: Research Consultant (Modeling and Quant. Analyst) - (Remote, US)	0	13.28	24-09-2026

Классификация: . Схожих патентов: 0. Схожих новостей: 10. Тональность: 0. Информативность: 12.8. Источник: discourse.pymc.io.