dec2bin = function(x, b=2){
  N <- length(x)
  xMax <- max(x)
  ndigits <- (floor(logb(xMax, base=2))+1)
  Base.b <- array(NA, dim=c(N, ndigits))
  for(i in 1:ndigits)
    {
    Base.b[, ndigits-i+1] <- (x %% b)
    x <- (x %/% b)
    }
  if(N ==1) Base.b[1, ] else Base.b
  }

model.eigen = function(delta.vec, X.mat)
{
X = X.mat[,delta.vec==1]
XtX = t(X)%*%X
eigen(XtX)
}

marginal.prob.weight = function(weight, lambda, model, X, y, vo, model.eigens)
{
n = length(y)
vec1 = apply(model,1,paste,collapse="")
prior.prob=numeric(length(vec1))
y.given.delta=numeric(length(vec1))
for(i in seq(along=vec1)){
 	delta.vec = model[i,]
	absgamma = sum(delta.vec) 
	TempX = X[,delta.vec==1]
      D.model = model.eigens[[i]]$values
      P.model = model.eigens[[i]]$vectors
      n.pars = length(D.model)
	g.lambda = sum(D.model^(lambda))/(n*sum(D.model^(-1)))
	prior.prob[i] = weight^(absgamma)*(1-weight)^(p-absgamma)*sqrt(prod(D.model^(-lambda)))
	Sbeta = P.model%*%diag(D.model^lambda,nrow=n.pars,ncol=n.pars)%*%t(P.model)/g.lambda
	I.mat = diag(rep(1,n)) 
	term1 = I.mat + TempX%*%Sbeta%*%t(TempX)
	y.given.delta[i] = dmvt(as.vector(y),rep(0,n),term1,n+vo,log=FALSE)
	
} 

un.norm.post = y.given.delta*prior.prob
marginal.prob = log(sum(un.norm.post)/sum(prior.prob))
marginal.prob
}

marginal.prob.lambda = function(lambda, model, X, y, vo, model.eigens)
{
eps = .0000001
optimize(marginal.prob.weight, c(eps,1-eps), maximum=TRUE, lambda=lambda, model=model, X=X, y=y, vo=vo, model.eigens=model.eigens)$objective
}

fixed.lambda.opt = function(X, y, lambda, vo=.01)
{
require(mvtnorm)
n=nrow(X)
p=ncol(X)
X=scale(X)
X=scale(X,scale=rep(sqrt(n-1),p))[1:n,] 
y=y-mean(y)

model=dec2bin(1:(2^p-1))
model.eigens = apply(model,1,model.eigen, X.mat=X)

eps = .0000001
opt.weight.set = optimize(marginal.prob.weight, c(eps,1-eps), maximum=TRUE, lambda=lambda, model=model, X=X, y=y, vo=vo, model.eigens=model.eigens)
opt.set = NULL
opt.set$opt.weight = opt.weight.set$maximum
opt.set$opt.value = opt.weight.set$objective
return(opt.set)
}

full.marginal.opt = function(X, y, lambda.bounds=c(-3,3), vo=.01)
{
require(mvtnorm)
n=nrow(X)
p=ncol(X)
X=scale(X)
X=scale(X,scale=rep(sqrt(n-1),p))[1:n,] 
y=y-mean(y)

model=dec2bin(1:(2^p-1))
model.eigens = apply(model,1,model.eigen, X.mat=X)

opt.lambda.set = optimize(marginal.prob.lambda, c(-3,3), maximum=TRUE, model=model, X=X, y=y, vo=vo, model.eigens=model.eigens)
opt.set = NULL
opt.set$opt.lambda = opt.lambda.set$maximum
opt.set$opt.weight = fixed.lambda.opt(X, y, opt.set$opt.lambda, vo)$opt.weight
opt.set$opt.value = opt.lambda.set$objective

return(opt.set)
}

post.probs = function(X, y, lambda, weight, vo=.01)
{
require(mvtnorm)
n=nrow(X)
p=ncol(X)
X=scale(X)
X=scale(X,scale=rep(sqrt(n-1),p))[1:n,] 
y=y-mean(y)

model=dec2bin(1:(2^p-1))
model.eigens = apply(model,1,model.eigen, X.mat=X)

vec1 = apply(model,1,paste,collapse="")
prior.prob=numeric(length(vec1))
y.given.delta=numeric(length(vec1))
for(i in seq(along=vec1)){
 	delta.vec = model[i,]
	absgamma = sum(delta.vec) 
	TempX = X[,delta.vec==1]
      D.model = model.eigens[[i]]$values
      P.model = model.eigens[[i]]$vectors
      n.pars = length(D.model)
	g.lambda = sum(D.model^(lambda))/(n*sum(D.model^(-1)))
	prior.prob[i] = weight^(absgamma)*(1-weight)^(p-absgamma)*sqrt(prod(D.model^(-lambda)))
	Sbeta = P.model%*%diag(D.model^lambda,nrow=n.pars,ncol=n.pars)%*%t(P.model)/g.lambda
	I.mat = diag(rep(1,n)) 
	term1 = I.mat + TempX%*%Sbeta%*%t(TempX)
	y.given.delta[i] = dmvt(as.vector(y),rep(0,n),term1,n+vo,log=FALSE)
	
} 
post.summary = NULL
un.norm.post = y.given.delta*prior.prob
post.summary$post.prob = un.norm.post/sum(un.norm.post)

post.summary$models = vec1
post.summary$lambda =lambda
post.summary$weight = weight
post.summary$log.marginal = log(sum(un.norm.post)/sum(prior.prob))

m = length(post.summary$models)
ord=order(post.summary$post.prob, decreasing=TRUE)
prmatrix(cbind( post.summary$models[ord], round(post.summary$post.prob[ord],6)), rowlab=rep("",m), collab=c( "model", "posterior"),quote=F)
cat("----------------------------------------------------------
------------------------------------------------------------------",fill=T)

return(post.summary)
}