Bayesian Inference (`wraquant.bayes`)¶

Returns:

ndarray of shape (n_samples, n_pred) – posterior predictive samples. Take np.percentile(result, [2.5, 97.5], axis=0) for a 95% prediction interval.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = np.column_stack([np.ones(100), rng.normal(size=100)])
>>> y = X @ np.array([1.0, 2.0]) + rng.normal(0, 0.5, 100)
>>> X_new = np.column_stack([np.ones(5), rng.normal(size=5)])
>>> ppc = posterior_predictive(y, X, X_new=X_new)
>>> pi = np.percentile(ppc, [2.5, 97.5], axis=0)
>>> print(f"95% PI for first prediction: [{pi[0, 0]:.2f}, {pi[1, 0]:.2f}]")

bayesian_linear_regression(y, X, prior_mean=None, prior_cov=None, a_0=1.0, b_0=1.0, alpha=0.05)[source]¶

Enhanced Bayesian linear regression with Normal-Inverse-Gamma prior.

This is the “textbook” conjugate Bayesian regression where both the coefficients and the noise variance are unknown:

y | X, beta, sigma^2 ~ N(X beta, sigma^2 I) beta | sigma^2 ~ N(m_0, sigma^2 V_0) sigma^2 ~ InvGamma(a_0, b_0)

The posterior is available in closed form as:

beta | sigma^2, y ~ N(m_n, sigma^2 V_n) sigma^2 | y ~ InvGamma(a_n, b_n) beta | y ~ Student-t(2 a_n, m_n, (b_n / a_n) V_n)

When to use this instead of OLS: Use this when you want full uncertainty quantification (credible intervals on coefficients and on sigma^2), when your data set is small and prior information matters, or when you need the marginal likelihood for model comparison (Bayes factors).

Parameters:

y (ndarray) – Response vector of shape (n,).
X (ndarray) – Design matrix of shape (n, k). Include an intercept column yourself if you want one.
prior_mean (ndarray | None, default: None) – Prior mean for beta, shape (k,). Defaults to zeros (agnostic prior).
prior_cov (ndarray | None, default: None) – Prior covariance scale for beta, V_0 of shape (k, k). The actual prior covariance is sigma^2 * prior_cov. Defaults to 100 * I (weakly informative).
a_0 (float, default: 1.0) – Shape parameter for the Inverse-Gamma prior on sigma^2. Default 1.0.
b_0 (float, default: 1.0) – Scale parameter for the Inverse-Gamma prior on sigma^2. Default 1.0.
alpha (float, default: 0.05) – Significance level for credible intervals. Default 0.05 gives 95 % intervals.

Return type:

BayesianLinearRegressionResult

Returns:

BayesianLinearRegressionResult with closed-form posterior analytics, credible intervals, and log marginal likelihood.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = np.column_stack([np.ones(100), rng.normal(size=100)])
>>> y = X @ np.array([1.0, 2.0]) + rng.normal(0, 0.5, 100)
>>> result = bayesian_linear_regression(y, X)
>>> print(result.posterior_mean)  # close to [1, 2]

bayesian_factor_model(X, n_factors=2, n_samples=1000, rng_seed=42)[source]¶

Bayesian factor model via Gibbs sampling (Bayesian PCA).

Estimates a latent factor model of the form:

X = F @ Lambda^T + epsilon

where F is the (n_obs, n_factors) matrix of latent factor scores and Lambda is the (n_variables, n_factors) loading matrix. Both are given conjugate Gaussian priors and sampled via Gibbs.

When to use this: Use a Bayesian factor model instead of standard PCA when you need uncertainty estimates on the loadings and explained variance, especially with short financial time series where the number of observations is not much larger than the number of assets.

Parameters:

X (ndarray) – Data matrix of shape (n_obs, n_variables). Columns are typically asset returns.
n_factors (int, default: 2) – Number of latent factors to estimate.
n_samples (int, default: 1000) – Number of Gibbs samples to draw (after a built-in burn-in of n_samples // 2).
rng_seed (int, default: 42) – Random seed for reproducibility.

Return type:

BayesianFactorModelResult

Returns:

BayesianFactorModelResult with posterior summaries for loadings, scores, and explained variance.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> F = rng.normal(size=(200, 2))
>>> L = rng.normal(size=(5, 2))
>>> X = F @ L.T + rng.normal(0, 0.3, (200, 5))
>>> result = bayesian_factor_model(X, n_factors=2)
>>> print(result.loadings_mean.shape)  # (5, 2)

bayesian_changepoint(data, hazard=0.01, mu_0=0.0, kappa_0=1.0, alpha_0=1.0, beta_0=1.0, threshold=0.3)[source]¶

Bayesian online changepoint detection (Adams & MacKay, 2007).

This algorithm processes observations one at a time and maintains a probability distribution over the run length – how long since the last changepoint. At each new observation the run-length distribution is updated in O(t) time, making the algorithm efficient for streaming data.

The underlying model assumes each segment has data drawn from a Normal distribution with unknown mean and variance, using a Normal-Inverse-Gamma conjugate prior that is re-initialised at each changepoint.

When to use this: Use this for detecting structural breaks in financial time series (regime changes, volatility shifts, mean reversions). Unlike classical CUSUM or Bai-Perron tests, this gives a full posterior probability of a changepoint at each time step rather than a binary yes/no answer.

Parameters:

data (ndarray) – Univariate time series of shape (T,).
hazard (float, default: 0.01) – Constant hazard rate = 1 / (expected run length). hazard = 0.01 means changepoints happen roughly every 100 time steps.
mu_0 (float, default: 0.0) – Prior mean for the Normal component.
kappa_0 (float, default: 1.0) – Prior precision scaling (how confident in mu_0).
alpha_0 (float, default: 1.0) – Shape parameter for the Inverse-Gamma noise prior.
beta_0 (float, default: 1.0) – Scale parameter for the Inverse-Gamma noise prior.
threshold (float, default: 0.3) – Minimum posterior probability to flag a changepoint.

Return type:

BayesianChangepointResult

Returns:

BayesianChangepointResult with run-length probabilities and changepoint posterior.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([rng.normal(0, 1, 100),
...                        rng.normal(5, 1, 100)])
>>> result = bayesian_changepoint(data, hazard=1/50)
>>> cps = result.most_likely_changepoints
>>> # Should find a changepoint near index 100

bayesian_portfolio_bl(returns, views=None, view_confidences=None, P=None, tau=0.05, risk_aversion=2.5, n_samples=5000, rng_seed=42)[source]¶

Black-Litterman model with full Bayesian posterior sampling.

The Black-Litterman model combines a market equilibrium prior (implied by CAPM or equal-weighted) with investor views to produce a posterior distribution of expected returns. This implementation goes beyond the standard BL point estimate: it samples from the full posterior to give uncertainty-aware optimal weights.

When to use this: Use BL when you have subjective views about expected returns or relative performance and want to combine them with market-implied priors. The Bayesian extension is especially useful for constructing robust portfolios with weight confidence intervals.

Parameters:

returns (ndarray) – Historical return matrix, shape (n_periods, n_assets).
views (ndarray | None, default: None) – View vector q, shape (n_views,). If None, the posterior equals the prior (equilibrium returns).
view_confidences (ndarray | None, default: None) – Diagonal of the view uncertainty matrix Omega, shape (n_views,). Smaller values mean more confident views. If None, defaults to tau * diag(P @ Sigma @ P^T).
P (ndarray | None, default: None) – Pick matrix mapping assets to views, shape (n_views, n_assets). Required if views is given.
tau (float, default: 0.05) – Scaling factor for the uncertainty in the prior mean. Typical range 0.01 – 0.10.
risk_aversion (float, default: 2.5) – Risk-aversion parameter delta for the quadratic utility. Default 2.5.
n_samples (int, default: 5000) – Number of posterior samples for weight uncertainty.
rng_seed (int, default: 42) – Random seed for reproducibility.

Return type:

BayesianPortfolioBLResult

Returns:

BayesianPortfolioBLResult with posterior mean/covariance of returns and weight uncertainty.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = rng.normal(0.001, 0.02, (252, 3))
>>> # View: asset 0 will outperform asset 1 by 2 % annualised
>>> P = np.array([[1, -1, 0]])
>>> q = np.array([0.02 / 252])
>>> result = bayesian_portfolio_bl(returns, views=q, P=P)

bayesian_volatility(returns, n_samples=2000, burn_in=1000, rng_seed=42)[source]¶

Bayesian stochastic volatility model via MCMC.

Estimates a time-varying volatility path using the standard stochastic volatility (SV) model:

y_t = exp(h_t / 2) * epsilon_t, epsilon_t ~ N(0, 1) h_t = mu + phi * (h_{t-1} - mu) + eta_t, eta_t ~ N(0, sigma_eta^2)

where h_t is the log-variance at time t, mu is the long-run mean of log-variance, phi is the persistence (typically close to 1 for financial data), and sigma_eta is the volatility-of-volatility.

The model is estimated using a Metropolis-within-Gibbs sampler: the log-volatility path h is sampled block-wise, and the parameters (mu, phi, sigma_eta) are sampled from their conditionals.

When to use this: Use this instead of GARCH when you believe volatility evolves as a latent (unobserved) state rather than being a deterministic function of past returns. The Bayesian approach gives full uncertainty bands on the volatility path.

Parameters:

returns (ndarray) – Return series of shape (T,). Should be de-meaned (or close to zero-mean).
n_samples (int, default: 2000) – Number of MCMC samples to keep after burn-in.
burn_in (int, default: 1000) – Number of initial samples to discard.
rng_seed (int, default: 42) – Random seed.

Return type:

BayesianVolatilityResult

Returns:

BayesianVolatilityResult with posterior volatility path and parameter posteriors.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> vol = np.exp(0.5 * np.cumsum(rng.normal(0, 0.1, 200)))
>>> returns = vol * rng.normal(size=200)
>>> result = bayesian_volatility(returns, n_samples=500, burn_in=200)
>>> print(result.vol_mean.shape)  # (200,)

bayesian_cointegration(y, x, n_samples=5000, rng_seed=42)[source]¶

Bayesian cointegration test between two time series.

Tests whether y and x are cointegrated by:

Estimating the cointegrating regression y = alpha + beta * x + e using Bayesian linear regression (to get a posterior for beta).
Testing the residuals e_t for stationarity using a Bayesian version of the ADF test: Delta e_t = rho * e_{t-1} + u_t. If rho < 0, the residuals are mean-reverting (cointegrated).

The posterior probability of cointegration is estimated as P(rho < 0 | data).

When to use this: Use this for pairs trading (testing if two assets share a long-run equilibrium) or for constructing cointegrated portfolios. The Bayesian approach gives a probability rather than a p-value, which is more natural for decision-making.

Parameters:

y (ndarray) – First time series, shape (T,).
x (ndarray) – Second time series, shape (T,).
n_samples (int, default: 5000) – Number of posterior samples for the ADF coefficient.
rng_seed (int, default: 42) – Random seed.

Return type:

BayesianCointegrationResult

Returns:

BayesianCointegrationResult with the posterior probability of cointegration and the cointegrating vector.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> x = np.cumsum(rng.normal(size=300))
>>> y = 0.8 * x + rng.normal(0, 0.5, 300)  # cointegrated
>>> result = bayesian_cointegration(y, x)
>>> print(f"P(cointegrated) = {result.prob_cointegrated:.2f}")

bayesian_regime_inference(returns, n_regimes=2, n_samples=2000, seed=None)[source]¶

Bayesian regime detection using conjugate priors.

Alternative to frequentist HMM (wraquant.regimes.base.detect_regimes) that provides full posterior uncertainty on regime assignments and transition probabilities. Uses a Gibbs sampler with conjugate Normal-Inverse-Gamma priors for each regime’s mean and variance, and a Dirichlet prior for transition probabilities.

This bridges bayes and regimes by returning a RegimeResult (the standard container from wraquant.regimes.base), so the output can be fed directly into downstream regime-aware portfolio or risk functions.

The generative model is:

s_t | s_{t-1} ~ Categorical(transition_matrix[s_{t-1}]) r_t | s_t ~ N(mu_{s_t}, sigma^2_{s_t})

with priors:

mu_k ~ N(0, prior_var) sigma^2_k ~ InvGamma(alpha_0, beta_0) transition[k] ~ Dirichlet(1, …, 1)

Parameters:

returns (np.ndarray) – Return series, shape (T,).
n_regimes (int, default: 2) – Number of regimes to detect (default 2).
n_samples (int, default: 2000) – Number of Gibbs samples to draw (after a burn-in of n_samples // 2). More samples give smoother posterior estimates.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

states – MAP (most probable) regime at each time step.
probabilities – Posterior regime probabilities (T, K).
transition_matrix – Posterior mean transition matrix.
means – Posterior mean returns per regime.
covariances – Posterior variance per regime.
statistics – Per-regime summary statistics.

Return type:

RegimeResult

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> # Two-regime returns: low vol and high vol
>>> regime = np.concatenate([np.zeros(200), np.ones(200)]).astype(int)
>>> mu = np.array([0.001, -0.002])
>>> sigma = np.array([0.01, 0.03])
>>> returns = rng.normal(mu[regime], sigma[regime])
>>> result = bayesian_regime_inference(returns, n_regimes=2, seed=42)
>>> result.n_regimes
2
>>> result.states.shape
(400,)

Notes

The Gibbs sampler alternates between: 1. Sampling regime states via forward-filtering backward-sampling. 2. Sampling regime parameters (mu, sigma^2) from conjugate posteriors. 3. Sampling transition probabilities from Dirichlet posteriors.

References

Chib (1996), “Calculating posterior distributions and modal estimates in Markov mixture models”
Hamilton (1989), “A new approach to the economic analysis of nonstationary time series”

See also

wraquant.regimes.base.detect_regimes: Frequentist regime detection. wraquant.regimes.base.RegimeResult: Standard result container. bayesian_changepoint: Changepoint detection (different model).

model_comparison(y, X_list, model_names=None)[source]¶

Bayesian model comparison via marginal likelihood and information criteria.

Fits each candidate model using bayesian_linear_regression and computes:

Log marginal likelihood (for Bayes factor computation).
WAIC (Widely Applicable Information Criterion) – a Bayesian analogue of AIC that accounts for the effective number of parameters.
LOO-CV (Leave-One-Out Cross-Validation) approximation via importance sampling.
Bayes factor relative to the best model.

Models are ranked by log marginal likelihood (higher is better).

When to use this: Use this when choosing between competing regression specifications (e.g., which factors belong in a return model). The marginal likelihood automatically penalises model complexity (Occam’s razor), unlike AIC/BIC which use fixed penalties.

Parameters:

y (np.ndarray) – Response vector of shape (n,).
X_list (list[np.ndarray]) – List of design matrices, one per model. Each has shape (n, k_m) where k_m can differ.
model_names (list[str] | None, default: None) – Optional names for each model. Defaults to ['model_0', 'model_1', ...].

Return type:

pd.DataFrame

Returns:

pd.DataFrame with columns log_marginal_likelihood, waic, loo_cv, bayes_factor, rank, sorted by rank (best first).

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> x1 = rng.normal(size=n)
>>> x2 = rng.normal(size=n)
>>> noise = rng.normal(size=n)
>>> y = 1.0 + 2.0 * x1 + rng.normal(0, 0.5, n)
>>> X_good = np.column_stack([np.ones(n), x1])
>>> X_bad = np.column_stack([np.ones(n), x2])
>>> df = model_comparison(y, [X_good, X_bad], ["good", "bad"])
>>> assert df.index[0] == "good"  # good model ranked first

metropolis_hastings(log_posterior, initial, n_samples=10000, proposal_std=1.0, burn_in=1000, thin=1, rng_seed=42)[source]¶

Random-walk Metropolis-Hastings sampler.

The Metropolis-Hastings algorithm is the foundational MCMC method for sampling from an arbitrary (unnormalised) posterior distribution. At each step it proposes a new point from a symmetric Gaussian proposal, accepts or rejects based on the posterior ratio, and thereby generates a Markov chain whose stationary distribution is the target posterior.

Use this when the posterior does not have a convenient closed form and you need posterior samples for uncertainty quantification, credible intervals, or posterior predictive checks.

Tuning: Aim for an acceptance rate of 20–50%. If the rate is too low, reduce proposal_std; if too high, increase it.

Parameters:

log_posterior (Callable[[ndarray], float]) – Function mapping a parameter vector (1-D array) to the log posterior density (up to a normalising constant).
initial (ndarray) – Initial parameter vector of shape (d,).
n_samples (int, default: 10000) – Total number of samples to draw (before burn-in and thinning).
proposal_std (float | ndarray, default: 1.0) – Standard deviation(s) for the isotropic Gaussian proposal. A scalar applies the same std to all dimensions; an array allows per-dimension tuning.
burn_in (int, default: 1000) – Number of initial samples to discard as warm-up.
thin (int, default: 1) – Keep every thin-th sample to reduce autocorrelation.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

Dictionary with keys:

samples – np.ndarray of shape (n_kept, d), posterior samples after burn-in and thinning.
acceptance_rate – float, fraction of proposals accepted.
log_posteriors – np.ndarray, log posterior values for each kept sample.

Return type:

Example

>>> import numpy as np
>>> log_p = lambda q: -0.5 * np.sum((q - 2.0) ** 2)
>>> result = metropolis_hastings(log_p, np.zeros(2),
...     n_samples=5000, burn_in=500, proposal_std=1.0)
>>> result['samples'].shape[1]
2

See also

hamiltonian_monte_carlo: Gradient-based sampler with higher: acceptance rates in high dimensions.

slice_sampler: Tuning-free univariate sampler. convergence_diagnostics: Assess MCMC convergence.

gibbs_sampler(conditionals, initial, n_samples=10000, burn_in=1000, thin=1, rng_seed=42)[source]¶

Gibbs sampler for a model specified via full conditional distributions.

Gibbs sampling is a special case of MCMC where each parameter is updated by drawing from its full conditional distribution (the distribution of that parameter given all others and the data). Because every draw is accepted, Gibbs sampling avoids the tuning burden of Metropolis-Hastings and is efficient when the full conditionals have known forms (e.g., conjugate models).

Use Gibbs when you can analytically derive each conditional (common in Bayesian linear regression, mixture models, and hierarchical models).

Parameters:

conditionals (Sequence[Callable[[ndarray, Generator], float]]) – A list of functions, one per parameter. Each callable has signature fn(current_params, rng) -> float and returns a single draw from the full conditional of that parameter.
initial (ndarray) – Initial parameter vector of shape (d,).
n_samples (int, default: 10000) – Total number of samples to draw (before burn-in and thinning).
burn_in (int, default: 1000) – Number of initial samples to discard.
thin (int, default: 1) – Keep every thin-th sample.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

Posterior samples of shape (n_kept, d).

Return type:

Raises:

ValueError – If the number of conditionals does not match the number of parameters.

Example

>>> import numpy as np
>>> # Gibbs for bivariate normal with rho=0
>>> cond_0 = lambda p, rng: rng.normal(0, 1)
>>> cond_1 = lambda p, rng: rng.normal(0, 1)
>>> samples = gibbs_sampler([cond_0, cond_1], np.zeros(2),
...     n_samples=2000, burn_in=200)
>>> samples.shape[1]
2

See also

metropolis_hastings: General-purpose MCMC when conditionals: are not available.

hamiltonian_monte_carlo: Gradient-based MCMC.

nuts_diagnostic(samples, chains=None)[source]¶

Compute NUTS-style diagnostic statistics (ESS, R-hat, etc.).

Provides the minimum set of diagnostics that every MCMC analysis should report. Can be used on output from any sampler, not just NUTS.

Interpretation:

ESS (Effective Sample Size): How many independent samples your chain is worth, accounting for autocorrelation. If you drew 10,000 samples but ESS = 500, your estimates are only as precise as 500 independent draws. Rule of thumb: ESS > 400 for reliable posterior summaries, ESS > 1000 for tail quantiles.
R-hat: Convergence diagnostic (see gelman_rubin). Only meaningful with multiple chains.
mean/std: Posterior summaries. MCSE = std / sqrt(ESS) gives the precision of the posterior mean estimate.

When to use:

After any MCMC run, before using the samples for inference.
To decide whether to run the sampler longer.

Parameters:

samples (ndarray) – Posterior samples. If 2D (n_samples, n_params), treated as a single chain. If 3D (n_chains, n_samples, n_params), R-hat is computed across chains.
chains (ndarray | None, default: None) – Optional multi-chain samples for R-hat computation. If provided, must be 3D (n_chains, n_samples, n_params).

Returns:

ess (ndarray) – Effective sample size per parameter.
r_hat (ndarray) – R-hat per parameter (NaN if single chain).
mean (ndarray) – Posterior mean per parameter.
std (ndarray) – Posterior standard deviation.

Return type:

dict[str, ndarray]

Example

>>> import numpy as np
>>> samples = np.random.default_rng(0).normal(size=(5000, 3))
>>> diag = nuts_diagnostic(samples)
>>> print(f"ESS: {diag['ess']}")  # Should be ~5000 for iid

trace_summary(samples, param_names=None, quantiles=(0.025, 0.25, 0.5, 0.75, 0.975), chains=None)[source]¶

Compute a publication-ready summary table for MCMC samples.

This is the standard output table that every Bayesian analysis should report. It combines posterior summaries (mean, std, quantiles) with convergence diagnostics (ESS, R-hat) in a single DataFrame.

Interpretation:

mean and std: point estimate and uncertainty.
2.5% and 97.5%: bounds of the 95% credible interval. If this interval excludes 0, the parameter is “significant” in a Bayesian sense.
50% (median): more robust than mean for skewed posteriors.
ESS: effective sample size. If ESS < 400, the quantile estimates (especially the tails) may be unreliable.
R-hat: convergence diagnostic. R-hat > 1.05 is a red flag; do not trust the results.

Parameters:

samples (ndarray) – Posterior samples of shape (n_samples, n_params) or (n_chains, n_samples, n_params).
param_names (Optional[Sequence[str]], default: None) – Parameter names. Defaults to ['param_0', 'param_1', ...].
quantiles (tuple[float, ...], default: (0.025, 0.25, 0.5, 0.75, 0.975)) – Quantiles to compute. Default includes the 95% credible interval endpoints and median.
chains (ndarray | None, default: None) – Optional multi-chain samples for R-hat. If None and samples is 3D, R-hat is computed from the 3D array.

Returns:

mean, std, quantiles, ESS, R-hat. One row per parameter.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> samples = rng.normal(size=(5000, 3))
>>> summary = trace_summary(samples, param_names=['alpha', 'beta', 'sigma'])
>>> print(summary)

gelman_rubin(chains)[source]¶

Compute the Gelman-Rubin (R-hat) convergence diagnostic.

R-hat compares the between-chain variance to the within-chain variance. If chains have converged to the same stationary distribution, these should be approximately equal (R-hat ~ 1). Large R-hat means the chains disagree, indicating they have not yet explored the same region of parameter space.

Interpretation:

R-hat < 1.01: Excellent convergence. Safe to use samples.
R-hat 1.01 - 1.05: Acceptable, but consider running longer.
R-hat 1.05 - 1.1: Marginal. Results may be unreliable.
R-hat > 1.1: NOT converged. Do not use these samples for inference. Run longer, tune the sampler, or reparameterise.

Requires at least 2 chains. Running multiple chains from different starting points is the gold standard for diagnosing convergence.

Parameters:: chains (ndarray | list[ndarray]) – Either a 3D array of shape (n_chains, n_samples, n_params) or a list of 2D arrays each of shape (n_samples, n_params).
Return type:: ndarray
Returns:: R-hat values for each parameter, shape (n_params,).

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> chains = rng.normal(size=(4, 2000, 3))
>>> r_hat = gelman_rubin(chains)
>>> print(f"R-hat: {r_hat}")  # Should be ~1.0

See also

convergence_diagnostics: Full battery including split R-hat,: ESS, and MCSE.

hamiltonian_monte_carlo(log_posterior, grad_log_posterior, initial, n_samples=5000, step_size=0.01, n_leapfrog=20, burn_in=1000, mass_matrix=None, rng_seed=42)[source]¶

Hamiltonian Monte Carlo sampler (pure numpy, no JAX/PyMC needed).

HMC uses Hamiltonian dynamics to propose distant points in parameter space that are likely to be accepted, dramatically reducing the random-walk behaviour of Metropolis-Hastings. This makes it much more efficient for correlated, high-dimensional posteriors.

The algorithm works by augmenting the parameter space with a “momentum” variable and simulating Hamiltonian dynamics using the leapfrog integrator. The total energy (Hamiltonian) is approximately conserved, leading to high acceptance rates even for large jumps.

When to use HMC: Use HMC when Metropolis-Hastings mixes poorly (high autocorrelation, low effective sample size) or when the posterior has strong correlations between parameters. HMC requires the gradient of the log-posterior, which is the main practical barrier.

Tuning: The two key parameters are step_size (too large causes rejections; too small wastes computation) and n_leapfrog (too few gives random-walk behaviour; too many wastes computation). A good heuristic is to aim for acceptance rates of 60–80 %.

Parameters:

log_posterior (Callable[[ndarray], float]) – Function mapping parameter vector to log posterior density (up to a normalising constant).
grad_log_posterior (Callable[[ndarray], ndarray]) – Function mapping parameter vector to the gradient of the log posterior (a vector of the same length).
initial (ndarray) – Initial parameter vector of shape (d,).
n_samples (int, default: 5000) – Number of samples to draw (before burn-in).
step_size (float, default: 0.01) – Leapfrog step size (epsilon). Typical range 0.001 – 0.1 depending on the scale of the problem.
n_leapfrog (int, default: 20) – Number of leapfrog steps per proposal (L).
burn_in (int, default: 1000) – Number of initial samples to discard.
mass_matrix (ndarray | None, default: None) – Diagonal mass matrix of shape (d,). If None, uses identity (unit mass for all parameters). Set this to the inverse of the marginal posterior variances for better performance.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

samples: np.ndarray of shape (n_samples, d) – posterior samples after burn-in.
acceptance_rate: float – fraction of proposals accepted.
log_posteriors: np.ndarray – log posterior values for each kept sample.

Return type:

Example

>>> import numpy as np
>>> # Sample from N(3, 1)
>>> log_p = lambda q: -0.5 * (q[0] - 3.0) ** 2
>>> grad_log_p = lambda q: np.array([-(q[0] - 3.0)])
>>> result = hamiltonian_monte_carlo(log_p, grad_log_p,
...     initial=np.array([0.0]), n_samples=5000, burn_in=500)
>>> print(f"Mean: {result['samples'].mean():.1f}")  # ~3.0

slice_sampler(log_posterior, initial=0.0, n_samples=5000, burn_in=500, w=1.0, rng_seed=42)[source]¶

Univariate slice sampler (Neal, 2003).

Slice sampling is an adaptive MCMC method that requires no tuning of a proposal distribution. It works by:

Drawing a horizontal “slice” under the (unnormalised) density at the current point.
Finding an interval around the current point that contains the slice.
Sampling uniformly from that interval, shrinking it if the proposed point falls outside the slice.

This makes it more robust than Metropolis-Hastings (no rejected samples, no proposal std to tune) while being simple to implement.

When to use this: Use slice sampling for univariate conditional distributions inside a Gibbs sampler, or for simple 1D posteriors where you don’t want to worry about tuning.

Parameters:

log_posterior (Callable[[float], float]) – Function mapping a scalar to the log (unnormalised) posterior density.
initial (float, default: 0.0) – Starting value.
n_samples (int, default: 5000) – Number of samples to keep (after burn-in).
burn_in (int, default: 500) – Number of initial samples to discard.
w (float, default: 1.0) – Initial bracket width. The algorithm adapts from this, so the exact value is not critical. A rough estimate of the posterior std is a good choice.
rng_seed (int, default: 42) – Random seed.

Return type:

Returns:

np.ndarray of shape (n_samples,) – posterior samples.

Example

>>> import numpy as np
>>> # Sample from N(2, 1)
>>> log_p = lambda x: -0.5 * (x - 2.0) ** 2
>>> samples = slice_sampler(log_p, initial=0.0, n_samples=5000)
>>> print(f"Mean: {samples.mean():.1f}")  # ~2.0

convergence_diagnostics(chains, param_names=None)[source]¶

Comprehensive MCMC convergence diagnostics for multi-chain output.

Computes a battery of diagnostics that every Bayesian analysis should report:

Split R-hat: The Gelman-Rubin statistic computed on split chains (each chain is split in half), which is more sensitive to non-stationarity than the original R-hat. Values below 1.01 are considered good; below 1.05 is acceptable.
Effective Sample Size (ESS): The number of effectively independent samples, accounting for autocorrelation. Low ESS means your inferences are based on fewer independent data points than you might think.
ESS per second: Not computed here (would need timing info), but ESS alone tells you if you need to run longer.
Autocorrelation time: The number of steps between effectively independent samples. tau = n / ESS.
MCSE (Monte Carlo Standard Error): The standard error of the posterior mean estimate due to finite sampling. MCSE = posterior_std / sqrt(ESS).

Interpreting the results: Your MCMC has converged if all of: (1) R-hat < 1.05, (2) ESS > 400, (3) trace plots show no trends. If R-hat > 1.1, you need more samples or better tuning.

Parameters:

chains (ndarray | list[ndarray]) – Multi-chain samples, either a 3D array of shape (n_chains, n_samples, n_params) or a list of 2D arrays.
param_names (Optional[Sequence[str]], default: None) – Optional parameter names. Defaults to ['param_0', ...].

Returns:

mean, std, r_hat, split_r_hat, ess, autocorrelation_time, mcse.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> chains = rng.normal(size=(4, 2000, 3))
>>> diag = convergence_diagnostics(chains)
>>> print(diag[['split_r_hat', 'ess']])

Models¶

Pure numpy/scipy Bayesian methods for quantitative finance.

Includes conjugate Bayesian linear regression, Bayesian Sharpe ratio estimation, Bayesian portfolio allocation, Bayesian VaR, credible intervals, Bayes factors, and posterior predictive sampling.

bayesian_regression(y, X, prior_mean=None, prior_cov=None)[source]¶

Conjugate Bayesian linear regression with known noise variance estimate.

The model is y = X @ beta + eps, eps ~ N(0, sigma^2 I), with a normal prior on beta. Sigma^2 is estimated from OLS residuals (plug-in approach).

Interpretation:

posterior_mean is the Bayesian point estimate of beta. With a weakly informative prior, it will be close to OLS.
posterior_cov encodes parameter uncertainty. The diagonal gives the variance of each coefficient; off-diagonals indicate which coefficients are correlated (common in multicollinear factor models).
log_marginal_likelihood is used for model comparison via Bayes factors. Higher values indicate better fit after penalising complexity.

When to use:

Small sample sizes (< 100 obs) where OLS standard errors are unreliable.
When you have informative priors from economic theory.
When you need to compare nested and non-nested models via Bayes factors (see model_comparison).

Parameters:

y (ndarray) – Response vector (n_obs,).
X (ndarray) – Design matrix (n_obs, n_features). Include an intercept column if desired.
prior_mean (ndarray | None, default: None) – Prior mean for beta (n_features,). Defaults to zeros (agnostic prior). Set to economic priors if available (e.g., CAPM beta prior of 1.0).
prior_cov (ndarray | None, default: None) – Prior covariance for beta (n_features, n_features). Defaults to 10 * I (weakly informative). Smaller values express stronger prior conviction.

Returns:

posterior_mean (ndarray) – Posterior mean of beta.
posterior_cov (ndarray) – Posterior covariance of beta.
sigma2 (float) – Estimated noise variance.
log_marginal_likelihood (float) – For model comparison.

Return type:

BayesianRegressionResult

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = np.column_stack([np.ones(100), rng.normal(size=100)])
>>> y = X @ np.array([0.5, 1.2]) + rng.normal(0, 0.3, 100)
>>> result = bayesian_regression(y, X)
>>> print(f"Beta: {result.posterior_mean}")  # close to [0.5, 1.2]

See also

bayesian_linear_regression: Full Normal-Inverse-Gamma conjugate: model that also treats sigma^2 as unknown.
model_comparison: Compare multiple regression models via Bayes: factors.

bayesian_sharpe(returns, prior_mu=0.0, prior_sigma=1.0, n_samples=10000, rng_seed=42)[source]¶

Estimate the Sharpe ratio with full Bayesian uncertainty.

Uses a normal-inverse-gamma conjugate prior for the mean and variance of returns, then computes the posterior distribution of mu / sigma.

Interpretation:

If the 95% credible interval includes 0, the strategy may not have genuine alpha.
Width of the interval indicates estimation uncertainty. A Sharpe of 1.0 with CI [0.2, 1.8] is very different from 1.0 with CI [0.9, 1.1].
prob_positive > 0.95 suggests genuine positive risk-adjusted returns.
prob_positive < 0.8 means you cannot confidently distinguish the strategy from random noise.

When to use:

Track records shorter than 3 years.
When you need confidence in the Sharpe estimate before allocating capital.
For comparing strategies with different track lengths.
When the classical Sharpe looks “too good” on a short sample.

Red flags:

Very wide CI (> 1.5 units): not enough data to judge.
prob_positive between 0.5 and 0.8: inconclusive.
Posterior mean much higher than OLS Sharpe: check the prior.

Parameters:

returns (ndarray) – Return series (daily, weekly, or monthly – the Sharpe ratio is in the same frequency as the returns).
prior_mu (float, default: 0.0) – Prior mean for the return mean. Default 0 (agnostic).
prior_sigma (float, default: 1.0) – Prior std for the return mean. Default 1 (weakly informative). Set lower to express stronger conviction.
n_samples (int, default: 10000) – Number of posterior samples. 10,000 is usually sufficient for stable credible intervals.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

posterior_mean (float) – Posterior mean Sharpe ratio.
posterior_std (float) – Posterior standard deviation.
ci_lower (float) – 2.5th percentile (lower bound).
ci_upper (float) – 97.5th percentile (upper bound).
prob_positive (float) – P(Sharpe > 0). Values > 0.95 suggest genuine positive risk-adjusted returns.
samples (ndarray) – Raw posterior samples for custom analysis (e.g., P(Sharpe > 1)).

Return type:

BayesianSharpeResult

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> daily_returns = rng.normal(0.0003, 0.01, 252)  # ~1 year
>>> result = bayesian_sharpe(daily_returns)
>>> print(f"Sharpe: {result.posterior_mean:.2f} "
...       f"[{result.ci_lower:.2f}, {result.ci_upper:.2f}]")
>>> print(f"P(Sharpe > 0): {result.prob_positive:.1%}")

See also

risk.metrics.sharpe_ratio: Classical point estimate. bayesian_regression: Full Bayesian regression for factor models.

bayesian_portfolio(returns, prior_cov_scale=1.0, n_samples=5000, rng_seed=42)[source]¶

Portfolio allocation with full parameter uncertainty.

Interpretation:

weights_mean is the “average” optimal portfolio across all plausible parameter values. It is typically more diversified than the plug-in MVO solution.
weights_std quantifies allocation uncertainty. Large std means the data does not strongly support a precise allocation to that asset.
If weights_std for an asset is comparable to weights_mean, the allocation is essentially noise – consider removing that asset or using stronger priors.
weight_samples can be used to compute the probability that a weight exceeds a threshold (e.g., P(w_i > 0.3)).

When to use:

Short estimation windows (< 5 years of monthly data).
When MVO gives extreme, unintuitive weights.
When you want to quantify how much you can trust the optimal allocation.

Red flags:

weights_std >> weights_mean: allocation is noise.
Negative expected_return under posterior mean weights: the prior is dragging estimates toward zero.

Parameters:

returns (ndarray) – Return matrix (n_periods, n_assets).
prior_cov_scale (float, default: 1.0) – Scale factor for the prior covariance. Larger values give a more diffuse (less informative) prior, letting the data speak. Default is 1.0.
n_samples (int, default: 5000) – Number of posterior portfolio weight samples. More samples give smoother weight distributions.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

weights_mean (ndarray) – Mean posterior weights.
weights_std (ndarray) – Weight uncertainty per asset.
expected_return (float) – Portfolio return under posterior mean weights.
expected_risk (float) – Portfolio risk under posterior mean weights.
weight_samples (ndarray) – (n_samples, n_assets) raw weight draws for custom analysis.

Return type:

BayesianPortfolioResult

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = rng.normal(0.001, 0.02, (252, 4))
>>> result = bayesian_portfolio(returns, n_samples=2000)
>>> print(f"Weights: {result.weights_mean}")
>>> print(f"Weight uncertainty: {result.weights_std}")

See also

bayesian_portfolio_bl: Black-Litterman model with subjective: views and Bayesian uncertainty.

opt.mean_variance: Classical plug-in MVO.

bayesian_var(returns, confidence=0.95, n_posterior=10000, rng_seed=42)[source]¶

Value-at-Risk with full Bayesian parameter uncertainty.

Interpretation:

var_mean is the posterior expected VaR – your best single estimate.
var_std quantifies estimation uncertainty. If var_std is large relative to var_mean, you need more data or a more informative prior.
ci_upper is the “conservative VaR”: use this for prudent risk management. ci_upper >> var_mean signals high model risk.
var_samples can be used to compute P(VaR > threshold), answering “what is the probability that our true VaR exceeds the regulatory threshold?”

When to use:

Short return histories (< 2 years daily).
When you need to report model risk alongside VaR.
For stress testing: use the upper CI as a conservative VaR.

Red flags:

var_std > 0.5 * var_mean: VaR estimate is unreliable.
ci_upper > 2 * var_mean: extreme model uncertainty.
Negative var_mean: possible with very noisy priors.

Parameters:

returns (ndarray) – Array of observed returns (daily, weekly, etc.).
confidence (float, default: 0.95) – Confidence level (e.g., 0.95 for 95% VaR, 0.99 for 99% VaR). Higher confidence = more extreme losses = more estimation uncertainty.
n_posterior (int, default: 10000) – Number of posterior samples. 10,000 is usually sufficient.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

var_mean (float) – Posterior mean VaR (positive = loss).
var_std (float) – Posterior std of VaR.
ci_lower (float) – 2.5th percentile (optimistic VaR).
ci_upper (float) – 97.5th percentile (conservative VaR – use this for prudent risk management).
var_samples (ndarray) – Raw posterior VaR samples.

Return type:

BayesianVaRResult

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = rng.normal(-0.0001, 0.015, 252)
>>> result = bayesian_var(returns, confidence=0.99)
>>> print(f"VaR: {result.var_mean:.4f} "
...       f"[{result.ci_lower:.4f}, {result.ci_upper:.4f}]")

See also

risk.metrics.value_at_risk: Classical historical VaR. risk.evt: Extreme value theory VaR for tail estimation.

credible_interval(samples, alpha=0.05)[source]¶

Compute the Highest Posterior Density (HPD) credible interval.

Interpretation:

There is a (1 - alpha) posterior probability that the true parameter lies in this interval. This is the statement most practitioners actually want (unlike frequentist confidence intervals).
A narrow HPD interval = precise estimate.
An HPD interval that includes 0 = cannot rule out zero effect.
For symmetric posteriors, HPD equals the equal-tailed interval. For skewed posteriors (common for variances, Sharpe ratios), HPD is always shorter.

Parameters:

samples (ndarray) – Posterior samples (1D array).
alpha (float, default: 0.05) – Significance level. Default 0.05 (95% credible interval).

Return type:

tuple[float, float]

Returns:

(lower, upper) bounds of the HPD interval.

Example

>>> import numpy as np
>>> samples = np.random.default_rng(0).normal(2.0, 0.5, 10000)
>>> lo, hi = credible_interval(samples, alpha=0.05)
>>> print(f"95% HPD: [{lo:.2f}, {hi:.2f}]")  # ~[1.0, 3.0]

bayes_factor(log_likelihood_1, log_likelihood_2)[source]¶

Compute the Bayes factor comparing two models.

BF = exp(log_likelihood_1 - log_likelihood_2).

Interpretation (Kass & Raftery, 1995 scale):

BF < 1/100: decisive evidence for model 2.
BF 1/100 to 1/10: strong evidence for model 2.
BF 1/10 to 1/3: moderate evidence for model 2.
BF 1/3 to 1: weak evidence for model 2.
BF 1 to 3: weak evidence for model 1.
BF 3 to 10: moderate evidence for model 1.
BF 10 to 100: strong evidence for model 1.
BF > 100: decisive evidence for model 1.

When to use:

Comparing nested or non-nested regression models.
Deciding whether an extra factor belongs in a return model.
Testing whether a regime-switching model is justified vs a single-regime model.

Parameters:

log_likelihood_1 (float) – Log marginal likelihood of model 1. Obtain from bayesian_linear_regression(...).log_marginal_likelihood.
log_likelihood_2 (float) – Log marginal likelihood of model 2.

Return type:

float

Returns:

Bayes factor (BF > 1 favors model 1, BF < 1 favors model 2).

Example

>>> # Compare a 3-factor model to a 1-factor model
>>> bf = bayes_factor(log_ml_3factor, log_ml_1factor)
>>> print(f"BF = {bf:.1f}")  # BF > 10 = strong evidence

posterior_predictive(y, X, prior_mean=None, prior_cov=None, n_samples=1000, rng_seed=42, X_new=None)[source]¶

Generate posterior predictive samples for Bayesian linear regression.

Interpretation:

Each row of the output is a plausible future outcome vector under a different posterior draw of beta and sigma.
The spread across rows at a given prediction point is the prediction interval. This is always wider than a confidence interval on the mean because it includes both parameter uncertainty and noise.
If prediction intervals are very wide relative to the response scale, the model has low predictive power.

When to use:

Forecasting returns or factor exposures with uncertainty.
Generating scenario paths for stress testing.
Checking model calibration: are observed values typically within the 95% predictive interval?

Parameters:

y (ndarray) – Response vector (n_obs,).
X (ndarray) – Design matrix (n_obs, n_features).
prior_mean (ndarray | None, default: None) – Prior mean for beta. Defaults to zeros.
prior_cov (ndarray | None, default: None) – Prior covariance for beta. Defaults to 10 * I.
n_samples (int, default: 1000) – Number of posterior predictive samples.
rng_seed (int, default: 42) – Random seed for reproducibility.
X_new (ndarray | None, default: None) – New design matrix for predictions. If None, uses X (in-sample predictions, useful for calibration checks).

Return type:

Returns:

ndarray of shape (n_samples, n_pred) – posterior predictive samples. Take np.percentile(result, [2.5, 97.5], axis=0) for a 95% prediction interval.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = np.column_stack([np.ones(100), rng.normal(size=100)])
>>> y = X @ np.array([1.0, 2.0]) + rng.normal(0, 0.5, 100)
>>> X_new = np.column_stack([np.ones(5), rng.normal(size=5)])
>>> ppc = posterior_predictive(y, X, X_new=X_new)
>>> pi = np.percentile(ppc, [2.5, 97.5], axis=0)
>>> print(f"95% PI for first prediction: [{pi[0, 0]:.2f}, {pi[1, 0]:.2f}]")

bayesian_linear_regression(y, X, prior_mean=None, prior_cov=None, a_0=1.0, b_0=1.0, alpha=0.05)[source]¶

Enhanced Bayesian linear regression with Normal-Inverse-Gamma prior.

This is the “textbook” conjugate Bayesian regression where both the coefficients and the noise variance are unknown:

y | X, beta, sigma^2 ~ N(X beta, sigma^2 I) beta | sigma^2 ~ N(m_0, sigma^2 V_0) sigma^2 ~ InvGamma(a_0, b_0)

The posterior is available in closed form as:

beta | sigma^2, y ~ N(m_n, sigma^2 V_n) sigma^2 | y ~ InvGamma(a_n, b_n) beta | y ~ Student-t(2 a_n, m_n, (b_n / a_n) V_n)

Parameters:

y (ndarray) – Response vector of shape (n,).
X (ndarray) – Design matrix of shape (n, k). Include an intercept column yourself if you want one.
prior_mean (ndarray | None, default: None) – Prior mean for beta, shape (k,). Defaults to zeros (agnostic prior).
prior_cov (ndarray | None, default: None) – Prior covariance scale for beta, V_0 of shape (k, k). The actual prior covariance is sigma^2 * prior_cov. Defaults to 100 * I (weakly informative).
a_0 (float, default: 1.0) – Shape parameter for the Inverse-Gamma prior on sigma^2. Default 1.0.
b_0 (float, default: 1.0) – Scale parameter for the Inverse-Gamma prior on sigma^2. Default 1.0.
alpha (float, default: 0.05) – Significance level for credible intervals. Default 0.05 gives 95 % intervals.

Return type:

BayesianLinearRegressionResult

Returns:

BayesianLinearRegressionResult with closed-form posterior analytics, credible intervals, and log marginal likelihood.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> X = np.column_stack([np.ones(100), rng.normal(size=100)])
>>> y = X @ np.array([1.0, 2.0]) + rng.normal(0, 0.5, 100)
>>> result = bayesian_linear_regression(y, X)
>>> print(result.posterior_mean)  # close to [1, 2]

bayesian_factor_model(X, n_factors=2, n_samples=1000, rng_seed=42)[source]¶

Bayesian factor model via Gibbs sampling (Bayesian PCA).

Estimates a latent factor model of the form:

X = F @ Lambda^T + epsilon

where F is the (n_obs, n_factors) matrix of latent factor scores and Lambda is the (n_variables, n_factors) loading matrix. Both are given conjugate Gaussian priors and sampled via Gibbs.

Parameters:

X (ndarray) – Data matrix of shape (n_obs, n_variables). Columns are typically asset returns.
n_factors (int, default: 2) – Number of latent factors to estimate.
n_samples (int, default: 1000) – Number of Gibbs samples to draw (after a built-in burn-in of n_samples // 2).
rng_seed (int, default: 42) – Random seed for reproducibility.

Return type:

BayesianFactorModelResult

Returns:

BayesianFactorModelResult with posterior summaries for loadings, scores, and explained variance.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> F = rng.normal(size=(200, 2))
>>> L = rng.normal(size=(5, 2))
>>> X = F @ L.T + rng.normal(0, 0.3, (200, 5))
>>> result = bayesian_factor_model(X, n_factors=2)
>>> print(result.loadings_mean.shape)  # (5, 2)

bayesian_changepoint(data, hazard=0.01, mu_0=0.0, kappa_0=1.0, alpha_0=1.0, beta_0=1.0, threshold=0.3)[source]¶

Bayesian online changepoint detection (Adams & MacKay, 2007).

Parameters:

data (ndarray) – Univariate time series of shape (T,).
hazard (float, default: 0.01) – Constant hazard rate = 1 / (expected run length). hazard = 0.01 means changepoints happen roughly every 100 time steps.
mu_0 (float, default: 0.0) – Prior mean for the Normal component.
kappa_0 (float, default: 1.0) – Prior precision scaling (how confident in mu_0).
alpha_0 (float, default: 1.0) – Shape parameter for the Inverse-Gamma noise prior.
beta_0 (float, default: 1.0) – Scale parameter for the Inverse-Gamma noise prior.
threshold (float, default: 0.3) – Minimum posterior probability to flag a changepoint.

Return type:

BayesianChangepointResult

Returns:

BayesianChangepointResult with run-length probabilities and changepoint posterior.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> data = np.concatenate([rng.normal(0, 1, 100),
...                        rng.normal(5, 1, 100)])
>>> result = bayesian_changepoint(data, hazard=1/50)
>>> cps = result.most_likely_changepoints
>>> # Should find a changepoint near index 100

bayesian_portfolio_bl(returns, views=None, view_confidences=None, P=None, tau=0.05, risk_aversion=2.5, n_samples=5000, rng_seed=42)[source]¶

Black-Litterman model with full Bayesian posterior sampling.

Parameters:

returns (ndarray) – Historical return matrix, shape (n_periods, n_assets).
views (ndarray | None, default: None) – View vector q, shape (n_views,). If None, the posterior equals the prior (equilibrium returns).
view_confidences (ndarray | None, default: None) – Diagonal of the view uncertainty matrix Omega, shape (n_views,). Smaller values mean more confident views. If None, defaults to tau * diag(P @ Sigma @ P^T).
P (ndarray | None, default: None) – Pick matrix mapping assets to views, shape (n_views, n_assets). Required if views is given.
tau (float, default: 0.05) – Scaling factor for the uncertainty in the prior mean. Typical range 0.01 – 0.10.
risk_aversion (float, default: 2.5) – Risk-aversion parameter delta for the quadratic utility. Default 2.5.
n_samples (int, default: 5000) – Number of posterior samples for weight uncertainty.
rng_seed (int, default: 42) – Random seed for reproducibility.

Return type:

BayesianPortfolioBLResult

Returns:

BayesianPortfolioBLResult with posterior mean/covariance of returns and weight uncertainty.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> returns = rng.normal(0.001, 0.02, (252, 3))
>>> # View: asset 0 will outperform asset 1 by 2 % annualised
>>> P = np.array([[1, -1, 0]])
>>> q = np.array([0.02 / 252])
>>> result = bayesian_portfolio_bl(returns, views=q, P=P)

bayesian_volatility(returns, n_samples=2000, burn_in=1000, rng_seed=42)[source]¶

Bayesian stochastic volatility model via MCMC.

Estimates a time-varying volatility path using the standard stochastic volatility (SV) model:

y_t = exp(h_t / 2) * epsilon_t, epsilon_t ~ N(0, 1) h_t = mu + phi * (h_{t-1} - mu) + eta_t, eta_t ~ N(0, sigma_eta^2)

The model is estimated using a Metropolis-within-Gibbs sampler: the log-volatility path h is sampled block-wise, and the parameters (mu, phi, sigma_eta) are sampled from their conditionals.

Parameters:

returns (ndarray) – Return series of shape (T,). Should be de-meaned (or close to zero-mean).
n_samples (int, default: 2000) – Number of MCMC samples to keep after burn-in.
burn_in (int, default: 1000) – Number of initial samples to discard.
rng_seed (int, default: 42) – Random seed.

Return type:

BayesianVolatilityResult

Returns:

BayesianVolatilityResult with posterior volatility path and parameter posteriors.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> vol = np.exp(0.5 * np.cumsum(rng.normal(0, 0.1, 200)))
>>> returns = vol * rng.normal(size=200)
>>> result = bayesian_volatility(returns, n_samples=500, burn_in=200)
>>> print(result.vol_mean.shape)  # (200,)

bayesian_cointegration(y, x, n_samples=5000, rng_seed=42)[source]¶

Bayesian cointegration test between two time series.

Tests whether y and x are cointegrated by:

Estimating the cointegrating regression y = alpha + beta * x + e using Bayesian linear regression (to get a posterior for beta).
Testing the residuals e_t for stationarity using a Bayesian version of the ADF test: Delta e_t = rho * e_{t-1} + u_t. If rho < 0, the residuals are mean-reverting (cointegrated).

The posterior probability of cointegration is estimated as P(rho < 0 | data).

Parameters:

y (ndarray) – First time series, shape (T,).
x (ndarray) – Second time series, shape (T,).
n_samples (int, default: 5000) – Number of posterior samples for the ADF coefficient.
rng_seed (int, default: 42) – Random seed.

Return type:

BayesianCointegrationResult

Returns:

BayesianCointegrationResult with the posterior probability of cointegration and the cointegrating vector.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> x = np.cumsum(rng.normal(size=300))
>>> y = 0.8 * x + rng.normal(0, 0.5, 300)  # cointegrated
>>> result = bayesian_cointegration(y, x)
>>> print(f"P(cointegrated) = {result.prob_cointegrated:.2f}")

bayesian_regime_inference(returns, n_regimes=2, n_samples=2000, seed=None)[source]¶

Bayesian regime detection using conjugate priors.

The generative model is:

s_t | s_{t-1} ~ Categorical(transition_matrix[s_{t-1}]) r_t | s_t ~ N(mu_{s_t}, sigma^2_{s_t})

with priors:

mu_k ~ N(0, prior_var) sigma^2_k ~ InvGamma(alpha_0, beta_0) transition[k] ~ Dirichlet(1, …, 1)

Parameters:

returns (np.ndarray) – Return series, shape (T,).
n_regimes (int, default: 2) – Number of regimes to detect (default 2).
n_samples (int, default: 2000) – Number of Gibbs samples to draw (after a burn-in of n_samples // 2). More samples give smoother posterior estimates.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

states – MAP (most probable) regime at each time step.
probabilities – Posterior regime probabilities (T, K).
transition_matrix – Posterior mean transition matrix.
means – Posterior mean returns per regime.
covariances – Posterior variance per regime.
statistics – Per-regime summary statistics.

Return type:

RegimeResult

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> # Two-regime returns: low vol and high vol
>>> regime = np.concatenate([np.zeros(200), np.ones(200)]).astype(int)
>>> mu = np.array([0.001, -0.002])
>>> sigma = np.array([0.01, 0.03])
>>> returns = rng.normal(mu[regime], sigma[regime])
>>> result = bayesian_regime_inference(returns, n_regimes=2, seed=42)
>>> result.n_regimes
2
>>> result.states.shape
(400,)

Notes

References

Chib (1996), “Calculating posterior distributions and modal estimates in Markov mixture models”
Hamilton (1989), “A new approach to the economic analysis of nonstationary time series”

See also

wraquant.regimes.base.detect_regimes: Frequentist regime detection. wraquant.regimes.base.RegimeResult: Standard result container. bayesian_changepoint: Changepoint detection (different model).

model_comparison(y, X_list, model_names=None)[source]¶

Bayesian model comparison via marginal likelihood and information criteria.

Fits each candidate model using bayesian_linear_regression and computes:

Log marginal likelihood (for Bayes factor computation).
WAIC (Widely Applicable Information Criterion) – a Bayesian analogue of AIC that accounts for the effective number of parameters.
LOO-CV (Leave-One-Out Cross-Validation) approximation via importance sampling.
Bayes factor relative to the best model.

Models are ranked by log marginal likelihood (higher is better).

Parameters:

y (np.ndarray) – Response vector of shape (n,).
X_list (list[np.ndarray]) – List of design matrices, one per model. Each has shape (n, k_m) where k_m can differ.
model_names (list[str] | None, default: None) – Optional names for each model. Defaults to ['model_0', 'model_1', ...].

Return type:

pd.DataFrame

Returns:

pd.DataFrame with columns log_marginal_likelihood, waic, loo_cv, bayes_factor, rank, sorted by rank (best first).

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> x1 = rng.normal(size=n)
>>> x2 = rng.normal(size=n)
>>> noise = rng.normal(size=n)
>>> y = 1.0 + 2.0 * x1 + rng.normal(0, 0.5, n)
>>> X_good = np.column_stack([np.ones(n), x1])
>>> X_bad = np.column_stack([np.ones(n), x2])
>>> df = model_comparison(y, [X_good, X_bad], ["good", "bad"])
>>> assert df.index[0] == "good"  # good model ranked first

MCMC¶

MCMC sampling utilities using pure numpy/scipy.

Includes Metropolis-Hastings sampler, Gibbs sampler, NUTS diagnostics (ESS, R-hat), trace summary tables, and the Gelman-Rubin convergence diagnostic.

metropolis_hastings(log_posterior, initial, n_samples=10000, proposal_std=1.0, burn_in=1000, thin=1, rng_seed=42)[source]¶

Random-walk Metropolis-Hastings sampler.

Use this when the posterior does not have a convenient closed form and you need posterior samples for uncertainty quantification, credible intervals, or posterior predictive checks.

Tuning: Aim for an acceptance rate of 20–50%. If the rate is too low, reduce proposal_std; if too high, increase it.

Parameters:

log_posterior (Callable[[ndarray], float]) – Function mapping a parameter vector (1-D array) to the log posterior density (up to a normalising constant).
initial (ndarray) – Initial parameter vector of shape (d,).
n_samples (int, default: 10000) – Total number of samples to draw (before burn-in and thinning).
proposal_std (float | ndarray, default: 1.0) – Standard deviation(s) for the isotropic Gaussian proposal. A scalar applies the same std to all dimensions; an array allows per-dimension tuning.
burn_in (int, default: 1000) – Number of initial samples to discard as warm-up.
thin (int, default: 1) – Keep every thin-th sample to reduce autocorrelation.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

Dictionary with keys:

samples – np.ndarray of shape (n_kept, d), posterior samples after burn-in and thinning.
acceptance_rate – float, fraction of proposals accepted.
log_posteriors – np.ndarray, log posterior values for each kept sample.

Return type:

Example

>>> import numpy as np
>>> log_p = lambda q: -0.5 * np.sum((q - 2.0) ** 2)
>>> result = metropolis_hastings(log_p, np.zeros(2),
...     n_samples=5000, burn_in=500, proposal_std=1.0)
>>> result['samples'].shape[1]
2

See also

hamiltonian_monte_carlo: Gradient-based sampler with higher: acceptance rates in high dimensions.

slice_sampler: Tuning-free univariate sampler. convergence_diagnostics: Assess MCMC convergence.

gibbs_sampler(conditionals, initial, n_samples=10000, burn_in=1000, thin=1, rng_seed=42)[source]¶

Gibbs sampler for a model specified via full conditional distributions.

Use Gibbs when you can analytically derive each conditional (common in Bayesian linear regression, mixture models, and hierarchical models).

Parameters:

conditionals (Sequence[Callable[[ndarray, Generator], float]]) – A list of functions, one per parameter. Each callable has signature fn(current_params, rng) -> float and returns a single draw from the full conditional of that parameter.
initial (ndarray) – Initial parameter vector of shape (d,).
n_samples (int, default: 10000) – Total number of samples to draw (before burn-in and thinning).
burn_in (int, default: 1000) – Number of initial samples to discard.
thin (int, default: 1) – Keep every thin-th sample.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

Posterior samples of shape (n_kept, d).

Return type:

Raises:

ValueError – If the number of conditionals does not match the number of parameters.

Example

>>> import numpy as np
>>> # Gibbs for bivariate normal with rho=0
>>> cond_0 = lambda p, rng: rng.normal(0, 1)
>>> cond_1 = lambda p, rng: rng.normal(0, 1)
>>> samples = gibbs_sampler([cond_0, cond_1], np.zeros(2),
...     n_samples=2000, burn_in=200)
>>> samples.shape[1]
2

See also

metropolis_hastings: General-purpose MCMC when conditionals: are not available.

hamiltonian_monte_carlo: Gradient-based MCMC.

nuts_diagnostic(samples, chains=None)[source]¶

Compute NUTS-style diagnostic statistics (ESS, R-hat, etc.).

Provides the minimum set of diagnostics that every MCMC analysis should report. Can be used on output from any sampler, not just NUTS.

Interpretation:

ESS (Effective Sample Size): How many independent samples your chain is worth, accounting for autocorrelation. If you drew 10,000 samples but ESS = 500, your estimates are only as precise as 500 independent draws. Rule of thumb: ESS > 400 for reliable posterior summaries, ESS > 1000 for tail quantiles.
R-hat: Convergence diagnostic (see gelman_rubin). Only meaningful with multiple chains.
mean/std: Posterior summaries. MCSE = std / sqrt(ESS) gives the precision of the posterior mean estimate.

When to use:

After any MCMC run, before using the samples for inference.
To decide whether to run the sampler longer.

Parameters:

samples (ndarray) – Posterior samples. If 2D (n_samples, n_params), treated as a single chain. If 3D (n_chains, n_samples, n_params), R-hat is computed across chains.
chains (ndarray | None, default: None) – Optional multi-chain samples for R-hat computation. If provided, must be 3D (n_chains, n_samples, n_params).

Returns:

ess (ndarray) – Effective sample size per parameter.
r_hat (ndarray) – R-hat per parameter (NaN if single chain).
mean (ndarray) – Posterior mean per parameter.
std (ndarray) – Posterior standard deviation.

Return type:

dict[str, ndarray]

Example

>>> import numpy as np
>>> samples = np.random.default_rng(0).normal(size=(5000, 3))
>>> diag = nuts_diagnostic(samples)
>>> print(f"ESS: {diag['ess']}")  # Should be ~5000 for iid

trace_summary(samples, param_names=None, quantiles=(0.025, 0.25, 0.5, 0.75, 0.975), chains=None)[source]¶

Compute a publication-ready summary table for MCMC samples.

This is the standard output table that every Bayesian analysis should report. It combines posterior summaries (mean, std, quantiles) with convergence diagnostics (ESS, R-hat) in a single DataFrame.

Interpretation:

mean and std: point estimate and uncertainty.
2.5% and 97.5%: bounds of the 95% credible interval. If this interval excludes 0, the parameter is “significant” in a Bayesian sense.
50% (median): more robust than mean for skewed posteriors.
ESS: effective sample size. If ESS < 400, the quantile estimates (especially the tails) may be unreliable.
R-hat: convergence diagnostic. R-hat > 1.05 is a red flag; do not trust the results.

Parameters:

samples (ndarray) – Posterior samples of shape (n_samples, n_params) or (n_chains, n_samples, n_params).
param_names (Optional[Sequence[str]], default: None) – Parameter names. Defaults to ['param_0', 'param_1', ...].
quantiles (tuple[float, ...], default: (0.025, 0.25, 0.5, 0.75, 0.975)) – Quantiles to compute. Default includes the 95% credible interval endpoints and median.
chains (ndarray | None, default: None) – Optional multi-chain samples for R-hat. If None and samples is 3D, R-hat is computed from the 3D array.

Returns:

mean, std, quantiles, ESS, R-hat. One row per parameter.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> samples = rng.normal(size=(5000, 3))
>>> summary = trace_summary(samples, param_names=['alpha', 'beta', 'sigma'])
>>> print(summary)

gelman_rubin(chains)[source]¶

Compute the Gelman-Rubin (R-hat) convergence diagnostic.

Interpretation:

R-hat < 1.01: Excellent convergence. Safe to use samples.
R-hat 1.01 - 1.05: Acceptable, but consider running longer.
R-hat 1.05 - 1.1: Marginal. Results may be unreliable.
R-hat > 1.1: NOT converged. Do not use these samples for inference. Run longer, tune the sampler, or reparameterise.

Requires at least 2 chains. Running multiple chains from different starting points is the gold standard for diagnosing convergence.

Parameters:: chains (ndarray | list[ndarray]) – Either a 3D array of shape (n_chains, n_samples, n_params) or a list of 2D arrays each of shape (n_samples, n_params).
Return type:: ndarray
Returns:: R-hat values for each parameter, shape (n_params,).

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> chains = rng.normal(size=(4, 2000, 3))
>>> r_hat = gelman_rubin(chains)
>>> print(f"R-hat: {r_hat}")  # Should be ~1.0

See also

convergence_diagnostics: Full battery including split R-hat,: ESS, and MCSE.

hamiltonian_monte_carlo(log_posterior, grad_log_posterior, initial, n_samples=5000, step_size=0.01, n_leapfrog=20, burn_in=1000, mass_matrix=None, rng_seed=42)[source]¶

Hamiltonian Monte Carlo sampler (pure numpy, no JAX/PyMC needed).

Parameters:

log_posterior (Callable[[ndarray], float]) – Function mapping parameter vector to log posterior density (up to a normalising constant).
grad_log_posterior (Callable[[ndarray], ndarray]) – Function mapping parameter vector to the gradient of the log posterior (a vector of the same length).
initial (ndarray) – Initial parameter vector of shape (d,).
n_samples (int, default: 5000) – Number of samples to draw (before burn-in).
step_size (float, default: 0.01) – Leapfrog step size (epsilon). Typical range 0.001 – 0.1 depending on the scale of the problem.
n_leapfrog (int, default: 20) – Number of leapfrog steps per proposal (L).
burn_in (int, default: 1000) – Number of initial samples to discard.
mass_matrix (ndarray | None, default: None) – Diagonal mass matrix of shape (d,). If None, uses identity (unit mass for all parameters). Set this to the inverse of the marginal posterior variances for better performance.
rng_seed (int, default: 42) – Random seed for reproducibility.

Returns:

samples: np.ndarray of shape (n_samples, d) – posterior samples after burn-in.
acceptance_rate: float – fraction of proposals accepted.
log_posteriors: np.ndarray – log posterior values for each kept sample.

Return type:

Example

>>> import numpy as np
>>> # Sample from N(3, 1)
>>> log_p = lambda q: -0.5 * (q[0] - 3.0) ** 2
>>> grad_log_p = lambda q: np.array([-(q[0] - 3.0)])
>>> result = hamiltonian_monte_carlo(log_p, grad_log_p,
...     initial=np.array([0.0]), n_samples=5000, burn_in=500)
>>> print(f"Mean: {result['samples'].mean():.1f}")  # ~3.0

slice_sampler(log_posterior, initial=0.0, n_samples=5000, burn_in=500, w=1.0, rng_seed=42)[source]¶

Univariate slice sampler (Neal, 2003).

Slice sampling is an adaptive MCMC method that requires no tuning of a proposal distribution. It works by:

Drawing a horizontal “slice” under the (unnormalised) density at the current point.
Finding an interval around the current point that contains the slice.
Sampling uniformly from that interval, shrinking it if the proposed point falls outside the slice.

This makes it more robust than Metropolis-Hastings (no rejected samples, no proposal std to tune) while being simple to implement.

When to use this: Use slice sampling for univariate conditional distributions inside a Gibbs sampler, or for simple 1D posteriors where you don’t want to worry about tuning.

Parameters:

log_posterior (Callable[[float], float]) – Function mapping a scalar to the log (unnormalised) posterior density.
initial (float, default: 0.0) – Starting value.
n_samples (int, default: 5000) – Number of samples to keep (after burn-in).
burn_in (int, default: 500) – Number of initial samples to discard.
w (float, default: 1.0) – Initial bracket width. The algorithm adapts from this, so the exact value is not critical. A rough estimate of the posterior std is a good choice.
rng_seed (int, default: 42) – Random seed.

Return type:

Returns:

np.ndarray of shape (n_samples,) – posterior samples.

Example

>>> import numpy as np
>>> # Sample from N(2, 1)
>>> log_p = lambda x: -0.5 * (x - 2.0) ** 2
>>> samples = slice_sampler(log_p, initial=0.0, n_samples=5000)
>>> print(f"Mean: {samples.mean():.1f}")  # ~2.0

convergence_diagnostics(chains, param_names=None)[source]¶

Comprehensive MCMC convergence diagnostics for multi-chain output.

Computes a battery of diagnostics that every Bayesian analysis should report:

Split R-hat: The Gelman-Rubin statistic computed on split chains (each chain is split in half), which is more sensitive to non-stationarity than the original R-hat. Values below 1.01 are considered good; below 1.05 is acceptable.
Effective Sample Size (ESS): The number of effectively independent samples, accounting for autocorrelation. Low ESS means your inferences are based on fewer independent data points than you might think.
ESS per second: Not computed here (would need timing info), but ESS alone tells you if you need to run longer.
Autocorrelation time: The number of steps between effectively independent samples. tau = n / ESS.
MCSE (Monte Carlo Standard Error): The standard error of the posterior mean estimate due to finite sampling. MCSE = posterior_std / sqrt(ESS).

Interpreting the results: Your MCMC has converged if all of: (1) R-hat < 1.05, (2) ESS > 400, (3) trace plots show no trends. If R-hat > 1.1, you need more samples or better tuning.

Parameters:

chains (ndarray | list[ndarray]) – Multi-chain samples, either a 3D array of shape (n_chains, n_samples, n_params) or a list of 2D arrays.
param_names (Optional[Sequence[str]], default: None) – Optional parameter names. Defaults to ['param_0', ...].

Returns:

mean, std, r_hat, split_r_hat, ess, autocorrelation_time, mcse.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> chains = rng.normal(size=(4, 2000, 3))
>>> diag = convergence_diagnostics(chains)
>>> print(diag[['split_r_hat', 'ess']])

Integrations¶

External package wrappers for Bayesian analysis.

Functions in this module require the bayes optional dependency group (PyMC, ArviZ, NumPyro) and are guarded by @requires_extra('bayes').

pymc_regression(y, X, samples=2000, chains=2, target_accept=0.9, random_seed=42)[source]¶

Bayesian linear regression using PyMC.

Fits the model y ~ Normal(X @ beta, sigma) with weakly informative priors.

Parameters:

y (ndarray | Series) – Response variable.
X (ndarray | DataFrame) – Design matrix. An intercept column is added automatically.
samples (int, default: 2000) – Number of posterior samples per chain.
chains (int, default: 2) – Number of MCMC chains.
target_accept (float, default: 0.9) – Target acceptance rate for NUTS.
random_seed (int, default: 42) – Random seed for reproducibility.

Returns:

trace: PyMC InferenceData object, coefficients_mean: np.ndarray of posterior mean coefficients, coefficients_std: np.ndarray of posterior std coefficients, sigma_mean: float — posterior mean of noise std, model: PyMC Model object.

Return type:

arviz_summary(trace, var_names=None, hdi_prob=0.94)[source]¶

Generate a summary table from a trace using ArviZ.

Parameters:

trace (Any) – ArviZ InferenceData object or a dict of arrays.
var_names (list[str] | None, default: None) – Variables to include. If None, includes all.
hdi_prob (float, default: 0.94) – Probability mass for the HDI interval. Default is 0.94.

Returns:

Summary table with mean, sd, HDI, ESS, R-hat.

Return type:

numpyro_regression(y, X, samples=2000, warmup=500, chains=1, rng_seed=0)[source]¶

Bayesian linear regression using NumPyro.

Fits the model y ~ Normal(X @ beta, sigma) using NUTS sampling.

Parameters:

y (ndarray | Series) – Response variable.
X (ndarray | DataFrame) – Design matrix. An intercept column is added automatically.
samples (int, default: 2000) – Number of posterior samples.
warmup (int, default: 500) – Number of warmup (burn-in) samples.
chains (int, default: 1) – Number of MCMC chains.
rng_seed (int, default: 0) – Random seed for reproducibility.

Returns:

samples: dict of posterior samples by parameter name, coefficients_mean: np.ndarray of posterior mean coefficients, coefficients_std: np.ndarray of posterior std coefficients, sigma_mean: float — posterior mean of noise std.

Return type:

bambi_regression(formula, data, family='gaussian', draws=1000, chains=2, seed=None)[source]¶

Fit a Bayesian regression using Bambi’s formula interface.

Bambi provides R-style formula syntax for Bayesian models built on PyMC.

Parameters:

formula (str) – Model formula string (e.g., "y ~ x1 + x2").
data (DataFrame) – DataFrame with variables referenced in formula.
family (str, default: 'gaussian') – Likelihood family ("gaussian", "bernoulli", "poisson").
draws (int, default: 1000) – Number of posterior draws per chain.
chains (int, default: 2) – Number of MCMC chains.
seed (int | None, default: None) – Random seed.

Returns:

model: Bambi Model object, trace: InferenceData posterior samples, summary: pd.DataFrame ArviZ summary table.

Return type:

emcee_sample(log_prob_fn, n_walkers, n_dim, n_steps=1000, initial=None, seed=None)[source]¶

Run ensemble MCMC sampling using emcee.

Parameters:

log_prob_fn (Any) – Log probability function accepting a parameter array.
n_walkers (int) – Number of walkers (must be >= 2 * n_dim).
n_dim (int) – Number of parameters.
n_steps (int, default: 1000) – Number of MCMC steps.
initial (ndarray | None, default: None) – Initial walker positions of shape (n_walkers, n_dim). If None, drawn from a standard normal distribution.
seed (int | None, default: None) – Random seed.

Returns:

samples: array of shape (n_walkers * n_steps, n_dim), log_prob: array of log probabilities, acceptance_fraction: mean acceptance fraction.

Return type:

blackjax_nuts(log_prob_fn, initial_position, n_samples=1000, step_size=0.01, seed=None)[source]¶

Run NUTS sampling using BlackJAX (JAX-based).

Parameters:

log_prob_fn (Any) – Log probability function (must be JAX-compatible).
initial_position (Any) – Initial parameter values.
n_samples (int, default: 1000) – Number of samples.
step_size (float, default: 0.01) – NUTS step size.
seed (int | None, default: None) – Random seed integer.

Returns:

samples: array of posterior samples, divergences: number of divergent transitions.

Return type: