pybhatlib

Open-source Python library for matrix-based econometric inference: multivariate normal CDF evaluation with analytic gradients, and maximum-likelihood estimation of discrete-choice models (MNP, MORP, MDCEV, MNL).

The numerical core follows the GAUSS BHATLIB library by Chandra Bhat and colleagues (UT Austin / Aptech Systems) and implements Bhat's (2018) analytic MVNCD approximation.

Installation

pip install pybhatlib                  # core (NumPy + SciPy + Numba)
pip install "pybhatlib[torch]"         # add PyTorch backend (optional GPU)

From source (for development, or to run the bundled examples and tutorials):

git clone https://github.com/UMN-Choi-Lab/pybhatlib.git
cd pybhatlib
pip install -e .                  # core (NumPy + SciPy + Numba)
pip install -e ".[torch]"         # add PyTorch backend (optional GPU)
pip install -e ".[dev]"           # add pytest, ruff, mypy
pip install -e ".[all]"           # everything (torch + dev)

Quick Start — Multinomial Probit (MNP)

Sample data: examples/data/TRAVELMODE.csv (3 modes — DA / SR / TR, 1125 observations).

from pybhatlib.models.mnp import MNPModel, MNPControl

model = MNPModel(
    data="examples/data/TRAVELMODE.csv",
    alternatives=["Alt1_ch", "Alt2_ch", "Alt3_ch"],
    availability="none",
    spec={
        "CON_SR": {"Alt1_ch": "sero", "Alt2_ch": "uno",     "Alt3_ch": "sero"},
        "CON_TR": {"Alt1_ch": "sero", "Alt2_ch": "sero",    "Alt3_ch": "uno"},
        "IVTT":   {"Alt1_ch": "IVTT_DA", "Alt2_ch": "IVTT_SR", "Alt3_ch": "IVTT_TR"},
        "OVTT":   {"Alt1_ch": "OVTT_DA", "Alt2_ch": "OVTT_SR", "Alt3_ch": "OVTT_TR"},
        "COST":   {"Alt1_ch": "COST_DA", "Alt2_ch": "COST_SR", "Alt3_ch": "COST_TR"},
    },
    control=MNPControl(iid=True),
)
results = model.fit()
results.summary()

Quick Start — Multivariate Ordered Response Probit (MORP)

from pybhatlib.models.morp import MORPModel, MORPControl

model = MORPModel(
    data=df,                                     # DataFrame or CSV path
    dep_vars=["satisfaction", "recommendation"], # ordinal outcome columns
    spec={
        "income":    {"satisfaction": "income",    "recommendation": "income"},
        "age":       {"satisfaction": "age",       "recommendation": "age"},
        "education": {"satisfaction": "education", "recommendation": "education"},
    },
    n_categories=[3, 3],
    control=MORPControl(iid=True, seed=42),
)
results = model.fit()
results.summary()

A runnable end-to-end example (with morp_ate, morp_predict, morp_predict_category) is at examples/tutorials/t05b_morp_ate_predict.ipynb.

Features

Models

• Multinomial Probit (MNP) — IID, flexible covariance, heteroscedastic-only, random coefficients, mixture-of-normals
• Multivariate Ordered Response Probit (MORP) — multiple ordinal outcomes with shared covariance; per-outcome spec mapping
• Multiple Discrete-Continuous Extreme Value (MDCEV) — traditional (Bhat 2008) and linear outside-good utility specifications, selected via MDCEVControl.utility
• Multinomial Logit (MNL)
• Linear regression (LR) — single continuous outcome, closed-form OLS with classical / White-robust / BHHH standard errors
• Mixed / simulation-based models — MixMNLModel, MNPKerCPModel, MORPFlexModel, MDCEVMixedModel on a shared maximum-simulated-likelihood engine (Halton draws, panel data, normal / log-normal / Yeo-Johnson coefficient transforms). New and under active refinement.

Numerical core

• vecup — vecdup, matdupfull, LDLT decomposition, truncated MVN moments
• matgradient — gradcovcor, gomegxomegax, spherical / Cholesky parameterizations
• gradmvn — Bhat (2018) MVNCD analytic approximation with analytic gradients

Estimation

• Multiple SE estimators (se_method="bhhh" | "hessian" | "sandwich"); BHHH is the default to match GAUSS BHATLIB's _max_CovPar=2. All three are computed at fit time and exposed as se_bhhh / se_hessian / se_sandwich; results.summary() prints a side-by-side diagnostic block (a large Hessian/BHHH divergence is a misspecification signal).
• MNPControl.active_mask to freeze a subset of parameters at their starting values without recoding the model.
• verbose levels: 0 silent, 1 summary, 2 per-iteration NLL, 3 per-iteration parameter / gradient / relative-gradient table.

Backend

• NumPy by default; optional PyTorch backend with GPU support (install with [torch]). All numerical functions take an optional xp kwarg.

Post-estimation

• Average Treatment Effects (ATE) with scenario-matrix support, forecasting, per-category MORP probability prediction.
• MORP ATEs from supplied coefficients without re-fitting: MORPResults.from_estimates(beta, thresholds, correlation) rebuilds a results object from natural-space estimates, and morp_ate_from_params / morp_joint_probs (mean joint category-combination probabilities, the GAUSS ate1.csv equivalent) compute effects directly from them.

Verification

pybhatlib reproduces Table 1 of the BHATLIB GAUSS reference implementation using the TRAVELMODE dataset (3 modes — DA, SR, TR; 1125 observations):

Models (a)–(c) reproduce the GAUSS estimates and BHHH standard errors to ≤0.001 on every parameter (verified end-to-end against GAUSS 26.1.1 + MaxLik 5.0.9). Model (d) is multi-modal — the Python optimum is a slightly better local mode than the GAUSS one.

For MORP, iid=False now uses GAUSS BHATLIB's unit-variance identification by default (MORPControl.fix_scales=True): the latent-utility variances are fixed at 1 and only correlations are estimated. The full-covariance MORP_DINING / MORP_WALK models reproduce the GAUSS mean log-likelihoods (−4.6598 / −3.7591) and correlation matrices. summary() reports the actual threshold cut-points (with delta-method standard errors) and a gradient column, matching GAUSS's output. See docs/plans/MORP_BHATLIB_PARITY.md.

The driving notebooks are t04h_bhatlib_table1.ipynb and t04i_bhat2018_table2.ipynb.

Tutorials

A complete tutorial series lives under examples/tutorials/ as Jupyter notebooks (each with a matching .py script under python_scripts/):

Contributing

See CONTRIBUTING.md for development setup, the PR-test workflow, and CODEOWNERS conventions. Key tests:

pytest tests/                 # full suite
pytest tests/ -m "not slow"   # skip integration tests
pytest tests/ -m torch        # PyTorch backend tests only

References

1. Bhat, C. R. (2018). New Matrix-Based Methods for the Analytic Evaluation of the Multivariate Cumulative Normal Distribution Function. Transportation Research Part B, 109: 238–256.
2. Bhat, C. R. (2008). The Multiple Discrete-Continuous Extreme Value (MDCEV) Model: Role of Utility Function Parameters, Identification Considerations, and Model Extensions. Transportation Research Part B, 42(3): 274–303.

License

MIT