vsva_homo = function(y,X,backward = TRUE,a=1+1e-4,b=1+1e-4,tol=1e-10)
{
# Details:	Variable Selection with the Variational Approximation ranking
#		algorithm for the homoscedastic linear regression models.
#		The search can combine both forward and backward moves (if backward=T).
# 
# References: "Variational approximation for heteroscedastic models and matching 
#		pursuit algorithms" by David Nott et al.
# Author:	Minh-Ngoc Tran (ngoctm@nus.edu.sg)
# Inputs:
# y		vector of responses
# X		design matrix for the mean model (you need to include the intercept, 
#   		if desired).
# backward	whether or not to do backward selection besides forward selection
# a, b	hyperparameters in the inverse gamma priors for sigma2
# tol		stopping tolerance for variational lower bound

# Outputs:
# 
# C		logical vector indicating variables selected in the mean model
# mub		vector of estimated coefficients in the mean model
# muv		log of standard error
#
#======================================================

p = dim(X)[2]
n = dim(X)[1]

lmprior = function(C) {						# model prior 
lmprior = logfac(sum(C))+logfac(p-sum(C))
return(lmprior)
}
logfac <- function(x) {	if (x==0) return(0) else return(sum(log(1:x))) }

C = c(TRUE,rep(FALSE,p-1)) # initialize the active set for mean model
fit = vbhetlm_homo(y,as.matrix(X[,C]),a,b,tol)
Lopt = fit$lb[length(fit$lb)]
stoprule = FALSE
while (!stoprule) {
	mub = fit$mub
	muv = fit$muv
	sv = fit$sv
	sig2b = fit$sig2b
	XC = as.matrix(X[,C])

	L = rep(0,p)
	for (j in 1:p) {
		if (C[j]) L[j] = -Inf
		else {
			xj = X[,j]
			ld = exp(muv-1/2*sv)
			sig2bjhat = (1/sig2b+sum(xj*xj)/ld)
			sig2bjhat = 1/sig2bjhat
			mubjhat = sum(xj*(y-XC%*%mub))/ld
			mubjhat = mubjhat*sig2bjhat
			Ctry = C; Ctry[j] = TRUE
			L[j] = .5*log(sig2bjhat/sig2b)+.5*(1/sig2bjhat-1/sig2b)*mubjhat^2 #+ lmprior(Ctry)
		}
	}
	for (j in 1:p) if (L[j]==max(L)) {jhat = j; break }
	Ctry = C
	Ctry[jhat] = TRUE
	fittry = vbhetlm_homo(y,as.matrix(X[,Ctry]),a,b,tol)
	Ltry = fittry$lb[length(fittry$lb)]
	if (Ltry+lmprior(Ctry)>Lopt+lmprior(C)) { C = Ctry; Lopt = Ltry; fit = fittry} else {stoprule = TRUE}		
}
mub = rep(0,p); mub[C]=as.numeric(fit$mub)
muv = fit$muv

if (backward) {
	stoprule = FALSE
	if (sum(C)<=1) stoprule = TRUE
	while (!stoprule) {
		mub = as.matrix(rep(0,p))
		mub[C] = as.numeric(fit$mub)
		muv = fit$muv
		sv = fit$sv
		sig2b = fit$sig2b
		L = rep(0,p)
		for (j in 1:p) {
			if (!C[j]) L[j] = Inf
			else {
				Cnew = C; Cnew[j] = FALSE
				XC = as.matrix(X[,Cnew])
				xj = X[,j]
				ld = exp(muv-1/2*sv)
				sig2bjhat = (1/sig2b+sum(xj*xj)/ld)
				sig2bjhat = 1/sig2bjhat
				mub.new = mub[Cnew]
				mubjhat = sum(xj*(y-XC%*%mub.new))/ld
				mubjhat = mubjhat*sig2bjhat
				L[j] = .5*log(sig2bjhat/sig2b)+.5*(1/sig2bjhat-1/sig2b)*mubjhat^2	#-lmprior(Cnew)
			}
		}
		for (j in 1:p) if (L[j]==min(L)) {jhat = j; break }
		Ctry = C; Ctry[jhat] = FALSE
		fittry = vbhetlm_homo(y,as.matrix(X[,Ctry]),a,b,tol)
		Ltry = fittry$lb[length(fittry$lb)]
		if (Ltry+lmprior(Ctry)>Lopt+lmprior(C)) { C = Ctry; Lopt = Ltry; fit = fittry} else {stoprule = TRUE}
		if (sum(C)<=1) stoprule = TRUE
	}
}
list(C=C,mub=mub,muv=muv)
}

vbhetlm_homo<-function(y,X,a=1+1e-5,b=1+1e-5,tol=1e-10) {
# maximize the lower bound in the homoscedastic linear regression model
# with simplified priors
#
# Inputs:
#
# y		vector of responses
# X		design matrix for the mean model (you need to include the intercept, 
#   		if desired).
# a, b	hyperparameters of the inverse gamma priors 
# tol		stopping tolerance for variational lower bound
# 
# Outputs:
# 
# mub		mean in variational normal posterior for mean coefficients
# sb		covariance matrix in variational normal posterior for mean coefficients
# lb		vector of lower bound values over iterations
# ...

n<-length(y)
p<-dim(X)[2]

vlb<-function(y,X,p,n,pmub,pmuv,psb,psv,mub,muv,sb,sv) {
# a helper function to evaluate the lower bound
  t1<-solve(psb,sb)
  lb<-(p+1)/2-n/2*log(2*pi)+1/2*log(det(t1))+1/2*log(sv/psv)
  lb<-lb-1/2*sum(diag(t1))-1/2*sv/psv
  v1<-mub-pmub
  lb<-lb-1/2*sum(v1*solve(psb,v1))
  lb = lb - 1/2*(muv-pmuv)^2/psv-n/2*muv
  v = (y-X%*%mub)^2+apply((X%*%sb)*X,1,sum)
  v = sum(v)
  ld<-exp(-muv+1/2*sv)
  lb<-lb-1/2*v*ld
  return(lb)
}

lmfit = lm(y ~ X-1)
r = residuals(lmfit)
r<-log(r*r)
muv<-mean(r)
sv<-var(r)
mub<-lmfit$coef
sb<-vcov(lmfit)/summary(lmfit)$sigma^2

sig2b = (b+.5*sum(mub*mub)+.5*sum(diag(sb)))/(a+1+p/2)
psb = sig2b*diag(p)
psv = (b+.5*muv*muv+.5*sv)/(a+1+1/2)
pmub = rep(0,p)
pmuv = 0

lb<-vlb(y,X,p,n,pmub,pmuv,psb,psv,mub,muv,sb,sv)
lball<-lb
lbold<-lb-abs(lb)

while((lb-lbold)>tol) {

  lbold<-lb
  dd<-exp(-muv+1/2*sv)
  sb = solve(dd*t(X)%*%X+solve(psb))
  mub = sb%*%(solve(psb,pmub)+dd*t(X)%*%y)
  sig2b = (b+.5*sum(mub*mub)+.5*sum(diag(sb)))/(a+1+p/2)
  psb = sig2b*diag(p)

  v<-(y-X%*%mub)
  v<-v*v+apply((X%*%sb)*X,1,sum)
  v = sum(v)
  muv<-(v-n)/(v+2/psv)
  sv = v/2*exp(-muv)+1/psv
  sv = 1/sv
  psv = (b+.5*muv*muv+.5*sv)/(a+1+1/2)

  lb<-vlb(y,X,p,n,pmub,pmuv,psb,psv,mub,muv,sb,sv)
  lball<-c(lball,lb)
}

list(mub=mub,muv=muv,sb=sb,sv=sv,sig2b=sig2b,sig2v=psv,lb=lball)

}



