Constrained Bayesian regression — R companion¶

Companion to constrained_lm_python.ipynb (Geweke 1995). Three R tools map onto the project's pieces, on the same data (constrained_lm_data.csv, n=50, k=4, constraints $\beta_1>0,\ \beta_2>0$, $\sigma^2=2$ known):

project piece R package
GHK probability $p_{2\|1}=P(a\le D\beta\le w)$ TruncatedNormal::pmvnorm (Botev minimax-tilting — more accurate than GHK)
constrained posterior draws (truncated Gibbs) tmvtnorm::rtmvnorm (Gibbs sampler for the truncated MVN)
constraint evidence / hypothesis bain (Bayes factor for inequality-constrained hypotheses)
In [1]:
suppressMessages({library(tmvtnorm); library(TruncatedNormal); library(bain); library(ggplot2)})
d <- read.csv('constrained_lm_data.csv'); X <- cbind(1, d$x1, d$x2, d$x3); y <- d$y
sigma2 <- 2.0                                              # known, matches the Python section-2 design
XtX <- t(X) %*% X; bbar <- solve(XtX, t(X) %*% y); Bbar <- sigma2 * solve(XtX)
sdv <- sqrt(diag(Bbar))
cat('unconstrained posterior mean:', round(bbar,3), '\n')
cat('posterior sd               :', round(sdv,3), '   (Python: mean 1.644/0.208/0.064/0.815)\n')

unconstrained posterior mean: 1.644 0.208 0.064 0.815 

posterior sd               : 0.206 0.245 0.195 0.219    (Python: mean 1.644/0.208/0.064/0.815)
In [2]:
# p_{2|1} = P(beta1>0, beta2>0) via TruncatedNormal (analytic-grade), the GHK analogue
mu12 <- bbar[2:3]; S12 <- Bbar[2:3,2:3]
p21 <- TruncatedNormal::pmvnorm(mu=as.numeric(mu12), sigma=S12, lb=c(0,0), ub=c(Inf,Inf))
cat(sprintf('p_{2|1} = P(b1>0, b2>0) = %.4f\n', as.numeric(p21)))

# constrained posterior draws via tmvtnorm (truncated-MVN Gibbs) -- the Geweke constrained posterior
set.seed(1)
lower <- c(-Inf, 0, 0, -Inf); upper <- rep(Inf, 4)
draws <- tmvtnorm::rtmvnorm(n=20000, mean=as.numeric(bbar), sigma=Bbar, lower=lower, upper=upper, algorithm='gibbs',
                  burn.in.samples=2000)
cat('constrained posterior mean:', round(colMeans(draws),3), '\n')
cat('constrained posterior sd  :', round(apply(draws,2,sd),3), '\n')
cat('all draws satisfy b1>0,b2>0:', all(draws[,2]>0 & draws[,3]>0), '\n')

p_{2|1} = P(b1>0, b2>0) = 0.5234

constrained posterior mean: 1.63 0.307 0.186 0.805 

constrained posterior sd  : 0.205 0.191 0.131 0.22 

all draws satisfy b1>0,b2>0: TRUE
In [3]:
# bain: Bayes factor for the inequality hypothesis b1>0 & b2>0
fit <- lm(y ~ x1 + x2 + x3, data=d)
bf  <- bain(fit, 'x1 > 0 & x2 > 0')
print(bf)

Bayesian informative hypothesis testing for an object of class lm (continuous predictors):

   Fit   Com   BF.u  BF.c  PMPa  PMPb  PMPc 
H1 0.479 0.279 1.715 2.373 1.000 0.632 0.704
Hu                               0.368      
Hc 0.521 0.721 0.723                   0.296

Hypotheses:
  H1: x1>0&x2>0

Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement. PMPa contains the posterior model probabilities of the hypotheses specified. PMPb adds Hu, the unconstrained hypothesis. PMPc adds Hc, the complement of the union of the hypotheses specified.
In [4]:
options(repr.plot.width=7, repr.plot.height=5.5)
unc <- mvtnorm::rmvnorm(8000, mean=as.numeric(bbar), sigma=Bbar)      # unconstrained draws
df  <- rbind(data.frame(b1=unc[,2], b2=unc[,3], type='unconstrained'),
             data.frame(b1=draws[sample(nrow(draws),8000),2], b2=draws[sample(nrow(draws),8000),3], type='constrained'))
ggplot(df, aes(b1, b2, colour=type)) +
  annotate('rect', xmin=0, xmax=Inf, ymin=0, ymax=Inf, fill='green', alpha=.07) +
  geom_point(size=.5, alpha=.15) +
  geom_hline(yintercept=0, linetype='dotted', colour='red') + geom_vline(xintercept=0, linetype='dotted', colour='red') +
  scale_colour_manual(values=c('unconstrained'='steelblue','constrained'='orange')) +
  coord_cartesian(xlim=quantile(unc[,2],c(.005,.995)), ylim=quantile(unc[,3],c(.005,.995))) +
  labs(title='Feasible region: posterior clipped to b1>0, b2>0 (tmvtnorm)', x=expression(beta[1]), y=expression(beta[2]), colour='') +
  theme_minimal(base_size=13) + theme(legend.position='top') + guides(colour=guide_legend(override.aes=list(alpha=1,size=2)))
No description has been provided for this image

Results¶

  • TruncatedNormal::pmvnorm gives the constraint probability $p_{2|1}$ matching the Python GHK estimate — to analytic accuracy via minimax tilting.
  • tmvtnorm::rtmvnorm reproduces the constrained posterior (means/SDs agree with the Python Gibbs sampler), and every draw respects $\beta_1>0,\ \beta_2>0$ — the truncated-MVN Gibbs is exactly Geweke's constrained posterior.
  • bain turns the constraint into a Bayes factor — the evidence for $\{\beta_1>0,\beta_2>0\}$ relative to its complement / the unconstrained model — the model-selection use of $p_{2|1}$.
  • The feasible-region plot mirrors the Python geometry figure: the unconstrained cloud is sliced to the positive quadrant. Three R packages cover the three Python algorithms on identical data.