Bayesian Tobit — Censored Gaussian Regression

Python · R  ·  Download Gibbs sampler

Model

When an outcome is only observed within a known region — a sensor that saturates, a value below a detection limit, or hours worked that pile up at 0 — ordinary regression is biased. The censored-Gaussian (Tobit) model fixes it: a latent yi∗∼N(xi′β,σ2)y^*_i \sim N(x_i'\beta, \sigma^2) is observed only within a known region, and the likelihood uses the density for exact observations and CDF mass for censored ones. This is the classic econometrics Tobit (Tobin 1958) and the PyMC "Censored Data" example. The new ingredient over the survival examples is left and interval censoring (not just right), handled in one framework.

yi∗∼N(xi′β,σ2),yi=observed only within a known region of yi∗y^*_i \sim N(x_i'\beta, \sigma^2), \qquad y_i = \text{observed only within a known region of } y^*_i
CensoringWhat we knowLikelihood contribution
exacty∗=yy^* = yϕ((y−x′β)/σ)/σ\phi((y-x'\beta)/\sigma)/\sigma — density
lefty∗≤cy^* \le c (Tobit at 0)Φ((c−x′β)/σ)\Phi((c-x'\beta)/\sigma)
righty∗≥cy^* \ge c1−Φ((c−x′β)/σ)1 - \Phi((c-x'\beta)/\sigma)
intervala<y∗≤ba < y^* \le bΦ((b−x′β)/σ)−Φ((a−x′β)/σ)\Phi((b-x'\beta)/\sigma) - \Phi((a-x'\beta)/\sigma)

Sampler — data augmentation (Chib 1992)

Rather than work with the mixed density/CDF likelihood directly, data augmentation (Chib 1992) makes it trivial: impute each censored yi∗y^*_i from a truncated normal on its region, and the completed data is an ordinary linear regression with conjugate updates. It is the same augmentation idea as the ZIP latent indicator and the Albert–Chib probit latent utility — here the latent is the unobserved value itself.

BlockDrawFull conditional
(1)ycens∗∣⋅y^*_{\text{cens}} \mid \cdotTruncated Normal on each observation's censoring region
(2)β∣⋅\beta \mid \cdotNormal conjugate (completed-data linear regression)
(3)σ2∣⋅\sigma^2 \mid \cdotInverse-Gamma conjugate

Notebooks

Section 1 validates on synthetic data: the Tobit recovers the true slopes (2.0, −1.5) where naive OLS attenuates them toward 0 (≈1.25, −0.91) — the classic censoring bias — and the imputed latent density matches the true latent across the whole range, including the censored region below 0 that censoring collapses onto the boundary. The same sampler recovers the truth under interval censoring too. Section 2 fits the classic Mroz (1987) labor-supply data: annual hours worked for 753 married women, left-censored at 0 for the 43% who did not work. The Tobit estimates match the AER::tobit benchmark (education ≈+81\approx +81 hrs/yr, experience strongly positive but concave, a child under 6 ≈−900\approx -900 hours, age negative), with every OLS coefficient visibly shrunk toward zero by ignoring the censoring. Cross-checked three ways — from-scratch augmentation Gibbs, PyMC (pm.Censored), and R (AER::tobit / survival::survreg).

Python Tobit from scratch — data-augmentation Gibbs + censored-normal MLE, latent imputation, Mroz hours (Python) Open → Python PyMC cross-check — pm.Censored, plus a data-augmentation-Gibbs vs marginalised-NUTS mixing comparison (Python) Open → R R cross-check — AER::tobit / survival::survreg (Gaussian, left-censored) (R) Open →

Downloads

tobit_gibbs.py From-scratch Tobit — augmentation Gibbs + censored-normal MLE, all censoring types (NumPy / scipy) mroz_hours.csv Mroz (1987) — 753 married women, annual hours worked (left-censored at 0)

Gibbs Sampler — Source Code

"""
Censored Gaussian regression (the Tobit model) -- from-scratch data-augmentation Gibbs + censored-normal
MLE. The PyMC "Censored Data" example, and the classic econometrics Tobit (Tobin 1958). Handles all
censoring types at once: exact, left, right, interval.

Model
-----
  y*_i ~ N(x_i' beta, sigma^2)   (a latent Gaussian outcome)
  we observe y*_i only within a known region:
     exact   : y_i = y*_i                          -> density   phi((y-xb)/sigma)/sigma
     left     : y*_i <= hi_i  (e.g. Tobit at 0)    -> Phi((hi-xb)/sigma)
     right    : y*_i >= lo_i                        -> 1 - Phi((lo-xb)/sigma)
     interval : lo_i < y*_i <= hi_i                 -> Phi((hi-xb)/sigma) - Phi((lo-xb)/sigma)

Data-augmentation Gibbs (Chib 1992): impute each censored y*_i from a TRUNCATED normal on its region,
then the completed data is an ordinary linear regression -> conjugate Normal/Inverse-Gamma updates.
  1. y*_cens | beta, sigma  ~ TruncNormal(xb, sigma; lo, hi)
  2. beta | y*, sigma       ~ Normal  (conjugate; prior beta ~ N(0, prior_sd^2 I))
  3. sigma^2 | y*, beta     ~ Inverse-Gamma(a0, b0)
"""
import numpy as np
from scipy.stats import truncnorm, norm
from scipy.optimize import minimize


def simulate_tobit(n=2000, beta=(1.0, 2.0, -1.5), sigma=2.0, L=0.0, seed=0):
    """Latent Gaussian regression, then LEFT-censor at L (the classic Tobit)."""
    rng = np.random.default_rng(seed); beta = np.asarray(beta, float); k = len(beta)
    X = np.column_stack([np.ones(n)] + [rng.normal(size=n) for _ in range(k - 1)])
    ystar = X @ beta + sigma * rng.standard_normal(n)
    y = np.maximum(ystar, L)                                   # observed; censored ones sit at L
    exact = ystar > L
    lo = np.where(exact, y, -np.inf); hi = np.where(exact, y, L)
    return y, X, lo, hi, exact, ystar


def _cens_negloglik(par, X, lo, hi, exact):
    k = X.shape[1]; beta = par[:k]; sigma = np.exp(par[k]); xb = X @ beta
    ll = 0.0
    if exact.any():
        z = (hi[exact] - xb[exact]) / sigma                   # for exact, lo==hi==y
        ll += np.sum(norm.logpdf(z) - np.log(sigma))
    c = ~exact
    if c.any():
        Phi_hi = norm.cdf((hi[c] - xb[c]) / sigma)            # Phi(inf)=1
        Phi_lo = norm.cdf((lo[c] - xb[c]) / sigma)            # Phi(-inf)=0
        ll += np.sum(np.log(np.clip(Phi_hi - Phi_lo, 1e-300, None)))
    return -ll


def tobit_mle(y, X, lo, hi, exact):
    k = X.shape[1]
    b0 = np.linalg.lstsq(X, y, rcond=None)[0]
    x0 = np.concatenate([b0, [np.log(np.std(y) + 1e-6)]])
    res = minimize(_cens_negloglik, x0, args=(X, lo, hi, exact), method='Nelder-Mead',
                   options=dict(maxiter=40000, xatol=1e-7, fatol=1e-7))
    beta = res.x[:k]; sigma = np.exp(res.x[k])
    return beta, sigma


def tobit_gibbs(y, X, lo, hi, exact, R=12000, burn=3000, seed=0, prior_sd=1e4, a0=0.01, b0=0.01):
    rng = np.random.default_rng(seed)
    n, k = X.shape; cens = ~exact
    XtX = X.T @ X; P0 = np.eye(k) / prior_sd ** 2
    beta = np.linalg.lstsq(X, y, rcond=None)[0]; sigma = np.std(y) + 1e-6
    ystar = y.astype(float).copy()
    keep = R - burn; B = np.zeros((keep, k)); S = np.zeros(keep); YS = np.zeros(n)
    for it in range(R):
        # 1. impute censored latents from truncated normals
        mu = X @ beta
        a = (lo[cens] - mu[cens]) / sigma; b = (hi[cens] - mu[cens]) / sigma
        ystar[cens] = truncnorm.rvs(a, b, loc=mu[cens], scale=sigma, random_state=rng)
        # 2. beta | y*, sigma   (conjugate Normal)
        prec = XtX / sigma ** 2 + P0
        cov = np.linalg.inv(prec); m = cov @ (X.T @ ystar) / sigma ** 2
        beta = m + np.linalg.cholesky(cov) @ rng.standard_normal(k)
        # 3. sigma^2 | y*, beta  (Inverse-Gamma)
        ssr = np.sum((ystar - X @ beta) ** 2)
        sigma = np.sqrt(1.0 / rng.gamma(a0 + n / 2.0, 1.0 / (b0 + 0.5 * ssr)))
        if it >= burn:
            B[it - burn] = beta; S[it - burn] = sigma; YS += ystar
    return dict(beta=B, sigma=S, ystar=YS / keep)   # ystar = posterior-mean imputed latent per obs


def summary(draws, names):
    d = np.atleast_2d(draws.T).T if draws.ndim == 1 else draws
    if draws.ndim == 1:
        return [(names if isinstance(names, str) else names[0], draws.mean(), draws.std(),
                 *np.percentile(draws, [2.5, 97.5]))]
    return [(names[j], d[:, j].mean(), d[:, j].std(), *np.percentile(d[:, j], [2.5, 97.5]))
            for j in range(d.shape[1])]

References