Binary Probit — bayesm::rbprobitGibbs (R)

Albert & Chib (1993) Gibbs sampler — R vs Python comparison

Same synthetic dataset as binary_probit_python.ipynb (data loaded from binary_probit_data.csv; run the Python notebook first).

Structure:

Section	Topic
Section 1	Run rbprobitGibbs, posterior summary and plots
Section 2	Outlier detection: CPO vs Pearson (MLE)

References:

• Albert & Chib (1993), JASA 88(422), 669–679 — data-augmentation Gibbs sampler
• Rossi, Allenby & McCulloch (2005) — Bayesian Statistics and Marketing (bayesm)

# Set library path so IRkernel finds user-installed packages (e.g. bayesm)
.libPaths(c(file.path(Sys.getenv("USERPROFILE"), "R", "win-library", "4.6"),
            .libPaths()))
library(bayesm)

DATA_PATH <- "C:/Users/user/project/binary_probit_data.csv"

dat        <- read.csv(DATA_PATH)
y          <- dat$y
X          <- as.matrix(dat[, c("x0","x1","x2","x3")])
is_outlier <- as.logical(dat$is_outlier)
n          <- nrow(X)
k          <- ncol(X)
beta_true  <- c(0.5, 1.2, -0.8, 0.5)
n_out      <- sum(is_outlier)

cat(sprintf("n=%d  k=%d  outliers=%d  y=1 rate=%.3f\n", n, k, n_out, mean(y)))
cat("bayesm version:", as.character(packageVersion("bayesm")), "\n")
cat("rbprobitGibbs found:", exists("rbprobitGibbs", where="package:bayesm"), "\n")

n=400  k=4  outliers=25  y=1 rate=0.590
bayesm version: 3.1.7 
rbprobitGibbs found: TRUE 

---

1. Run rbprobitGibbs and compare with Python Gibbs

Prior: $\beta \sim N(\bar{\beta}_0,\, A_0^{-1})$ where $A_0$ is the precision matrix. Using $A_0 = 0.01 \cdot I$ — equivalent to covariance $B_0 = 100 \cdot I$ in Python.

Note: rbprobitGibbs returns only β draws — no latent $z_i$. CPO is therefore computed from β alone (harmonic mean over β draws).

Prior <- list(betabar = rep(0, k), A = 0.01 * diag(k))   # A = precision matrix
Mcmc  <- list(R = 15000, keep = 1, nprint = 3000)
Data  <- list(y = y, X = X)

t0  <- proc.time()
out <- rbprobitGibbs(Data=Data, Prior=Prior, Mcmc=Mcmc)
cat(sprintf("rbprobitGibbs done in %.1f sec\n", (proc.time()-t0)["elapsed"]))

betadraw <- out$betadraw[3001:15000, ]   # drop first 3000 as burn-in
cat(sprintf("Kept draws: %d\n", nrow(betadraw)))

 
Starting Gibbs Sampler for Binary Probit Model
   with  400  observations
Table of y Values
y
  0   1 
164 236 
 
Prior Parms: 
betabar
[1] 0 0 0 0
A
     [,1] [,2] [,3] [,4]
[1,] 0.01 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00
[3,] 0.00 0.00 0.01 0.00
[4,] 0.00 0.00 0.00 0.01
 
MCMC parms: 
R=  15000  keep=  1  nprint=  3000
 
 MCMC Iteration (est time to end - min) 
 3000 (0.0)
 6000 (0.0)
 9000 (0.0)
 12000 (0.0)
 15000 (0.0)
 Total Time Elapsed: 0.02 
rbprobitGibbs done in 0.7 sec
Kept draws: 12000

param_names <- c("beta[0]", "beta[1]", "beta[2]", "beta[3]")

cat("\nbayesm posterior summary:\n")
cat(sprintf("%-10s  %7s  %7s  %7s  %7s  %7s  %7s\n",
            "param", "true", "mean", "sd", "q025", "q975", "covered"))
for (j in seq_len(k)) {
    b   <- betadraw[, j]
    lo  <- quantile(b, 0.025)
    hi  <- quantile(b, 0.975)
    cov <- beta_true[j] >= lo & beta_true[j] <= hi
    cat(sprintf("%-10s  %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %7s\n",
                param_names[j], beta_true[j], mean(b), sd(b),
                lo, hi, ifelse(cov, "TRUE", "FALSE")))
}

bayesm posterior summary:
param          true     mean       sd     q025     q975  covered
beta[0]      0.5000   0.2828   0.0708   0.1438   0.4238    FALSE
beta[1]      1.2000   0.6252   0.0740   0.4829   0.7700    FALSE
beta[2]     -0.8000  -0.4725   0.0802  -0.6300  -0.3169    FALSE
beta[3]      0.5000   0.3547   0.0751   0.2101   0.5046     TRUE

par(mfrow=c(2,2), mar=c(4,3,2.5,1))
for (j in seq_len(k)) {
    hist(betadraw[, j], breaks=60, freq=FALSE,
         main=param_names[j], xlab=expression(beta),
         col="steelblue", border=NA)
    abline(v=beta_true[j], col="black", lwd=2, lty=2)
    abline(v=mean(betadraw[, j]), col="red", lwd=1.5)
    legend("topright",
           legend=c(sprintf("true %.2f", beta_true[j]),
                    sprintf("mean %.4f", mean(betadraw[, j]))),
           col=c("black","red"), lty=c(2,1), lwd=c(2,1.5),
           cex=0.75, bty="n")
}
mtext("bayesm rbprobit — posterior distributions  (dashed = true, red = posterior mean)",
      outer=TRUE, side=3, line=-1.2, cex=0.85)

---

2. Outlier detection: CPO vs Pearson (MLE)

CPO via harmonic mean estimator (Albert & Chib 1995): $$\text{CPO}_i \approx \left(\frac{1}{R}\sum_r \frac{1}{p(y_i \mid \beta^{(r)})}\right)^{-1}$$

MLE via glm(family=binomial(link="probit")) for Pearson residuals $r_i = y_i - \hat{p}_i$.

# CPO via harmonic mean
mu_mat  <- X %*% t(betadraw)          # n x R  linear predictor
p_mat   <- pnorm(mu_mat)              # n x R  Phi(x'beta^r)

lik_mat           <- p_mat            # p(y_i=1 | beta^r) for y_i=1 rows
lik_mat[y == 0, ] <- (1 - p_mat)[y == 0, ]   # p(y_i=0 | beta^r) for y_i=0 rows
lik_mat <- pmax(lik_mat, 1e-10)

cpo  <- 1.0 / rowMeans(1.0 / lik_mat)
lpml <- sum(log(pmax(cpo, 1e-15)))
cat(sprintf("LPML = %.2f\n", lpml))

LPML = -214.66

# MLE via glm
df_mle   <- data.frame(y=y, X)
fit_mle  <- glm(y ~ 0 + x0 + x1 + x2 + x3,
                data=df_mle, family=binomial(link="probit"))
p_hat    <- fitted(fit_mle)
pearson  <- y - p_hat

cat("MLE beta:\n")
print(round(coef(fit_mle), 4))
cat(sprintf("Log-lik: %.2f\n", logLik(fit_mle)[1]))

MLE beta:
     x0      x1      x2      x3 
 0.2831  0.6219 -0.4679  0.3534 
Log-lik: -210.33

# Detection table
cat("\nDetection summary -- top-k vs", n_out, "planted outliers\n")
cat(strrep("-", 60), "\n")
cat(sprintf("%9s  %14s  %14s\n", "k flagged", "CPO  TP/Prec", "|Pearson| TP/Prec"))
cat(strrep("-", 60), "\n")

for (k_flag in c(n_out, 2*n_out, 3*n_out)) {
    cpo_idx <- order(cpo)[1:k_flag]
    prs_idx <- order(abs(pearson), decreasing=TRUE)[1:k_flag]
    tp_c <- sum(is_outlier[cpo_idx])
    tp_p <- sum(is_outlier[prs_idx])
    cat(sprintf("  k=%3d      %2d/%d (%.2f)      %2d/%d (%.2f)\n",
                k_flag,
                tp_c, n_out, tp_c/k_flag,
                tp_p, n_out, tp_p/k_flag))
}

Detection summary -- top-k vs 25 planted outliers
------------------------------------------------------------ 
k flagged    CPO  TP/Prec  |Pearson| TP/Prec
------------------------------------------------------------ 
  k= 25      13/25 (0.52)      13/25 (0.52)
  k= 50      20/25 (0.40)      20/25 (0.40)
  k= 75      23/25 (0.31)      23/25 (0.31)

# Scatter: CPO vs |Pearson|
plot(abs(pearson), cpo,
     col  = ifelse(is_outlier, "red", "steelblue"),
     pch  = ifelse(is_outlier, 4, 19),
     lwd  = ifelse(is_outlier, 2, 1),
     cex  = ifelse(is_outlier, 1.4, 0.7),
     xlab = "|Pearson residual|  (MLE-based)",
     ylab = "CPO  (bayesm rbprobitGibbs)",
     main = paste("CPO vs |Pearson residual|",
                  sprintf("  LPML = %.1f", lpml)))
abline(h = sort(cpo)[n_out],              col="purple", lty=2, lwd=1.5)
abline(v = sort(abs(pearson), dec=TRUE)[n_out], col="orange", lty=2, lwd=1.5)
legend("topright",
       legend=c("Normal obs", "Planted outlier",
                "CPO bottom-25 cutoff", "|Pearson| top-25 cutoff"),
       col=c("steelblue","red","purple","orange"),
       pch=c(19,4,NA,NA), lty=c(NA,NA,2,2), lwd=c(1,2,1.5,1.5),
       pt.cex=c(0.8,1.2,NA,NA), cex=0.8, bty="n")

---

Summary

	bayesm::rbprobitGibbs (R)	Python binary_probit_gibbs
Algorithm	Albert & Chib (1993)	Albert & Chib (1993)
Prior	$A_0$ (precision matrix)	$B_0$ (covariance matrix)
Latent $z_i$ draws	not returned	returned → AC (1995) residuals
CPO	from $\beta$ draws only	from $\beta + z$ draws
$\bar{e}_i$ (latent residual)	not available	directly available

Equivalence: both samplers implement the same AC (1993) two-block Gibbs. Posterior means and CPO rankings should be near-identical (MC noise only).

Key difference: the Python implementation stores $z^{(r)}$ draws, enabling the full Albert & Chib (1995) latent residual diagnostics ($\bar{e}_i$, posterior SD of $e_i$). With rbprobitGibbs, only $\beta$ draws are returned.