gbd <-
function(n=1000,p=1,intercept=TRUE,betas=rep(0,p+intercept),mu=rep(0,p),Sigma=diag(1,p),seed=1){
	require(MASS)
	set.seed(seed)
	X<-mvrnorm(n,mu,Sigma)
	Xm<-X
	if(intercept){Xm<-cbind(rep(1,n),X)}
	lp<-Xm %*% betas
	pv<-expit(lp)
	Y<-rbinom(n,1,pv)
	return(list(Y=Y,X=X))
}