## 2022-08-11 comment in response to an email question
## To replicate the table in the article, use the bsand_mcmc 
## function and set the priors using the following options:

## prior on B: 
### Jeffreys: jeff=TRUE, plugin=FALSE
### plugin: jeff=FALSE, plugin=TRUE

## prior on beta
### informative: beta0=beta,n0=n 
### uniform: beta0=rep(0,2),n0=Inf 

source("bsand_mcmc.r")

#####
n<-500
beta<-c(1,1)
p<-length(beta)

width_b<-width_m<-cover_b<-cover_m<-0
for(s in 1:10000)
{
  set.seed(s)
  X<-cbind(1,rexp(n))
  y<-rnorm(n,X%*%beta,X%*%beta)
  mle<-lm(y~ -1+X)$coef

  ### sandwich cov
  iA<-solve(t(X)%*%X)
  l.beta<-X*c(y-X%*%mle) ; B<-t(l.beta) %*% l.beta
  C.mle<- iA%*%B%*%iA
  Qm<-matrix(qnorm(c(.025,.975)),p,2,byrow=TRUE)
  cim<-mle[1:p]+diag(sqrt(diag(C.mle[1:p,1:p]))) %*% Qm
  cover_m<-cover_m+  ( cim[,1]< beta )  * ( cim[,2]>beta )
  width_m<-width_m + cim[,2]-cim[,1]
  ###

  ### full bayes
  BETA<-bsand_mcmc(y,X,n0=Inf,jeff=TRUE) 
#  BETA<-bmod_mcmc(y,X,n0=Inf,jeff=TRUE)
  cib<-t(apply(BETA,2,quantile,prob=c(.025,.975)))
  cover_b<-cover_b+  ( cib[,1]< beta )  * ( cib[,2]>beta )
  width_b<-width_b + cib[,2]-cib[,1]


  cat(n,s," ",round(cover_b/s,3)," " , round(width_b/s,3),"    ",
        round(cover_m/s,3)," " , round(width_m/s,3),"\n")
}