statspai.mcmc

mcmc

Bayesian econometrics by MCMC (statspai.mcmc).

NumPy-only samplers for the standard regression models, convergence diagnostics for any chain, and marginal likelihoods. The PyMC-backed estimators for causal designs live in :mod:statspai.bayes.

BMAResult dataclass

Bases: ResultProtocolMixin

Result of :func:bma.

Attributes:

Name	Type	Description
table	DataFrame	One row per coefficient: pip (posterior inclusion probability), post_mean and post_sd (averaged over models, zero where the coefficient is excluded), cond_mean and cond_sd (given inclusion) and sign_certainty (posterior probability, given inclusion, that the coefficient has the sign of cond_mean, from the model-wise normal approximations).
models	DataFrame	The models averaged over, best first: one 0/1 column per term, n_terms, the criterion (bic or log_marglik), r2 (Gaussian) or deviance, and post_prob.
params	Series	table['post_mean'].
std_errors	Series	table['post_sd'].
n_obs	int
n_models	int	Number of models retained.
method, family	str

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame(rng.normal(size=(200, 4)), columns=list("abcd"))
>>> df["y"] = 1 + df["a"] - 0.5 * df["c"] + rng.normal(size=200)
>>> out = sp.bma("y ~ a + b + c + d", df)
>>> isinstance(out, sp.BMAResult)
True
>>> bool(out.table.loc["a", "pip"] > 0.99)
True

pip property

pip: Series

Posterior inclusion probabilities.

predict

predict(data: Optional[DataFrame] = None) -> ndarray

Model-averaged prediction of the mean.

Each retained model's prediction, on the scale of the outcome, weighted by its posterior probability.

plot

plot() -> Any

Posterior inclusion probabilities, largest first.

BayesCausalResult

Bases: ResultProtocolMixin

Posterior of a Bayesian tree causal model (:func:sp.bart_causal).

Attributes:

Name	Type	Description
estimate, se	float	Posterior mean and standard deviation of the estimand (the sample average of the individual effects over the target units).
ci	tuple	Equal-tailed credible interval at level.
draws	Series	Posterior draws of the estimand (all chains, in chain order).
cate	DataFrame	Per-row posterior mean, sd and interval of tau(x_i).
diagnostics	dict	rhat and ess of the estimand across chains, mean leaves per tree, the propensity model, number of chains.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> x = rng.normal(size=n)
>>> z = rng.binomial(1, 1 / (1 + np.exp(-x)))
>>> y = x + 1.0 * z + rng.normal(scale=0.5, size=n)
>>> df = pd.DataFrame({'y': y, 'z': z, 'x': x})
>>> r = sp.bart_causal(df, 'y', 'z', ['x'], method='bcf', draws=200,
...                    burnin=200, chains=1, seed=1)
>>> bool(r.ci[0] < 1.0 < r.ci[1])
True

tau_draws property

tau_draws: ndarray

draws x n posterior draws of the individual effects.

BayesR2Result dataclass

Bases: ResultProtocolMixin

Result of :func:bayes_r2 and :func:loo_r2.

Attributes:

Name	Type	Description
kind	str	'model', 'residual' or 'loo'.
draws	ndarray	Posterior draws of R-squared (kind='loo': Bayesian bootstrap draws around the point estimate).
estimate	float	Median of the draws; for kind='loo' the plug-in value 1 - var(y - yhat_loo) / var(y).
mean, sd, lower, upper	float	Mean, standard deviation and central interval of the draws.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=100)})
>>> df["y"] = df["x"] + rng.normal(size=100)
>>> fit = sp.bayes_regress("y ~ x", df, draws=1000, burnin=300, seed=1)
>>> bool(0.3 < sp.bayes_r2(fit).estimate < 0.7)
True

PPCResult dataclass

Bases: ResultProtocolMixin

Result of :func:ppc.

Attributes:

Name	Type	Description
stat	str	Name of the test statistic.
observed	float	Its value on the data.
replicated	ndarray	Its value on each replicated data set.
p_value	float	Share of replications with a statistic at least as large as the observed one. Values near 0 or 1 mean the model does not reproduce this feature of the data.
y, y_rep	ndarray	The observed outcome and the replications (one row each).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=80)})
>>> df["y"] = 1 + 2 * df["x"] + rng.normal(size=80)
>>> fit = sp.bayes_regress("y ~ x", df, draws=500, burnin=200, seed=1)
>>> out = sp.ppc(fit, stat="max", seed=1)
>>> bool(0.02 < out.p_value < 0.98)
True

plot

plot(kind: str = 'stat', n_overlay: int = 50) -> Any

Draw the check.

kind='stat': histogram of the replicated statistic with the observed value marked. kind='density': the distribution of the observed outcome over those of n_overlay replications.

BayesFactorResult dataclass

Bases: ResultProtocolMixin

Result of :func:bayes_factor.

Attributes:

Name	Type	Description
table	DataFrame	One row per model: log_marglik, prior_prob, post_prob and log_bf (log Bayes factor against the first model, so positive values favour the row's model).
log_bf	float	Log Bayes factor of the first model against the second.
bf	float	The Bayes factor itself (inf when it overflows).
evidence	str	Kass and Raftery's label for 2 log BF and the model it favours.

Examples:

>>> import statspai as sp
>>> out = sp.bayes_factor(-100.0, -103.0)
>>> round(out.log_bf, 1)
3.0
>>> out.evidence
'strong evidence for model 1'

LOOResult dataclass

Bases: ResultProtocolMixin

Expected log predictive density of a model, with its pieces.

Returned by :func:loo, :func:waic and :func:kfold.

Attributes:

Name	Type	Description
kind	{'loo', 'waic', 'kfold'}
estimates	DataFrame	Rows elpd, p (effective number of parameters) and ic (-2 elpd); columns estimate and se. The standard errors measure how much the sum over observations would vary in another sample of the same size, not Monte Carlo error.
pointwise	DataFrame	One row per observation: elpd, p, ic and, for kind='loo', mcse_elpd, pareto_k and n_eff.
elpd, se_elpd, p, se_p, ic, se_ic	float	The entries of estimates.
mcse_elpd	float	Monte Carlo standard error of elpd (kind='loo'); nan when some Pareto shape exceeds the threshold, since the error of those terms is not estimable.
k_threshold	float	Pareto shapes above it mark observations whose leave-one-out density is not reliably estimated.
n_obs, n_draws	int
log_weights	ndarray or None	The smoothed log weights (kind='loo'), kept for :func:loo_predict.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y, mu = rng.normal(size=20), rng.normal(0, 0.2, size=(500, 1))
>>> ll = -0.5 * (y - mu) ** 2
>>> out = sp.loo(ll)
>>> out.kind, out.n_obs
('loo', 20)

pareto_k property

pareto_k: ndarray

Pareto shape by observation (kind='loo').

pareto_k_table

pareto_k_table() -> DataFrame

Counts of observations by reliability of their estimate.

bad_observations

bad_observations() -> ndarray

Positions of the observations whose Pareto shape is too high.

plot

plot() -> Any

Pareto shapes by observation, with the reliability threshold.

PSISResult dataclass

Bases: ResultProtocolMixin

Result of :func:psis.

Attributes:

Name	Type	Description
log_weights	(ndarray, shape(draws, observations))	Smoothed log importance weights, each column normalised to sum to one on the natural scale.
pareto_k	ndarray	Estimated shape of the generalized Pareto tail of the raw ratios.
n_eff	ndarray	Effective sample size of the smoothed weights, r_eff / sum(w ** 2).
k_threshold	float	See :func:k_threshold.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> out = sp.psis(np.random.default_rng(0).normal(size=(500, 2)))
>>> out.pareto_k.shape, out.weights().shape
((2,), (500, 2))

weights

weights() -> ndarray

The smoothed weights on the natural scale.

MCMCDiagnostic dataclass

Bases: ResultProtocolMixin

A convergence diagnostic for each parameter of a chain.

Attributes:

Name	Type	Description
name	str	Which diagnostic this is.
table	DataFrame	One row per parameter. The columns depend on the diagnostic.
settings	dict	The arguments the diagnostic was run with.
passed	bool or None	Whether every parameter passes at the diagnostic's own threshold, when the diagnostic defines one.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> draws = np.random.default_rng(0).normal(size=(2000, 2))
>>> out = sp.geweke_diag(draws)
>>> list(out.table.columns)
['z', 'p_value']
>>> isinstance(out, sp.MCMCDiagnostic)
True

GPResult dataclass

Bases: ResultProtocolMixin

Result of :func:gp_regress.

Attributes:

Name	Type	Description
params	Series	The hyperparameters: one length_scale[<regressor>] per regressor (or a single length_scale), signal_var, noise_var and the constant mean.
fitted	Series	Posterior mean of the regression function at the sample points.
loglik	float	Log marginal likelihood at the hyperparameters. When the constant is estimated this is the restricted likelihood, which is not comparable with a fit that fixes the mean.
kernel, formula, n_obs, level

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = np.sort(rng.uniform(0, 6, 80))
>>> df = pd.DataFrame({"x": x, "y": np.sin(x) + 0.2 * rng.normal(size=80)})
>>> fit = sp.gp_regress("y ~ x", df)
>>> isinstance(fit, sp.GPResult)
True
>>> fit.predict(pd.DataFrame({"x": [1.0, 2.0]})).shape
(2, 4)

log_marginal_likelihood

log_marginal_likelihood() -> float

Log marginal likelihood at the hyperparameters.

predict

predict(newdata: Optional[DataFrame] = None, noise: bool = False, level: Optional[float] = None) -> DataFrame

Posterior mean and interval of the regression function.

Parameters:

Name	Type	Description	Default
newdata	DataFrame	Points at which to predict. Default: the sample.	None
noise	bool	False gives the interval of the regression function, True that of a new observation (adds the noise variance).	False
level	float	Mass of the interval. Default: the level of the fit.	None

Returns:

Type	Description
DataFrame	Columns mean, sd, lower, upper.

expected_improvement

expected_improvement(newdata: DataFrame, minimize: bool = True, best: Optional[float] = None) -> Series

Expected improvement over the best value seen so far.

The criterion of efficient global optimisation (Jones, Schonlau and Welch 1998): large where the predicted value is good, or where it is uncertain. The next point to evaluate is the one that maximises it.

Parameters:

Name	Type	Description	Default
newdata	DataFrame	Candidate points.	required
minimize	bool	Whether small or large outcomes are wanted.	True
best	float	The value to improve on. Default: the smallest (largest) outcome in the sample.	None

Returns:

Type	Description
Series	E[max(best - Y(x), 0)] for each candidate (Y(x) - best when maximising), from the posterior of the regression function.

BayesMixedResult dataclass

Bases: BayesRegressResult

Posterior of a hierarchical model fitted by :func:bayes_mixed.

Everything of :class:BayesRegressResult, plus:

Attributes:

Name	Type	Description
random_effects	DataFrame	Posterior mean and standard deviation of every group's random effects (columns <term> and <term>_sd), indexed by group.
re_cov	DataFrame	Posterior mean of the covariance matrix D of the random effects.
n_groups	int

Notes

predict returns the fixed part x' beta (on the outcome scale of the family at a random effect of zero); add random_effects for a group-specific prediction.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> g = np.repeat(np.arange(30), 6)
>>> df = pd.DataFrame({"id": g, "x": rng.normal(size=180)})
>>> df["y"] = 1 + 0.5 * df["x"] + rng.normal(size=30)[g] + rng.normal(size=180)
>>> fit = sp.bayes_mixed("y ~ x", df, group="id", draws=500, burnin=200, seed=1)
>>> isinstance(fit, sp.BayesMixedResult)
True
>>> list(fit.params.index)
['Intercept', 'x', 'sigma2', 'var(Intercept)']
>>> fit.random_effects.shape
(30, 2)

cite

cite(format: str = 'keys') -> Any

The verified paper.bib key of the sampler.

BayesRegressResult dataclass

Bases: _PredictiveMethods, ResultProtocolMixin

Posterior of a model fitted by :func:bayes_regress.

Attributes:

Name	Type	Description
model	str	The likelihood ('normal', 'logit', ...).
params	Series	Posterior means.
std_errors	Series	Posterior standard deviations.
table	DataFrame	One row per parameter: posterior mean, standard deviation, Monte Carlo standard error of the mean (mcse), effective sample size, the equal-tailed credible interval, the median and the posterior probability of a positive value.
draws	DataFrame	The retained draws, chains stacked, one column per parameter.
chain	ndarray	Chain label of every row of draws.
acceptance_rate	float or None	Metropolis acceptance rate, when the sampler has such a step.
prior	dict	The prior actually used.
n_obs	int
level	float	Mass of the credible interval in table.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=200)})
>>> df["y"] = 1 + 0.5 * df["x"] + rng.normal(size=200)
>>> fit = sp.bayes_regress("y ~ x", df, draws=1000, burnin=200, seed=1)
>>> isinstance(fit, sp.BayesRegressResult)
True
>>> list(fit.params.index)
['Intercept', 'x', 'sigma2']
>>> fit.conf_int().shape
(3, 2)

method property

method: str

'Bayesian <model> regression (<sampler>)', the payload's name.

cite

cite(format: str = 'keys') -> Any

The verified paper.bib keys behind the fitted model.

conf_int

conf_int(level: Optional[float] = None, kind: str = 'equal') -> DataFrame

Credible intervals.

Parameters:

Name	Type	Description	Default
level	float	Posterior mass; defaults to the level the model was fitted with.	None
kind	('equal', 'hpd')	Equal-tailed (quantiles) or highest posterior density.	'equal'

prob

prob(expr: str) -> float

Posterior probability of a statement about the parameters.

expr is evaluated on the draws with the parameter names as variables, e.g. fit.prob("x > 0") or fit.prob("x1 > x2"). Names that are not valid Python identifiers are available through backticks, as in DataFrame.eval.

chain_list

chain_list() -> List[DataFrame]

The retained draws of each chain.

diagnostics

diagnostics() -> Dict[str, Any]

Convergence diagnostics of the retained draws.

Returns a dict with geweke, heidel and raftery (first chain) and, when several chains were run, gelman_rubin. A diagnostic that cannot be computed on this chain (for instance Raftery-Lewis on a short run) carries the reason instead.

predict

predict(data: Optional[DataFrame] = None, what: str = 'mean', level: Optional[float] = None) -> Union[ndarray, DataFrame]

Posterior of the expected outcome at the rows of data.

Parameters:

Name	Type	Description	Default
data	DataFrame	New data; the estimation sample when omitted.	None
what	('mean', 'interval', 'linear', 'draws', 'probabilities')	'mean': posterior mean of E[y | x] (the probability for logit / probit, the rate for count models, the latent mean for tobit, the conditional quantile for quantile regression, the latent index for ordered probit). 'interval': that mean with its credible interval. 'linear': posterior mean of the linear index. 'draws': the full matrix, one row per posterior draw. 'probabilities': for ordered probit, the posterior mean probability of each category.	'mean'
level	float	Mass of the interval for what='interval'.	None

log_marginal_likelihood

log_marginal_likelihood(method: Optional[str] = None, alpha: float = 0.05) -> float

Log marginal likelihood log p(y | model).

Parameters:

Name	Type	Description	Default
method	('exact', 'chib', 'gelfand-dey', 'laplace')	'exact': closed form, conjugate model only. 'chib': Chib (1995) from the Gibbs output; normal and probit models. 'gelfand-dey': Gelfand and Dey (1994) with Geweke's truncated normal weighting density; every model. 'laplace': Laplace approximation at the posterior mode; every model with a smooth posterior (not quantile). Default: exact where it exists, then Chib, then Gelfand-Dey.	'exact'
alpha	float	Tail mass truncated from the weighting density of the Gelfand-Dey estimator.	0.05

Returns:

Type	Description
float

Notes

A marginal likelihood is only meaningful under a proper prior and depends on it: with a prior variance of c on a coefficient the marginal likelihood falls like c^{-1/2} as c grows. marginal_likelihood_details has the numerical standard error of the simulation-based estimates.

marginal_likelihood_details

marginal_likelihood_details(method: Optional[str] = None, alpha: float = 0.05) -> Dict[str, Any]

Log marginal likelihood with its method and numerical error.

plot

plot(params: Optional[Sequence[str]] = None, kind: str = 'trace') -> Any

Trace and density plots of the draws.

Parameters:

Name	Type	Description	Default
params	list of str	Parameters to show; all by default (at most 12).	None
kind	('trace', 'density', 'acf')		'trace'

BARTResult dataclass

Bases: ResultProtocolMixin

Result of :func:bart.

Attributes:

Name	Type	Description
fitted	DataFrame	Posterior mean, standard deviation and interval of the regression function at the sample points (of the success probability for family='binary').
sigma	Series	Draws of the error standard deviation (absent for binary).
variable_importance	Series	Share of all splitting rules that use each regressor, averaged over draws. A screening device, not an effect size.
params	Series	sigma (posterior mean) and the average number of leaves per tree.
n_trees, formula, family, n_obs, level

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.uniform(-2, 2, 200)})
>>> df["y"] = np.sin(2 * df["x"]) + 0.2 * rng.normal(size=200)
>>> fit = sp.bart("y ~ x", df, n_trees=30, draws=200, burnin=100, seed=1)
>>> isinstance(fit, sp.BARTResult)
True
>>> fit.predict(pd.DataFrame({"x": [0.0, 1.0]})).shape
(2, 4)

predict_draws

predict_draws(newdata: DataFrame) -> ndarray

Posterior draws of the regression function at newdata: an array with one row per draw.

predict

predict(newdata: Optional[DataFrame] = None, level: Optional[float] = None) -> DataFrame

Posterior mean and interval of the regression function.

Parameters:

Name	Type	Description	Default
newdata	DataFrame	Default: the sample.	None
level	float		None

Returns:

Type	Description
DataFrame	Columns mean, sd, lower, upper. For family='binary' these describe the success probability.

bayes_arima

bayes_arima(y: Union[str, ndarray, Series], data: Optional[DataFrame] = None, order: Tuple[int, int, int] = (1, 0, 0), constant: Optional[bool] = None, mean_prior: Tuple[float, float] = (0.0, 1000000.0), sigma2_prior: Tuple[float, float] = (0.001, 0.001), horizon: int = 0, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95, tune: Optional[float] = None) -> BayesRegressResult

Bayesian ARIMA(p, d, q) with the exact Gaussian likelihood.

The series is differenced d times and the result is modelled as a stationary and invertible ARMA(p, q) process around a constant. The posterior is sampled by random-walk Metropolis.

Parameters:

Name	Type	Description	Default
y	(array, Series or str)	The series, or a column name of data.	required
data	DataFrame		None
order	(p, d, q)		(1, 0, 0)
constant	bool	Include a constant: the mean of the series when d = 0, the drift when d > 0. Default True when d = 0 and False otherwise.	None
mean_prior	(mean, variance)	Normal prior of the constant.	(0, 1e6)
sigma2_prior	(alpha0, delta0)	sigma2 ~ IG(alpha0 / 2, delta0 / 2).	(0.001, 0.001)
horizon	int	Number of periods to forecast. The posterior predictive mean and interval go to model_info['forecast'].	0
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
chains	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000
tune	float	Scale of the proposal relative to the curvature at the mode. Default 2.38 / sqrt(k).	None

Returns:

Type	Description
BayesRegressResult	model='arima' with parameters const, ar.L1 ..., ma.L1 ... and sigma2. The MA polynomial is 1 + ma.L1 L + ..., the convention of :func:statspai.arima. log_marginal_likelihood() is available for comparing orders with the same d.

Notes

The prior of the AR coefficients is uniform over the stationary region, and that of the MA coefficients is uniform over the invertible region. Both regions are bounded, so the prior is proper without any tuning constant.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.zeros(300)
>>> for t in range(1, 300):
...     y[t] = 0.6 * y[t - 1] + rng.normal()
>>> fit = sp.bayes_arima(y, order=(1, 0, 0), draws=1000, burnin=300, seed=1)
>>> list(fit.params.index)
['const', 'ar.L1', 'sigma2']
>>> bool(0.4 < fit.params["ar.L1"] < 0.8)
True

References

ramirezhassan2026introduction

bayes_bootstrap

bayes_bootstrap(formula: Optional[str] = None, data: Optional[DataFrame] = None, statistic: Optional[Callable[[DataFrame, ndarray], Any]] = None, draws: int = 2000, concentration: float = 1.0, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Bayesian bootstrap of a regression or of any weighted statistic.

Rubin's (1981) analogue of the bootstrap: instead of resampling rows, draw a weight vector from a flat Dirichlet distribution and recompute the statistic with those weights. The draws are a sample from the posterior of the statistic under a nonparametric model that puts all its mass on the observed rows, so no likelihood is assumed. For least squares the posterior standard deviations are close to heteroskedasticity-robust standard errors.

Parameters:

Name	Type	Description	Default
formula	str	'y ~ x1 + x2': weighted least squares coefficients are the statistic.	None
data	DataFrame		None
statistic	callable	Instead of a formula, statistic(data, weights) returning a number, a sequence or a dict / Series of named values. weights sum to one.	None
draws	int	Number of posterior draws.	2000
concentration	float	Dirichlet concentration of the weights. 1 is Rubin's Bayesian bootstrap; larger values pull the weights toward equality.	1
seed	int		None
level	float	Mass of the credible intervals.	0.95

Returns:

Type	Description
BayesRegressResult	With model='bayes_bootstrap'. The draws are independent, so no convergence diagnostic applies, and there is no marginal likelihood.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=200)})
>>> df["y"] = 1 + df["x"] + rng.normal(size=200)
>>> fit = sp.bayes_bootstrap("y ~ x", df, draws=500, seed=1)
>>> list(fit.params.index)
['Intercept', 'x']
>>> med = sp.bayes_bootstrap(
...     data=df, statistic=lambda d, w: np.average(d["y"], weights=w),
...     draws=500, seed=1)

References

rubin1981bayesian

bart_causal

bart_causal(data: DataFrame, y: str, treat: str, covariates: Sequence[str], *, method: str = 'bcf', estimand: str = 'ATE', propensity: Any = 'logit', draws: int = 1000, burnin: int = 1000, thin: int = 1, chains: int = 4, n_trees_control: Optional[int] = None, n_trees_moderate: int = 50, sd_control: float = 2.0, sd_moderate: float = 1.0, numcut: int = 100, min_leaf: int = 5, seed: Optional[int] = None, level: float = 0.95, cut_grid: str = 'auto') -> BayesCausalResult

Bayesian tree models for a binary treatment: BCF or BART.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Continuous outcome.	required
treat	str	0/1 treatment.	required
covariates	sequence of str	Numeric confounders / moderators.	required
method	('bcf', 'bart')	'bcf': Bayesian causal forest (separate prognostic and effect forests, effect forest strongly regularised). 'bart': one BART response surface of the outcome on the covariates, the treatment and (by default) the propensity score; effects are differences of its predictions (bartCause's default response model).	'bcf'
estimand	('ATE', 'ATT', 'ATC')	The sample average of the individual effects over all units, the treated or the controls.	'ATE'
propensity	('logit', 'none')	The propensity score used as a covariate of the prognostic forest ('bcf') or the response surface ('bart'): a logit on the covariates, none, a column name, or one score per row.	'logit'
draws	int	Posterior draws kept per chain, burn-in sweeps, thinning, and independent chains (each with its own seed; R-hat needs >= 2).	1000
burnin	int	Posterior draws kept per chain, burn-in sweeps, thinning, and independent chains (each with its own seed; R-hat needs >= 2).	1000
thin	int	Posterior draws kept per chain, burn-in sweeps, thinning, and independent chains (each with its own seed; R-hat needs >= 2).	1000
chains	int	Posterior draws kept per chain, burn-in sweeps, thinning, and independent chains (each with its own seed; R-hat needs >= 2).	1000
n_trees_control	int	Trees in the prognostic ('bart': the only) forest and in the effect forest. None: 200 for 'bcf' (R bcf), 75 for 'bart' (dbarts::bart2, which bartCause uses).	None, 50
n_trees_moderate	int	Trees in the prognostic ('bart': the only) forest and in the effect forest. None: 200 for 'bcf' (R bcf), 75 for 'bart' (dbarts::bart2, which bartCause uses).	None, 50
sd_control	float	Prior medians of the scale of the prognostic and effect functions, in standard deviations of y ('bcf').	2.0, 1.0
sd_moderate	float	Prior medians of the scale of the prognostic and effect functions, in standard deviations of y ('bcf').	2.0, 1.0
numcut	int	Maximum number of cut points per covariate.	100
min_leaf	int	Smallest number of observations in a leaf.	5
seed	int		None
level	float		0.95
cut_grid	('auto', 'quantile', 'uniform')	Candidate split points of a covariate with more than numcut + 1 distinct values: numcut sample quantiles, or numcut points equally spaced over its range (fewer distinct values: the midpoints between them). 'auto' is 'quantile' for 'bcf' and 'uniform' for 'bart' (dbarts's useQuantiles = FALSE). The grid is part of the prior: with quantile cuts and 200 trees the BART posterior mean of the ATE sat 0.21 posterior SDs from bartCause's, with its own prior 0.02.	'auto'

Returns:

Type	Description
BayesCausalResult

Notes

Evidence. The samplers are checked against known truth on simulated designs (tests/test_bart_causal.py: posterior intervals cover the true ATE and the BCF effect forest recovers effect heterogeneity), and against R bcf 2.0.2 / bartCause 1.0.10 by seed replication on a fixed data set: several seeds on each side, comparing the distributions of the posterior means. That is a stochastic screen (CLAUDE.md §5.1, tier S), not parity; the two sides implement the same model with different samplers.

sp.bcf (the random-forest / gradient-boosting approximation) is unchanged; this function is the posterior-sampling BCF.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(2)
>>> n = 200
>>> x = rng.normal(size=n)
>>> z = rng.binomial(1, 0.5, size=n)
>>> y = x + z * (1 + (x > 0)) + rng.normal(scale=0.5, size=n)
>>> df = pd.DataFrame({'y': y, 'z': z, 'x': x})
>>> r = sp.bart_causal(df, 'y', 'z', ['x'], method='bart', draws=150,
...                    burnin=150, chains=1, seed=3)
>>> r.cate.shape[0]
200

References

[@hahn2020bayesian]

bayes_r2

bayes_r2(fit: Any, kind: str = 'model', level: float = 0.95) -> BayesR2Result

Bayesian R-squared: the share of predictive variance explained.

For each posterior draw, Var(fit) / (Var(fit) + Var(res)) with Var(fit) the variance over observations of the expected outcome. Unlike the classical ratio it cannot exceed one when the prior or a small sample pulls the fit away from least squares, and it comes with a posterior distribution.

Parameters:

Name	Type	Description	Default
fit	fitted model	A result of sp.bayes_regress.	required
kind	('model', 'residual')	'model': Var(res) is the residual variance the model expects: sigma^2 for the Gaussian model, the mean of mu (1 - mu) for binary outcomes, of mu for Poisson and of mu + alpha mu^2 for the negative binomial. For Gaussian and binary models this is what R rstanarm::bayes_R2 reports. 'residual': Var(res) is the variance of y minus the expected outcome of that draw, the first definition in Gelman et al. (2019); it depends on the realised outcomes.	'model'
level	float	Mass of the reported interval.	0.95

Returns:

Type	Description
BayesR2Result

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=100)})
>>> df["y"] = df["x"] + rng.normal(size=100)
>>> fit = sp.bayes_regress("y ~ x", df, draws=1000, burnin=300, seed=1)
>>> out = sp.bayes_r2(fit)
>>> out.draws.shape
(1000,)

References

gelman2019rsquared

loo_r2

loo_r2(fit: Any, loo_result: Optional[LOOResult] = None, n_boot: int = 4000, seed: Optional[int] = None, level: float = 0.95) -> BayesR2Result

Leave-one-out R-squared.

1 - Var(y - yhat_loo) / Var(y), with yhat_loo the prediction of each observation from a posterior that did not see it (:func:loo_predict). It estimates the share of variance the model would explain in new data, so, unlike :func:bayes_r2, it goes down when a regressor adds noise and can be negative.

Parameters:

Name	Type	Description	Default
fit	fitted model	A result of sp.bayes_regress.	required
loo_result	LOOResult	The output of sp.loo(fit); computed when omitted.	None
n_boot	int	Bayesian bootstrap draws (Dirichlet weights on the observations) behind the interval. They describe sampling uncertainty in the two variances, not posterior uncertainty.	4000
seed	int		None
level	float		0.95

Returns:

Type	Description
BayesR2Result	estimate is the plug-in value and does not depend on the seed.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=100)})
>>> df["y"] = df["x"] + rng.normal(size=100)
>>> fit = sp.bayes_regress("y ~ x", df, draws=1000, burnin=300, seed=1)
>>> bool(sp.loo_r2(fit, seed=1).estimate < sp.bayes_r2(fit).estimate + 0.05)
True

References

gelman2019rsquared, vehtari2017practical

mad_sd

mad_sd(x: Any, axis: Optional[int] = 0) -> Any

Median absolute deviation scaled to a standard deviation.

1.4826 * median(|x - median(x)|): equal to the standard deviation for a normal sample and insensitive to the tails. Next to the median it is the robust summary of a set of simulation draws used throughout Gelman, Hill and Vehtari (2020).

Parameters:

Name	Type	Description	Default
x	array - like or DataFrame		required
axis	int or None	Axis to reduce; None flattens.	0

Returns:

Type	Description
(float, ndarray or Series)

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> round(float(sp.mad_sd(np.array([1.0, 2.0, 3.0, 4.0, 100.0]))), 4)
1.4826

References

gelman2020regression

ppc

ppc(fit: Any, stat: Union[str, Callable[[ndarray], float]] = 'mean', draws: Optional[int] = 1000, seed: Optional[int] = None) -> PPCResult

Posterior predictive check of a fitted Bayesian regression.

Simulates replicated data sets from the fitted model, one per posterior draw, and compares a statistic of each with the same statistic of the observed data. A model that cannot produce data looking like the data it was fitted to is wrong in that respect, whatever its coefficients say.

Parameters:

Name	Type	Description	Default
fit	fitted model	A result of sp.bayes_regress.	required
stat	str or callable	'mean', 'sd', 'var', 'min', 'max', 'median', 'skew', 'prop_zero' (share of zeros, the usual check of a count model), or a function of one data vector returning a number.	'mean'
draws	int	Number of replicated data sets (posterior draws used, evenly spaced). None uses every draw.	1000
seed	int		None

Returns:

Type	Description
PPCResult

Notes

Choose statistics the model does not fit directly. The mean of a Gaussian regression with an intercept is reproduced by construction and its check says nothing; the minimum, the maximum, the share of zeros or the skewness can fail. The p-value is a description of the discrepancy, not a calibrated test: it uses the data twice and is conservative.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=200)})
>>> df["y"] = rng.negative_binomial(1, 0.2, size=200)
>>> fit = sp.bayes_regress("y ~ x", df, model="poisson", draws=1000,
...                        burnin=300, seed=1)
>>> bool(sp.ppc(fit, stat="prop_zero", seed=1).p_value < 0.05)
True

References

gabry2019visualization, gelman2020regression

bayes_factor

bayes_factor(*models: Any, names: Optional[Sequence[str]] = None, prior_probs: Optional[Sequence[float]] = None, method: Optional[str] = None) -> BayesFactorResult

Bayes factors and posterior model probabilities.

Parameters:

Name	Type	Description	Default
*models	Any	Two or more fitted models that have a log_marginal_likelihood() method (the results of :func:statspai.bayes_regress), or their log marginal likelihoods as numbers.	()
names	list of str	Labels for the models.	None
prior_probs	list of float	Prior model probabilities; equal by default.	None
method	str	Passed to each model's log_marginal_likelihood.	None

Returns:

Type	Description
BayesFactorResult

Notes

Bayes factors require proper priors and are sensitive to them: making the prior of an extra coefficient more diffuse favours the smaller model without bound (Lindley's paradox). They compare models for the same outcome on the same observations; the function refuses fitted models with different numbers of observations.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=100), "z": rng.normal(size=100)})
>>> df["y"] = 1 + df["x"] + rng.normal(size=100)
>>> m1 = sp.bayes_regress("y ~ x", df, model="conjugate", seed=1)
>>> m2 = sp.bayes_regress("y ~ x + z", df, model="conjugate", seed=1)
>>> bool(sp.bayes_factor(m1, m2).log_bf > 0)
True

References

kass1995bayes

savage_dickey

savage_dickey(result: Any, param: str, value: float = 0.0) -> Dict[str, Any]

Savage-Dickey density ratio for a point restriction.

The Bayes factor of the restriction param = value against the fitted (unrestricted) model equals the posterior density at value divided by the prior density there (Dickey 1971), provided the prior of the other parameters under the restriction is their conditional prior given param = value. With the independent priors of :func:statspai.bayes_regress that holds when the prior covariance is diagonal. Under model='conjugate' the coefficients' prior is scaled by sigma2, so conditioning on param = value also updates sigma2: the ratio is then the Bayes factor against the restricted model whose prior is sigma2 ~ InvGamma((alpha0 + 1) / 2, (delta0 + (value - b0)^2 / B0) / 2), not InvGamma(alpha0 / 2, delta0 / 2).

Parameters:

Name	Type	Description	Default
result	BayesRegressResult		required
param	str	A coefficient of the model.	required
value	float		0

Returns:

Type	Description
dict	bf01 (restriction against unrestricted), log_bf01, posterior_density, prior_density, estimator (how the posterior density was obtained) and evidence.

Notes

The posterior density is exact for the conjugate model (Student-t), a Rao-Blackwell average of the full conditionals for the normal and probit models, and a kernel estimate from the draws otherwise. A kernel estimate is unreliable far in the tail; the function then reports the density of the normal approximation and says so.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=200), "z": rng.normal(size=200)})
>>> df["y"] = 1 + df["x"] + rng.normal(size=200)
>>> fit = sp.bayes_regress("y ~ x + z", df, model="conjugate", seed=1)
>>> bool(sp.savage_dickey(fit, "z")["bf01"] > 1)
True

References

dickey1971weighted, verdinelli1995computing

kfold

kfold(fit: Any, k: int = 10, folds: Optional[Any] = None, seed: Optional[int] = None, refit: Optional[Callable[[DataFrame], Any]] = None, data: Optional[DataFrame] = None) -> LOOResult

K-fold cross-validation of a Bayesian model by refitting.

Holds out each fold in turn, refits on the rest, and scores every held-out observation by its log predictive density averaged over the posterior draws of the fit that did not see it. Slower than :func:loo and free of its importance-sampling approximation: the choice when Pareto shapes are high, and the only valid one when whole groups must be left out together.

Parameters:

Name	Type	Description	Default
fit	fitted model	A result of sp.bayes_regress. It is refitted with the same settings on each training set.	required
k	int		10
folds	array-like of int	Fold of each observation, e.g. from :func:kfold_split with groups=. Overrides k and seed.	None
seed	int	Seed of the random fold assignment.	None
refit	callable	refit(train) -> fitted model for a model that does not carry its own refitting recipe; the returned object must have log_lik(data). Requires data.	None
data	DataFrame	The estimation data, when refit is given.	None

Returns:

Type	Description
LOOResult	With kind='kfold'. p is the within-sample log predictive density of the full fit minus the cross-validated one.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=60)})
>>> df["y"] = 1 + 2 * df["x"] + rng.normal(size=60)
>>> fit = sp.bayes_regress("y ~ x", df, draws=500, burnin=200, seed=1)
>>> out = sp.kfold(fit, k=5, seed=1)
>>> out.kind
'kfold'

References

vehtari2017practical

kfold_split

kfold_split(n: int, k: int = 10, seed: Optional[int] = None, groups: Optional[Any] = None, stratify: Optional[Any] = None) -> ndarray

Assign n observations to k folds of near-equal size.

Parameters:

Name	Type	Description	Default
n	int		required
k	int		10
seed	int		None
groups	array - like	Keep the observations of a group in one fold (clustered data).	None
stratify	array - like	Spread the levels of this variable evenly over folds.	None

Returns:

Type	Description
ndarray of int	Fold of each observation, 0 .. k - 1.

Examples:

>>> import statspai as sp
>>> sorted(set(sp.kfold_split(10, k=5, seed=1).tolist()))
[0, 1, 2, 3, 4]

loo

loo(x: Any, r_eff: Any = None, save_weights: bool = True) -> LOOResult

Leave-one-out cross-validation by Pareto smoothed importance sampling.

Estimates the expected log predictive density for new data, elpd = sum_i log p(y_i | y_{-i}), from one fit: the posterior draws are reweighted to stand in for the posterior without observation i. The Pareto shape of each observation's weights says whether that reweighting can be trusted.

Parameters:

Name	Type	Description	Default
x	fitted model or array, shape (draws, observations)	A result with a log_lik() method (sp.bayes_regress), or the pointwise log-likelihood matrix.	required
r_eff	float or array	Relative efficiency of the draws by observation. For a fitted model it is computed from the model's chains; for a bare matrix it defaults to 1 (independent draws).	None
save_weights	bool	Keep the smoothed log weights on the result, for :func:loo_predict and :func:loo_r2.	True

Returns:

Type	Description
LOOResult

Notes

p is the effective number of parameters, the difference between the within-sample log predictive density and elpd. A value far above the number of parameters in the model, or many high Pareto shapes, points to a misspecified model or a few very influential observations; :func:kfold is the remedy that always works.

The reliability threshold on the Pareto shape depends on the number of draws (see :func:k_threshold). The fixed 0.7 of older software and of Vehtari, Gelman and Gabry (2017) is its limit for long chains.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=80)})
>>> df["y"] = 1 + 2 * df["x"] + rng.normal(size=80)
>>> fit = sp.bayes_regress("y ~ x", df, draws=1000, burnin=300, seed=1)
>>> out = sp.loo(fit)
>>> bool(1.5 < out.p < 5)
True

References

vehtari2017practical, vehtari2024pareto

loo_compare

loo_compare(*models: Any, names: Optional[Sequence[str]] = None) -> DataFrame

Compare models on expected log predictive density.

Parameters:

Name	Type	Description	Default
*models	LOOResult, or a single dict of them	Results of :func:loo, :func:waic or :func:kfold computed on the same observations. Fitted models are accepted and passed through :func:loo first.	()
names	list of str	Labels; model1, model2, ... by default, or the keys of the dict.	None

Returns:

Type	Description
DataFrame	One row per model, best first. elpd_diff is the difference from the best model and se_diff its standard error, computed from the paired pointwise differences (much smaller than the standard errors of the two totals when the models' predictions are correlated, which they nearly always are). The remaining columns repeat each model's elpd, p and ic.

Notes

A difference of less than about 4 is small whatever its standard error; for larger ones compare elpd_diff with se_diff. The normal approximation behind se_diff is itself unreliable with fewer than about 100 observations or when the models are nearly identical.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y, mu = rng.normal(size=30), rng.normal(0, 0.2, size=(500, 1))
>>> a = -0.5 * (y - mu) ** 2
>>> out = sp.loo_compare(sp.loo(a), sp.loo(a - 0.1), names=["a", "b"])
>>> out.index[0]
'a'

References

vehtari2017practical

loo_predict

loo_predict(fit: Any, loo_result: Optional[LOOResult] = None, what: str = 'mean', draws: Optional[Any] = None) -> ndarray

Leave-one-out prediction of each observation from the other ones.

The posterior expectation of E[y_i | x_i] reweighted to the posterior that did not see observation i. Unlike fitted values, these do not benefit from having been fitted to the point they predict, so residuals from them are honest.

Parameters:

Name	Type	Description	Default
fit	fitted model	A result with predict(what='draws') (sp.bayes_regress).	required
loo_result	LOOResult	The output of sp.loo(fit); computed when omitted.	None
what	('mean', 'linear')	Scale of the prediction: the expected outcome or the linear index. 'mean' is what posterior_epred returns, so for model='tobit', 'oprobit' and 'ologit' it is the latent linear index and for model='quantile' the fitted conditional quantile, not the mean of the observed outcome.	'mean'
draws	(array, shape(draws, observations))	Predict this quantity instead, e.g. a matrix of your own posterior functionals.	None

Returns:

Type	Description
ndarray	One prediction per observation.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=60)})
>>> df["y"] = 1 + 2 * df["x"] + rng.normal(size=60)
>>> fit = sp.bayes_regress("y ~ x", df, draws=500, burnin=200, seed=1)
>>> sp.loo_predict(fit).shape
(60,)

References

vehtari2017practical

psis

psis(log_ratios: Any, r_eff: Any = None) -> PSISResult

Pareto smoothed importance sampling.

Stabilises importance weights by replacing the largest ones with the expected order statistics of a generalized Pareto distribution fitted to the upper tail. The fitted shape k is also the diagnostic: it measures how heavy the tail of the ratios is, and so how many draws a reliable estimate needs.

Parameters:

Name	Type	Description	Default
log_ratios	(array, shape(draws) or (draws, observations))	Log importance ratios. For leave-one-out cross-validation these are minus the pointwise log-likelihoods.	required
r_eff	float or array	Relative efficiency of the draws (effective sample size divided by the number of draws) for each column. Default 1, which is right for independent draws and optimistic for a slowly mixing chain. It sets the tail length and scales n_eff.	None

Returns:

Type	Description
PSISResult

Notes

The tail has ceil(min(0.2 S, 3 sqrt(S / r_eff))) draws. Smoothed weights are truncated at the largest raw weight. A column whose tail has fewer than five distinct values is left unsmoothed and gets k = inf.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> out = sp.psis(rng.normal(size=(1000, 3)))
>>> out.log_weights.shape
(1000, 3)
>>> bool(np.allclose(np.exp(out.log_weights).sum(axis=0), 1.0))
True

References

vehtari2024pareto, zhang2009new

waic

waic(x: Any) -> LOOResult

Widely applicable information criterion.

elpd_waic = sum_i [log mean_s p(y_i | theta_s) - var_s log p(y_i | theta_s)]: the within-sample log predictive density minus the posterior variance of each pointwise log-likelihood (Watanabe 2010). Asymptotically equal to leave-one-out cross-validation; in finite samples :func:loo is more robust and carries its own diagnostic, so prefer it.

Parameters:

Name	Type	Description	Default
x	fitted model or array, shape (draws, observations)		required

Returns:

Type	Description
LOOResult	With kind='waic'.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y, mu = rng.normal(size=20), rng.normal(0, 0.2, size=(500, 1))
>>> ll = -0.5 * (y - mu) ** 2
>>> round(sp.waic(ll).ic, 1) == round(-2 * sp.waic(ll).elpd, 1)
True

References

watanabe2010asymptotic, vehtari2017practical

gelman_rubin

gelman_rubin(chains: Any, confidence: float = 0.95, split: bool = False, autoburnin: bool = False) -> MCMCDiagnostic

Gelman and Rubin (1992) potential scale reduction factor.

Compares the variance between several chains started from dispersed values with the variance within them. A factor near 1 says the chains have mixed; the conventional alarm is above 1.1 (1.01 for the split version). Point estimate and upper confidence limit as in R coda::gelman.diag (with the Brooks and Gelman 1998 degrees of freedom correction), plus the multivariate factor.

Parameters:

Name	Type	Description	Default
chains	list of chains, array (chains x draws x parameters), or a fit	At least two chains of equal length. A fitted model that ran several chains is read chain by chain (its chain_list()).	required
confidence	float	Coverage of the upper limit.	0.95
split	bool	Split every chain in two before computing the factor, which also detects a single chain that drifts. Use this when only one chain is available.	False
autoburnin	bool	Discard the first half of every chain first. coda::gelman.diag does this by default; here the draws are assumed to be post burn-in already.	False

Returns:

Type	Description
MCMCDiagnostic	table columns psrf and upper; settings['mpsrf'] is the multivariate factor (None for one parameter).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(7)
>>> out = sp.gelman_rubin([rng.normal(size=(1000, 2)) for _ in range(3)])
>>> bool((out.table["psrf"] < 1.05).all())
True

References

gelman1992inference, brooks1998general

geweke_diag

geweke_diag(chain: ChainLike, frac1: float = 0.1, frac2: float = 0.5) -> MCMCDiagnostic

Geweke (1992) convergence diagnostic.

Compares the mean of the first frac1 of the chain with the mean of the last frac2. If the chain is stationary the two means are equal and the standardised difference is asymptotically standard normal. Standard errors use the spectral density at zero of each window, so autocorrelation is accounted for. Same statistic as R coda::geweke.diag.

Parameters:

Name	Type	Description	Default
chain	array-like, DataFrame or fitted Bayesian result		required
frac1	float	Fraction of the chain in the early window.	0.1
frac2	float	Fraction of the chain in the late window.	0.5

Returns:

Type	Description
MCMCDiagnostic	table has the z-score and its two-sided p-value by parameter; passed is True when every |z| < 1.96.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> x = np.random.default_rng(4).normal(size=3000)
>>> out = sp.geweke_diag(x)
>>> bool(abs(out.table["z"].iloc[0]) < 4)
True

References

geweke1992evaluating

heidel_diag

heidel_diag(chain: ChainLike, eps: float = 0.1, pvalue: float = 0.05) -> MCMCDiagnostic

Heidelberger and Welch (1983) stationarity and half-width tests.

Stationarity: a Cramer-von Mises test that the chain is a draw from a stationary process, applied to the whole chain and then, if it rejects, after discarding the first 10%, 20%, ... up to 50%. The half-width test then asks whether the part of the chain retained estimates the mean precisely enough: it passes when the half-width of the 95% interval for the mean, relative to the mean, is below eps. Same procedure as R coda::heidel.diag.

Parameters:

Name	Type	Description	Default
chain	array-like, DataFrame or fitted Bayesian result		required
eps	float	Target ratio of half-width to mean.	0.1
pvalue	float	Significance level of the stationarity test.	0.05

Returns:

Type	Description
MCMCDiagnostic	table columns: stationary (bool), start (first iteration retained, 1-based), p_value, halfwidth_passed, mean and halfwidth. For a chain that fails stationarity the last four are missing.

Notes

The half-width test is relative to the mean, so it is uninformative for a parameter whose posterior mean is near zero.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> x = 5 + np.random.default_rng(6).normal(size=4000)
>>> out = sp.heidel_diag(x)
>>> bool(out.table["stationary"].iloc[0])
True

References

heidelberger1983simulation

hpd_interval

hpd_interval(chain: ChainLike, prob: float = 0.95) -> DataFrame

Highest posterior density interval of each parameter.

The shortest interval that holds prob of the draws (Chen and Shao 1999). For a unimodal posterior it is the HPD region; for a symmetric one it coincides with the equal-tailed interval. Same rule as R coda::HPDinterval.

Parameters:

Name	Type	Description	Default
chain	array-like, DataFrame or fitted Bayesian result		required
prob	float	Posterior mass inside the interval.	0.95

Returns:

Type	Description
DataFrame	Columns lower and upper, one row per parameter.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> x = np.random.default_rng(2).exponential(size=20000)
>>> iv = sp.hpd_interval(x, prob=0.9)
>>> bool(iv["lower"].iloc[0] < 0.01)
True

References

chen1999monte

mcmc_ess

mcmc_ess(chain: ChainLike) -> Series

Effective sample size of each parameter of an MCMC chain.

The number of independent draws that would estimate the posterior mean as precisely as the autocorrelated chain does: n var(x) / S(0), with S(0) the spectral density at zero of an AR model fitted to the chain (order by AIC). Same definition as R coda::effectiveSize.

Parameters:

Name	Type	Description	Default
chain	array-like, DataFrame or fitted Bayesian result	Draws, one row per iteration and one column per parameter.	required

Returns:

Type	Description
Series	Effective sample size by parameter.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> x = np.random.default_rng(1).normal(size=5000)
>>> bool(4000 < sp.mcmc_ess(x).iloc[0] < 6000)
True

References

plummer2006coda

mcmc_summary

mcmc_summary(chain: ChainLike, quantiles: Sequence[float] = (0.025, 0.25, 0.5, 0.75, 0.975), hpd: Optional[float] = None) -> DataFrame

Posterior summary of an MCMC chain.

Mean, standard deviation, the naive standard error of the mean (sd / sqrt(n), valid only for independent draws), the time-series standard error (which accounts for autocorrelation), the effective sample size and quantiles. The first four columns and the quantiles are those of R summary(coda::mcmc(x)).

Parameters:

Name	Type	Description	Default
chain	array-like, DataFrame or fitted Bayesian result	Draws, one row per iteration and one column per parameter.	required
quantiles	sequence of float	Posterior quantiles to report.	(0.025, 0.25, 0.5, 0.75, 0.975)
hpd	float	If given, also report the highest posterior density interval with this mass (columns hpd_lower, hpd_upper).	None

Returns:

Type	Description
DataFrame	One row per parameter.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> draws = np.random.default_rng(3).normal(size=(4000, 2))
>>> list(sp.mcmc_summary(draws).columns[:5])
['mean', 'sd', 'naive_se', 'ts_se', 'ess']

References

plummer2006coda

raftery_diag

raftery_diag(chain: ChainLike, q: float = 0.025, r: float = 0.005, s: float = 0.95, converge_eps: float = 0.001) -> MCMCDiagnostic

Raftery and Lewis (1992) run-length diagnostic.

How long a chain must be to estimate the posterior quantile q to within +/- r with probability s. The chain is reduced to the indicator of being below its q quantile, thinned until that indicator behaves as a first-order Markov chain, and the two-state chain's transition probabilities give the burn-in M and the total run N (burn-in included). The dependence factor I = N / Nmin compares the total with the length an independent sample would need; values above 5 indicate strong autocorrelation. Same output as R coda::raftery.diag, which prints I to three significant digits.

Parameters:

Name	Type	Description	Default
chain	array-like, DataFrame or fitted Bayesian result		required
q	float	Quantile to be estimated.	0.025
r	float	Desired margin of error of the estimated cumulative probability.	0.005
s	float	Probability of attaining that margin.	0.95
converge_eps	float	Precision required for the estimated time to convergence.	0.001

Returns:

Type	Description
MCMCDiagnostic	table columns: thin (k), burnin (M), total (N), n_min (Nmin) and dependence_factor (I). passed is True when every dependence factor is at most 5.

Raises:

Type	Description
DataInsufficient	If the chain is shorter than Nmin, the length needed even with independent draws.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> x = np.random.default_rng(5).normal(size=5000)
>>> out = sp.raftery_diag(x)
>>> int(out.table["n_min"].iloc[0])
3746

References

raftery1992how

gp_regress

gp_regress(formula: str, data: DataFrame, kernel: str = 'rbf', ard: bool = True, length_scale: Any = None, signal_var: Optional[float] = None, noise_var: Optional[float] = None, mean: Optional[float] = None, optimize_hyper: bool = True, restarts: int = 4, seed: Optional[int] = None, level: float = 0.95, interpolate: bool = False, likelihood: str = 'reml') -> GPResult

Gaussian process regression.

A nonparametric regression with honest uncertainty bands: the regression function has a Gaussian process prior, and its posterior at any point is normal with a closed-form mean and variance.

Parameters:

Name	Type	Description	Default
formula	str	"y ~ x1 + x2". The regressors enter the kernel; the intercept of the formula is the constant mean of the process.	required
data	DataFrame		required
kernel	('rbf', 'matern52', 'matern32')	'rbf' (squared exponential) gives infinitely smooth functions, the Matern kernels rougher ones.	'rbf'
ard	bool	One length scale per regressor. With False a single length scale is shared, which only makes sense for regressors on a common scale.	True
length_scale	optional	Starting values, or the values used when optimize_hyper=False. Defaults: the standard deviation of each regressor, and half the variance of the outcome for each variance.	None
signal_var	optional	Starting values, or the values used when optimize_hyper=False. Defaults: the standard deviation of each regressor, and half the variance of the outcome for each variance.	None
noise_var	optional	Starting values, or the values used when optimize_hyper=False. Defaults: the standard deviation of each regressor, and half the variance of the outcome for each variance.	None
mean	float	A known constant mean. By default the constant is estimated with a flat prior and its uncertainty enters the bands.	None
optimize_hyper	bool	Maximise the marginal likelihood over the hyperparameters.	True
restarts	int	Additional random starting points; the marginal likelihood can have several local maxima.	4
seed	int	For the random starting points.	None
level	float		0.95
interpolate	bool	The outcome has no noise (the output of a deterministic computer model or simulator with fixed random numbers): the fit passes through the data and the bands shrink to zero at the sample points. The noise variance is fixed at a negligible value that keeps the kernel matrix invertible instead of being estimated.	False
likelihood	('reml', 'ml')	How the constant mean is treated when the hyperparameters are estimated. 'reml' integrates it out under a flat prior; 'ml' replaces it by its estimate (the profile likelihood, the convention of kriging software such as R's rkriging and DiceKriging). They differ in small samples; REML gives somewhat larger variances.	'reml'

Returns:

Type	Description
GPResult	params holds the hyperparameters, predict() the posterior mean and bands of the regression function at any points, expected_improvement() the criterion for choosing the next point when the function is being minimised or maximised.

Notes

The cost grows with the cube of the sample size; a few thousand observations are the practical limit.

The bands condition on the estimated hyperparameters and do not carry their uncertainty, so they are somewhat too narrow in small samples.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = np.sort(rng.uniform(0, 6, 80))
>>> df = pd.DataFrame({"x": x, "y": np.sin(x) + 0.2 * rng.normal(size=80)})
>>> fit = sp.gp_regress("y ~ x", df)
>>> list(fit.params.index)
['length_scale[x]', 'signal_var', 'noise_var', 'mean']
>>> band = fit.predict(pd.DataFrame({"x": [1.5]}))
>>> bool(band["lower"].iloc[0] < np.sin(1.5) < band["upper"].iloc[0])
True

References

rasmussen2005gaussian

bayes_ivreg

bayes_ivreg(formula: str, data: DataFrame, prior_mean: Any = 0.0, prior_var: Any = 1000.0, first_stage_prior_var: Any = 1000.0, sigma_prior: Optional[Tuple[float, Any]] = None, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Bayesian linear instrumental variables by Gibbs sampling.

The joint model of the first stage and the structural equation,

d = z' delta + v, y = beta d + x' gamma + eps, (v, eps) ~ N(0, Sigma),

with z the exogenous regressors and the instruments. The posterior of beta carries the uncertainty of the first stage and of the error covariance. With a strong instrument and weak priors it is close to normal around 2SLS; with a weak instrument it is wide and skewed, and it then depends on the priors.

Parameters:

Name	Type	Description	Default
formula	str	'y ~ x1 + x2 + (d ~ z1 + z2)', the syntax of sp.ivreg: one endogenous regressor d and its excluded instruments in parentheses.	required
data	DataFrame		required
prior_mean	Any	Normal prior of the structural coefficients, in the order exogenous terms then the endogenous regressor (scalar, vector or matrix as in :func:statspai.bayes_regress). Default variance 1000.	0.0
prior_var	Any	Normal prior of the structural coefficients, in the order exogenous terms then the endogenous regressor (scalar, vector or matrix as in :func:statspai.bayes_regress). Default variance 1000.	0.0
first_stage_prior_var	(float, array or matrix)	Prior variance of the first-stage coefficients (mean zero).	1000
sigma_prior	(df, scale)	Sigma ~ InvWishart(df, scale) with scale a number (times the identity) or a 2 x 2 matrix. Default (3, 0.02), which adds 0.02 to sums of squared residuals and so is weak on any ordinary scale. R bayesm::rivGibbs defaults to an identity scale, sigma_prior=(3, 1.0).	None
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
chains	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000

Returns:

Type	Description
BayesRegressResult	model='iv'. Parameters: the structural coefficients, then rho (the correlation of the two errors: zero means the regressor is exogenous), sigma2_y and sigma2_d (the error variances of the structural equation and of the first stage) and the first-stage coefficients, prefixed fs:. model_info['first_stage_F'] is the usual F statistic of the excluded instruments.

Notes

A statspai.AssumptionWarning is raised when the first-stage F is below 10. The posterior is still the correct posterior of this model, but it is far from normal, its mean may not exist, and it is sensitive to the priors; report the median and the interval, and how they move with prior_var.

fit.prob("rho > 0") is the posterior probability of positive selection; an interval for rho that excludes zero is the Bayesian counterpart of a Hausman test.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> z, v = rng.normal(size=n), rng.normal(size=n)
>>> d = z + v
>>> y = 1 + 0.5 * d + 0.8 * v + 0.6 * rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "d": d, "z": z})
>>> fit = sp.bayes_ivreg("y ~ (d ~ z)", df, draws=2000, burnin=500, seed=1)
>>> list(fit.params.index)[:3]
['Intercept', 'd', 'rho']
>>> bool(fit.prob("rho > 0") > 0.99)
True

References

rossi2005bayesian, ramirezhassan2026introduction

bayes_mixed

bayes_mixed(formula: str, data: DataFrame, group: str, random: Optional[Sequence[str]] = None, family: str = 'normal', random_intercept: bool = True, prior_mean: Any = 0.0, prior_var: Any = None, sigma2_prior: Tuple[float, float] = (0.001, 0.001), re_prior: Optional[Tuple[float, Any]] = None, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesMixedResult

Bayesian hierarchical model for longitudinal / panel data.

y_it = x_it' beta + w_it' b_i + e_it, where the random effects of group i are b_i ~ N(0, D). The fixed effects have a normal prior and D an inverse-Wishart prior. The random effects are assumed independent of the regressors (a random-effects model, not a fixed-effects one).

Parameters:

Name	Type	Description	Default
formula	str	Fixed part, 'y ~ x1 + x2'.	required
data	DataFrame		required
group	str	Column identifying the groups (units of the panel).	required
random	list of str	Columns with a random slope. They should also appear in the formula, so that the random slope is a deviation from the fixed one.	None
family	('normal', 'logit', 'poisson')	'normal': Gaussian errors with variance sigma2; blocked Gibbs sampler of Chib and Carlin (1999). 'logit' / 'poisson': the usual generalised linear mixed model, sampled by Metropolis steps for the fixed and the random effects.	'normal'
random_intercept	bool		True
prior_mean	Any	Normal prior of the fixed effects, as in :func:statspai.bayes_regress. Default variance 1000 (100 for logit).	0.0
prior_var	Any	Normal prior of the fixed effects, as in :func:statspai.bayes_regress. Default variance 1000 (100 for logit).	0.0
sigma2_prior	(alpha0, delta0)	sigma2 ~ InvGamma(alpha0 / 2, delta0 / 2), normal family.	(0.001, 0.001)
re_prior	(df, scale)	D ~ InvWishart(df, df * scale), whose mean is df * scale / (df - q - 1). scale may be a number, a vector of variances or a q x q matrix. Default df = q + 1 and scale = 0.02 / df (q the number of random effects): in the full conditional of D the prior then adds 0.02 to the sum of squares of each random effect, which is negligible unless the effects are tiny or the groups very few. R MCMCpack and the textbook use an identity scale, re_prior=(q + 2, 1.0), which asserts variances near one whatever the units of the outcome.	None
draws	int	As in :func:statspai.bayes_regress. Proposals of the Metropolis steps are tuned during the burn-in only.	10000
burnin	int	As in :func:statspai.bayes_regress. Proposals of the Metropolis steps are tuned during the burn-in only.	10000
thin	int	As in :func:statspai.bayes_regress. Proposals of the Metropolis steps are tuned during the burn-in only.	10000
chains	int	As in :func:statspai.bayes_regress. Proposals of the Metropolis steps are tuned during the burn-in only.	10000
seed	int	As in :func:statspai.bayes_regress. Proposals of the Metropolis steps are tuned during the burn-in only.	10000
level	int	As in :func:statspai.bayes_regress. Proposals of the Metropolis steps are tuned during the burn-in only.	10000

Returns:

Type	Description
BayesMixedResult	The reported parameters are the fixed effects, sigma2 (normal family) and the elements of D (var(<term>), cov(<term>,<term>)).

Notes

A prior on D with a scale far from the variance of the random effects is informative, and with few groups even a well-scaled prior matters. The fit measures the prior's share of each variance component (model_info['re_prior_share']) and warns above 25 percent; set re_prior and report how the results move.

R MCMCpack::MCMChlogit and MCMChpoisson add an observation-level normal error to the linear index. The models here do not; they are the ones sp.melogit / sp.mepoisson and Stata melogit / mepoisson fit by maximum likelihood.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> g = np.repeat(np.arange(40), 5)
>>> df = pd.DataFrame({"id": g, "x": rng.normal(size=200)})
>>> u = rng.normal(size=40)[g]
>>> df["y"] = 1 + 0.5 * df["x"] + u + rng.normal(size=200)
>>> fit = sp.bayes_mixed("y ~ x", df, group="id", draws=1000,
...                      burnin=300, seed=1)
>>> bool(fit.prob("x > 0") > 0.99)
True
>>> slopes = sp.bayes_mixed("y ~ x", df, group="id", random=["x"],
...                         draws=1000, burnin=300, seed=1)
>>> list(slopes.re_cov.columns)
['Intercept', 'x']

References

chib1999mcmc, ramirezhassan2026introduction

bayes_mixture

bayes_mixture(formula: str, data: DataFrame, components: Union[int, str] = 2, alpha: float = 1.0, alpha_prior: Optional[Tuple[float, float]] = None, prior_mean: Any = None, prior_scale: Any = None, sigma2_prior: Optional[Tuple[float, float]] = None, max_components: Optional[int] = None, grid: int = 200, draws: int = 5000, burnin: int = 1000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Bayesian mixture of normal regressions, finite or Dirichlet process.

Each observation belongs to an unobserved component with its own coefficients and error variance. "y ~ 1" gives a mixture of normals: a flexible density estimate and a model-based clustering.

Parameters:

Name	Type	Description	Default
formula	str	"y ~ 1" or "y ~ x1 + x2".	required
data	DataFrame		required
components	int or 'dp'	Number of components, or 'dp' for a Dirichlet process mixture in which the number of occupied components is inferred.	2
alpha	float	Finite mixture: the weights are Dirichlet(alpha, ..., alpha). Dirichlet process: the concentration; the prior expected number of clusters grows like alpha * log(n).	1
alpha_prior	(shape, rate)	Dirichlet process only: a Gamma prior for the concentration, which is then sampled and reported.	None
prior_mean	array	Prior mean b0 of every component's coefficients. Default: the pooled least-squares coefficients.	None
prior_scale	float or matrix	B0 in beta_k | sigma_k^2 ~ N(b0, sigma_k^2 B0). Default 100 * n * inv(X'X): a component mean can sit about ten component standard deviations away from the pooled fit.	None
sigma2_prior	(a0, d0)	sigma_k^2 ~ IG(a0 / 2, d0 / 2). Default a0 = 4 and d0 such that the prior mean of a component variance is a quarter of the pooled residual variance.	None
max_components	int	Dirichlet process only: cap on the number of occupied clusters (memory). Default min(n, 50).	None
grid	int	Points of the density estimate ("y ~ 1" only).	200
draws	int	As in :func:statspai.bayes_regress.	5000
burnin	int	As in :func:statspai.bayes_regress.	5000
thin	int	As in :func:statspai.bayes_regress.	5000
chains	int	As in :func:statspai.bayes_regress.	5000
seed	int	As in :func:statspai.bayes_regress.	5000
level	int	As in :func:statspai.bayes_regress.	5000

Returns:

Type

Description

BayesRegressResult

model='mixture'. For a finite mixture the table has, for each component, its weight, coefficients and variance; components are ordered by their first coefficient within every draw. For a Dirichlet process it has the number of occupied clusters (and alpha). model_info holds 'n_clusters': posterior distribution of the number of occupied components, 'similarity': the matrix of posterior probabilities that two observations share a component, 'cluster': a point estimate of the partition, the sampled partition closest to 'similarity' in squared error, 'density': posterior mean and band of the density on a grid ("y ~ 1" only).

Notes

The default prior depends on the data through the pooled fit. It is printed in result.prior; pass the three prior arguments to fix it.

Component labels are not identified. Ordering by the first coefficient is a convention that works when components differ in that coefficient; the similarity matrix, the partition and the density do not depend on labels.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.r_[rng.normal(-2, 0.5, 100), rng.normal(2, 0.5, 100)]
>>> fit = sp.bayes_mixture("y ~ 1", pd.DataFrame({"y": y}),
...                        draws=300, burnin=200, seed=1)
>>> list(fit.params.index)[:2]
['weight[1]', 'weight[2]']
>>> int(np.unique(fit.model_info["cluster"]).size)
2

References

neal2000markov, escobar1995bayesian

bayes_mnprobit

bayes_mnprobit(formula: str, data: DataFrame, prior_mean: Any = 0.0, prior_var: Any = 100.0, sigma_prior: Optional[Tuple[float, Any]] = None, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Bayesian multinomial probit.

An unordered outcome with three or more categories. Each category has a latent utility that is linear in the regressors with normal errors, correlated across categories, and the category with the largest utility is observed. Unlike the multinomial logit the model does not impose independence of irrelevant alternatives.

Parameters:

Name	Type	Description	Default
formula	str	"choice ~ x1 + x2". The regressors describe the decision maker; every non-base category gets its own coefficients.	required
data	DataFrame		required
prior_mean	float or array	Normal prior of the stacked coefficients on the unrestricted latent scale.	0 and 100
prior_var	float or array	Normal prior of the stacked coefficients on the unrestricted latent scale.	0 and 100
sigma_prior	(df, scale)	Inverse-Wishart prior of the covariance of the utility differences. Default (J + 2, (J + 2) I) for J categories.	None
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
chains	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000

Returns:

Type	Description
BayesRegressResult	model='mnprobit'. The first level of the outcome is the base. Parameters <level>:<term> are coefficients of the utility of that level relative to the base, in units of the standard deviation of the first utility difference. var(level) and cov(a,b) are the covariance of the utility differences relative to its first element.

Notes

Utilities are identified up to location and scale. Location is fixed by differencing against the base and scale by dividing each draw by the standard deviation of the first difference, the approach of McCulloch and Rossi (1994). The error covariance is weakly identified in most data: expect wide posteriors for it and a slow chain.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> x = rng.normal(size=n)
>>> u = np.column_stack([np.zeros(n), 0.8 * x, -0.5 * x]) + rng.normal(size=(n, 3))
>>> df = pd.DataFrame({"x": x, "choice": np.array(["a", "b", "c"])[u.argmax(1)]})
>>> fit = sp.bayes_mnprobit("choice ~ x", df, draws=500, burnin=200, seed=1)
>>> list(fit.params.index)[:4]
['b:Intercept', 'b:x', 'c:Intercept', 'c:x']

References

mcculloch1994exact, rossi2005bayesian

bayes_mvprobit

bayes_mvprobit(formulas: Sequence[str], data: DataFrame, prior_mean: Any = 0.0, prior_var: Any = 100.0, sigma_prior: Optional[Tuple[float, Any]] = None, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Bayesian multivariate probit.

Several binary outcomes, one probit equation each, with errors that are correlated across equations. The correlations say whether the outcomes move together beyond what the regressors explain.

Parameters:

Name	Type	Description	Default
formulas	list of str	One formula per binary outcome, ["y1 ~ x", "y2 ~ x + z"]. Outcomes must be coded 0 and 1.	required
data	DataFrame		required
prior_mean	float or array	Normal prior of the stacked coefficients on the unrestricted latent scale.	0 and 100
prior_var	float or array	Normal prior of the stacked coefficients on the unrestricted latent scale.	0 and 100
sigma_prior	(df, scale)	Inverse-Wishart prior of the unrestricted error covariance. Default (M + 3, (M + 3) I), centred on unit variances and no correlation.	None
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
chains	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000

Returns:

Type	Description
BayesRegressResult	model='mvprobit'. Parameters <outcome>:<term> are probit coefficients (error variance one) and corr(a,b) the error correlations.

Notes

The sampler does not restrict the error covariance to a correlation matrix. Coefficients and covariance are drawn on an arbitrary scale and each draw is normalised afterwards, which is the approach of Rossi, Allenby and McCulloch. The prior on the reported quantities is the one implied by the stated prior.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> x = rng.normal(size=n)
>>> e = rng.multivariate_normal([0, 0], [[1, 0.6], [0.6, 1]], size=n)
>>> df = pd.DataFrame({"x": x, "y1": (0.5 * x + e[:, 0] > 0) * 1,
...                    "y2": (-0.5 * x + e[:, 1] > 0) * 1})
>>> fit = sp.bayes_mvprobit(["y1 ~ x", "y2 ~ x"], df, draws=500, burnin=200, seed=1)
>>> list(fit.params.index)
['y1:Intercept', 'y1:x', 'y2:Intercept', 'y2:x', 'corr(y1,y2)']

References

rossi2005bayesian, albert1993bayesian

bayes_regress

bayes_regress(formula: str, data: DataFrame, model: str = 'normal', prior_mean: Any = None, prior_var: Any = None, sigma2_prior: Optional[Tuple[float, float]] = None, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95, tune: Optional[float] = None, quantile: float = 0.5, scale: Optional[float] = None, scale_prior: Tuple[float, float] = (0.001, 0.001), lower: Optional[float] = None, upper: Optional[float] = None, dof: float = 5.0, size_prior: Optional[Tuple[float, float]] = None, cut_prior_var: float = 1.0, inference: str = 'mcmc', prior: Optional[str] = None, offset: Any = None, exposure: Optional[str] = None, trials: Optional[str] = None) -> BayesRegressResult

Bayesian regression by MCMC.

Fits one of ten likelihoods with a normal prior on the coefficients and returns the posterior draws with their summary, convergence diagnostics and marginal likelihood.

Parameters:

Name	Type	Description	Default
formula	str	'y ~ x1 + x2'.	required
data	DataFrame		required
model	str	'normal' Gaussian linear model, independent priors beta ~ N(b0, B0), sigma2 ~ IG(alpha0/2, delta0/2); Gibbs. 'conjugate' Gaussian linear model with the conjugate prior beta | sigma2 ~ N(b0, sigma2 B0); closed form, exact marginal likelihood. 't' linear model with Student-t errors of dof degrees of freedom (outlier-robust); Gibbs. 'logit', 'poisson', 'negbin' Metropolis with two proposals shaped by the posterior curvature at the mode, a random walk and a heavy-tailed independence proposal centred at the mode (most draws are then nearly independent). For grouped data the logit takes trials= or an outcome written cbind(successes, failures). 'probit' Albert and Chib (1993) data augmentation. 'ologit' ordered logit (proportional odds), the model of sp.ologit and R polr; same reporting as 'oprobit'; Metropolis on all parameters. 'oprobit' ordered probit: Albert and Chib (1993) data augmentation, the cutpoints updated by a Metropolis step with the latent data integrated out (in the spirit of Cowles 1996). 'tobit' censored regression (lower=, upper=), Chib (1992). 'quantile' asymmetric Laplace likelihood at quantile=, Kozumi and Kobayashi (2011). 'mlogit' multinomial logit, first level of the outcome as the base; coefficients are named <level>:<term>; random-walk Metropolis.	'normal'
prior_mean	float or array	Prior mean of the coefficients (scalar or one per coefficient, in the order of the design matrix). 0 when not given; giving it selects prior='vague'.	None
prior_var	(float, array or matrix)	Prior variance: a scalar, a vector of variances or a full covariance matrix. Under prior='vague' the default is 1000 (1 for 'conjugate', where it multiplies sigma2, and for the scale-free binary and ordered models 100). That default is vague only for regressors on a moderate scale; a note is attached to the result when it is not.	None
sigma2_prior	(alpha0, delta0)	sigma2 ~ InvGamma(alpha0 / 2, delta0 / 2) for the normal, conjugate, t and tobit models. Default (0.001, 0.001); with prior='weakly_informative' the default is instead sigma ~ Exponential(1 / sd(y)).	None
draws	int	Draws kept per chain.	10000
burnin	int	Iterations discarded at the start of each chain.	2000
thin	int	Keep one draw in thin; each chain runs burnin + draws * thin iterations.	1
chains	int	Number of independent chains. With more than one, chains start from dispersed values and the Gelman-Rubin factor is reported.	1
seed	int		None
level	float	Mass of the credible intervals in the summary.	0.95
tune	float	Scale of the Metropolis proposal relative to the posterior curvature. Default 2.38 / sqrt(number of parameters moved).	None
quantile	float	Quantile for model='quantile'.	0.5
scale	float	Fix the scale of the asymmetric Laplace likelihood instead of estimating it (scale=1 is R MCMCpack::MCMCquantreg). With a fixed scale the posterior spread depends on the units of y.	None
scale_prior	(n0, s0)	sigma ~ InvGamma(n0 / 2, s0 / 2) on the estimated scale of the asymmetric Laplace likelihood.	(0.001, 0.001)
lower	float	Censoring points for model='tobit'. Observations at or beyond a point are treated as censored there.	None
upper	float	Censoring points for model='tobit'. Observations at or beyond a point are treated as censored there.	None
dof	float	Degrees of freedom for model='t'.	5
size_prior	(shape, rate)	Gamma prior on the size 1 / alpha of the negative binomial. Default (0.5, 0.1); (1, 1), an Exponential(1), with prior='weakly_informative'.	None
cut_prior_var	float	Prior variance of the log distance between consecutive cutpoints of the ordered probit (mean zero).	1
inference	('mcmc', 'vb')	'vb' (model='normal' only): mean-field variational Bayes, q(beta) q(sigma2) by coordinate ascent. Fast and deterministic; the means are close to the posterior means, the standard deviations are too small when coefficients and variance are dependent a posteriori. draws are then independent draws from the approximation and model_info['elbo'] is the evidence lower bound.	'mcmc'
prior	('vague', 'weakly_informative')	Omitted, it is 'weakly_informative' for the models that have one (normal, logit, probit, poisson, negbin) when prior_mean and prior_var are left at their defaults and inference='mcmc', and 'vague' otherwise. 'vague' uses prior_mean and prior_var as given. 'weakly_informative' scales independent normal priors to the data (models normal, logit, probit, poisson, negbin): standard deviation 2.5 sd(y) / sd(x) for each slope and 2.5 sd(y) for the intercept, centred at mean(y), with the regressors centred; sd(y) is replaced by 1 and mean(y) by 0 outside the Gaussian model. These are the defaults of R rstanarm. They rule out effects that are absurd on the scale of the data and little else, so the result does not depend on the units of the regressors, as it does under a fixed prior_var. prior_mean and prior_var must be left at their defaults. .. versionchanged:: 1.40.0 Omitted, it was 'vague' (a fixed N(0, 1000), N(0, 100) for logit and probit) before 1.40. Pass prior='vague' for the earlier numbers.	'vague'
offset	str or array	A term added to the linear index with coefficient one (logit, Poisson and negative binomial models): a column of data or an array. Predictions for new data need the column form.	None
exposure	str	Column whose logarithm is the offset of a count model.	None
trials	str	model='logit' on grouped data: column with the number of trials of each row, the outcome being the number of successes. posterior_epred is then the success probability and posterior_predict a number of successes.	None

Returns:

Type	Description
BayesRegressResult

Notes

Convergence is checked on every fit: a ConvergenceWarning is raised when the smallest effective sample size is below 100 or the split potential scale reduction factor exceeds 1.05. result.diagnostics() runs the Geweke, Heidelberger-Welch and Raftery-Lewis diagnostics.

The result carries the predictive side of the model: posterior_linpred (linear index), posterior_epred (expected outcome) and posterior_predict (new outcomes, with observation noise) return one row per draw, and log_lik the pointwise log-likelihoods that :func:statspai.loo, :func:statspai.waic and :func:statspai.kfold consume. See also :func:statspai.ppc, :func:statspai.bayes_r2 and :func:statspai.loo_r2.

The ordered probit reports the slopes and J - 1 cutpoints and has no intercept, as sp.oprobit. The prior N(prior_mean, prior_var) applies to the slopes and to minus the first cutpoint.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=300)})
>>> df["y"] = (0.3 + df["x"] + rng.normal(size=300) > 0).astype(int)
>>> fit = sp.bayes_regress("y ~ x", df, model="probit", draws=2000,
...                        burnin=500, seed=1,
...                        prior="weakly_informative")
>>> bool(fit.prob("x > 0") > 0.99)
True
>>> lml = fit.log_marginal_likelihood()

References

albert1993bayesian, cowles1996accelerating, chib1992bayes, chib1995marginal, gelfand1990sampling, gelfand1994bayesian, geweke1999using, kozumi2011gibbs, ramirezhassan2026introduction, gelman2008weakly, gelman2020regression

bayes_shrink

bayes_shrink(formula: str, data: DataFrame, prior: str = 'lasso', lam: Optional[float] = None, lam_prior: Tuple[float, float] = (1.0, 1.0), spike_sd: float = 0.02, slab_sd: float = 1.0, inclusion: float = 0.5, sigma2_prior: Tuple[float, float] = (0.001, 0.001), standardize: bool = True, draws: int = 10000, burnin: int = 2000, thin: int = 1, seed: Optional[int] = None, level: float = 0.95, global_scale: Optional[float] = None, p0: Optional[float] = None, slab_scale: Optional[float] = None, slab_df: float = 4.0, family: str = 'gaussian') -> BayesRegressResult

Linear regression with a shrinkage or a variable-selection prior.

For many candidate regressors, when the aim is prediction or a summary of which ones matter. For model probabilities over subsets see :func:statspai.bma.

Parameters:

Name	Type	Description	Default
formula	str	'y ~ x1 + x2 + ...'; the intercept is kept and not shrunk.	required
data	DataFrame		required
prior	('lasso', 'ssvs', 'horseshoe')	'lasso': beta_j | sigma ~ Laplace(0, sigma / lam). 'ssvs': beta_j ~ (1 - g_j) N(0, spike^2) + g_j N(0, slab^2) with g_j ~ Bernoulli(inclusion). 'horseshoe': beta_j | sigma ~ N(0, sigma^2 tau^2 lam_j^2) with half-Cauchy local scales lam_j ~ C+(0, 1) and global scale tau ~ C+(0, global_scale) (Carvalho, Polson and Scott 2010). Small coefficients are shrunk hard toward zero and large ones almost not at all, which neither the lasso nor a ridge prior does.	'lasso'
lam	float	Lasso penalty. When omitted it is estimated, with the prior lam^2 ~ Gamma(shape, rate) given by lam_prior.	None
lam_prior	(shape, rate)		(1, 1)
spike_sd	float	SSVS: standard deviations of the spike and of the slab, in units of sd(y) per standard deviation of the regressor. A slope smaller than a few spike_sd counts as zero.	0.02
slab_sd	float	SSVS: standard deviations of the spike and of the slab, in units of sd(y) per standard deviation of the regressor. A slope smaller than a few spike_sd counts as zero.	0.02
inclusion	float	SSVS: prior probability of the slab.	0.5
sigma2_prior	(alpha0, delta0)	sigma2 ~ InvGamma(alpha0 / 2, delta0 / 2).	(0.001, 0.001)
standardize	bool	Scale every regressor to unit standard deviation for the prior. Without it the prior treats a regressor in dollars and one in thousands of dollars differently.	True
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000
global_scale	float	Horseshoe: scale of the half-Cauchy prior on tau. Default 1 (the original horseshoe) unless p0 is given.	None
p0	float	Horseshoe: prior guess of the number of coefficients that are far from zero. Sets global_scale = p0 / (p - p0) / sqrt(n), the calibration of Piironen and Vehtari (2017); with many regressors and few expected signals this is far smaller than 1 and keeps the noise coefficients from adding up.	None
slab_scale	float	Horseshoe: regularize it (Piironen and Vehtari 2017). The local scales are replaced by c lam_j / sqrt(c^2 + tau^2 lam_j^2) with c^2 ~ InvGamma(slab_df / 2, slab_df * slab_scale^2 / 2), so a coefficient that escapes the shrinkage still has a Student-t prior of scale slab_scale instead of none at all. In units of the residual standard deviation per standard deviation of the regressor; 2.5 is a wide slab. This is the hs() prior of R rstanarm. It matters when the data identify a coefficient weakly, with collinear regressors or few observations per coefficient, where the plain horseshoe lets it drift.	None
slab_df	float	Degrees of freedom of the slab.	4
family	('gaussian', 'logit')	'logit': a binary outcome with the horseshoe (plain or regularized) on the slopes of a logistic regression, sampled by Polya-Gamma data augmentation. There is no residual variance; the scales are those of the log odds, and p0 is calibrated with a latent standard deviation of 2. Too many regressors for the number of observations separate a binary outcome, where the maximum-likelihood fit does not exist; the horseshoe keeps the fit finite, and with slab_scale keeps the separated coefficients from drifting.	'gaussian'

Returns:

Type

Description

BayesRegressResult

model='lasso' or 'ssvs'. Coefficients are on the original scale. For SSVS table['pip'] is the posterior inclusion probability; for the lasso lam is among the parameters when estimated. For the horseshoe tau is among the parameters, table['shrinkage'] is the posterior mean of each coefficient's shrinkage factor (0: untouched, 1: shrunk to zero) and model_info['m_eff'] the posterior mean of the effective number of unshrunk coefficients. The fit supports posterior_predict, log_lik and so :func:statspai.loo.

Notes

The Bayesian lasso shrinks but never sets a coefficient exactly to zero: its posterior mean is not sparse. Credible intervals after shrinkage are not confidence intervals for the unshrunk effect; for inference on one coefficient with many controls use sp.rlasso_effect or sp.dml.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame(rng.normal(size=(150, 6)), columns=list("abcdef"))
>>> df["y"] = 1 + 2 * df["a"] - df["c"] + rng.normal(size=150)
>>> fit = sp.bayes_shrink("y ~ a + b + c + d + e + f", df, prior="ssvs",
...                       draws=1000, burnin=300, seed=1)
>>> bool(fit.table.loc["a", "pip"] > 0.95)
True
>>> las = sp.bayes_shrink("y ~ a + b + c + d + e + f", df, draws=1000,
...                       burnin=300, seed=1)

References

park2008bayesian, george1993variable, carvalho2010horseshoe, makalic2016simple, piironen2017sparsity, polson2013bayesian

abc

abc(simulate: Callable[..., Any], observed: Any, prior: Any, method: str = 'rejection', n_sim: int = 100000, tol: Optional[float] = None, quantile: float = 0.01, adjust: str = 'none', scale: Any = 'mad', n_synthetic: int = 100, draws: int = 5000, burnin: int = 1000, names: Optional[Sequence[str]] = None, vectorized: bool = False, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Approximate Bayesian computation and Bayesian synthetic likelihood.

Posterior inference for a model that can be simulated but has no tractable likelihood. The data enter only through summary statistics.

Parameters:

Name	Type	Description	Default
simulate	callable	simulate(theta, rng) returns the summary statistics of one data set simulated at the parameter vector theta. With vectorized=True it receives an (N, d) array and returns an (N, s) array.	required
observed	array	The summary statistics of the data.	required
prior	list of scipy.stats distributions, or callable	One frozen distribution per parameter (independent priors), or prior(rng, size) returning a (size, d) array. method='synthetic' needs the first form.	required
method	('rejection', 'synthetic')		'rejection'
n_sim	int	Rejection: number of prior draws.	100000
tol	float	Rejection: accept a draw when the scaled distance between its summaries and the observed ones is at most tol. Default: keep the closest quantile of the draws.	None
quantile	float		0.01
adjust	('none', 'linear')	Rejection: 'linear' regresses the accepted parameters on their summaries with Epanechnikov weights and moves every draw to the observed summaries.	'none'
scale	('mad', 'sd')	Rejection: each summary is divided by its median absolute deviation (or standard deviation, or the given numbers) over the simulations before distances are taken. None leaves the summaries as they are.	'mad'
n_synthetic	int	Synthetic likelihood: simulations per parameter value.	100
draws	int	Synthetic likelihood: length of the Metropolis-Hastings chain.	5000
burnin	int	Synthetic likelihood: length of the Metropolis-Hastings chain.	5000
names	list of str	Parameter names. Default theta1, theta2, ...	None
vectorized	bool		False
seed	int		None
level	float		0.95

Returns:

Type	Description
BayesRegressResult	model='abc'. draws are the accepted (and adjusted) parameters, or the chain. acceptance_rate is the share of prior draws kept, or the Metropolis acceptance rate. model_info reports the tolerance and the scaling used.

Notes

Rejection ABC targets the posterior given that the summaries fall within tol of the observed ones, not the posterior given the data. It equals the posterior given the summaries as tol goes to zero, and the posterior given the data only if the summaries are sufficient. Check sensitivity to quantile.

The regression adjustment is exact when the parameters are linear in the summaries with a spread that does not depend on them, and an extrapolation otherwise.

The synthetic likelihood assumes normal summaries. Its posterior depends on n_synthetic through the noise of the estimated mean and covariance; a few times the number of summaries is the minimum.

Examples:

>>> import numpy as np
>>> from scipy import stats
>>> import statspai as sp
>>> y = np.random.default_rng(0).normal(1.0, 1.0, size=50)
>>> def simulate(theta, rng):
...     return rng.normal(theta[:, 0], 1 / np.sqrt(50))[:, None]
>>> fit = sp.abc(simulate, [y.mean()], [stats.norm(0, 3)], n_sim=20000,
...              vectorized=True, seed=1)
>>> bool(abs(fit.params["theta1"] - y.mean()) < 0.1)
True

References

beaumont2002approximate, wood2010statistical

bayes_sur

bayes_sur(formulas: Sequence[str], data: DataFrame, prior_mean: Any = 0.0, prior_var: Any = 1000.0, sigma_prior: Optional[Tuple[float, Any]] = None, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Bayesian seemingly unrelated regressions.

Several linear equations estimated jointly because their errors are correlated. The posterior of each coefficient uses the information in the other equations' errors, as feasible GLS does, and carries the uncertainty of the error covariance.

Parameters:

Name	Type	Description	Default
formulas	list of str	One formula per equation, ["y1 ~ x1 + x2", "y2 ~ x1 + z"]. The same regressors in every equation give the multivariate regression model.	required
data	DataFrame	Rows with a missing value in any equation are dropped from all.	required
prior_mean	Any	Normal prior of the stacked coefficients (equation by equation, in the order of formulas): scalars, a vector or a matrix.	0.0
prior_var	Any	Normal prior of the stacked coefficients (equation by equation, in the order of formulas): scalars, a vector or a matrix.	0.0
sigma_prior	(df, scale)	Sigma ~ InvWishart(df, scale), scale a number (times the identity) or an M x M matrix. Default (M + 1, 0.02), weak on any ordinary scale; R bayesm::rsurGibbs uses an identity scale.	None
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
chains	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000

Returns:

Type	Description
BayesRegressResult	model='sur'. Parameters <outcome>:<term> for the coefficients, then var(<outcome>), cov(<a>,<b>) and corr(<a>,<b>) for the error covariance.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> e = rng.multivariate_normal([0, 0], [[1, 0.7], [0.7, 1]], size=n)
>>> df = pd.DataFrame({"x": rng.normal(size=n), "z": rng.normal(size=n)})
>>> df["y1"] = 1 + 0.5 * df["x"] + e[:, 0]
>>> df["y2"] = -1 + 0.8 * df["z"] + e[:, 1]
>>> fit = sp.bayes_sur(["y1 ~ x", "y2 ~ z"], df, draws=1000, burnin=300, seed=1)
>>> list(fit.params.index)[:4]
['y1:Intercept', 'y1:x', 'y2:Intercept', 'y2:z']
>>> bool(fit.prob("`corr(y1,y2)` > 0") > 0.99)
True

References

zellner1962efficient, rossi2005bayesian

stochvol

stochvol(y: Union[str, ndarray, Series], data: Optional[DataFrame] = None, demean: bool = True, mu_prior: Tuple[float, float] = (0.0, 100.0), phi_prior: Tuple[float, float] = (5.0, 1.5), sigma_prior: float = 1.0, draws: int = 10000, burnin: int = 2000, thin: int = 1, chains: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BayesRegressResult

Stochastic volatility model for a return series.

The log variance h_t of the returns follows a stationary AR(1) with level mu, persistence phi and volatility of volatility sigma. Unlike GARCH the variance has its own shock, and the whole path of volatilities is estimated with its uncertainty.

Parameters:

Name	Type	Description	Default
y	(array, Series or str)	Returns (not prices), or a column name of data.	required
data	DataFrame		None
demean	bool	Subtract the sample mean first.	True
mu_prior	(mean, variance)	Normal prior of the level of the log variance.	(0, 100)
phi_prior	(a, b)	(phi + 1) / 2 ~ Beta(a, b): prior mean of phi 0.54, most mass on positive persistence.	(5, 1.5)
sigma_prior	float	sigma^2 ~ Gamma(1/2, 1 / (2 * sigma_prior)), equivalently +-sigma ~ N(0, sigma_prior).	1
draws	int	As in :func:statspai.bayes_regress.	10000
burnin	int	As in :func:statspai.bayes_regress.	10000
thin	int	As in :func:statspai.bayes_regress.	10000
chains	int	As in :func:statspai.bayes_regress.	10000
seed	int	As in :func:statspai.bayes_regress.	10000
level	int	As in :func:statspai.bayes_regress.	10000

Returns:

Type	Description
BayesRegressResult	model='stochvol', parameters mu, phi, sigma. model_info['volatility'] is a DataFrame with the posterior mean and interval of the standard deviation exp(h_t / 2) at every date.

Notes

The log chi-square error of log y_t^2 is approximated by a mixture of seven normals; the approximation error is small next to posterior uncertainty in samples of the usual size. Exact zeros in y are offset by a small constant before taking logs, and the result says so.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> h = np.zeros(300)
>>> for t in range(1, 300):
...     h[t] = -1 + 0.95 * (h[t - 1] + 1) + 0.2 * rng.normal()
>>> y = np.exp(h / 2) * rng.normal(size=300)
>>> fit = sp.stochvol(y, draws=500, burnin=300, seed=1)
>>> list(fit.params.index)
['mu', 'phi', 'sigma']
>>> fit.model_info["volatility"].shape
(300, 3)

References

kastner2014ancillarity, ramirezhassan2026introduction

bart

bart(formula: str, data: DataFrame, family: str = 'gaussian', n_trees: int = 200, k: float = 2.0, base: float = 0.95, power: float = 2.0, nu: float = 3.0, q: float = 0.9, numcut: int = 100, max_depth: int = 8, min_leaf: int = 5, draws: int = 1000, burnin: int = 500, thin: int = 1, seed: Optional[int] = None, level: float = 0.95) -> BARTResult

Bayesian additive regression trees.

A flexible regression that finds nonlinearities and interactions on its own and returns posterior bands for the regression function.

Parameters:

Name	Type	Description	Default
formula	str	"y ~ x1 + x2 + C(group)". Regressors enter the trees one column of the design matrix at a time; do not add polynomial or interaction terms.	required
data	DataFrame		required
family	('gaussian', 'binary')	'binary' is a probit link for a 0/1 outcome.	'gaussian'
n_trees	int		200
k	float	The prior puts the regression function within the range of the outcome with probability about 2 Phi(k) - 1. Larger values shrink more.	2
base	float	A node at depth d splits with prior probability base * (1 + d) ** -power.	0.95 and 2
power	float	A node at depth d splits with prior probability base * (1 + d) ** -power.	0.95 and 2
nu	float	Prior of the error variance: scaled inverse chi-square with nu degrees of freedom, placing probability q below the residual variance of a linear fit.	3 and 0.9
q	float	Prior of the error variance: scaled inverse chi-square with nu degrees of freedom, placing probability q below the residual variance of a linear fit.	3 and 0.9
numcut	int	Candidate cutpoints per regressor (midpoints of the observed values, or quantiles when there are more than numcut).	100
max_depth	int	Hard cap on the depth of a tree. The prior makes deeper trees vanishingly rare.	8
min_leaf	int	A split is allowed only if both children keep at least this many observations. Without it trees isolate single observations and the error variance is underestimated.	5
draws	int		1000
burnin	int		1000
thin	int		1000
seed	int		None
level	float		0.95

Returns:

Type	Description
BARTResult	fitted and predict() give the posterior mean and band of the regression function, variable_importance the share of splitting rules per regressor, sigma the error scale.

Notes

The sampler uses birth and death moves only. The bands are for the regression function, not for new observations; add sigma for those.

BART bands tend to be somewhat narrow where the function changes quickly and at the edge of the data. Treat them as approximate.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x1": rng.uniform(size=300), "x2": rng.uniform(size=300)})
>>> df["y"] = 10 * np.sin(np.pi * df.x1 * df.x2) + rng.normal(size=300)
>>> fit = sp.bart("y ~ x1 + x2", df, n_trees=50, draws=200, burnin=200, seed=1)
>>> truth = 10 * np.sin(np.pi * df.x1 * df.x2)
>>> bool(np.sqrt(np.mean((fit.fitted["mean"] - truth) ** 2)) < 1.0)
True

References

chipman2010bart