SSVS<-function(y,x,
               runs=20000,burn=5000,update=1000,
               a1=0.001,b1=0.001,prec.beta=0.1,inprob=0.5){

    ################            model        ###############:
    # y~N(int+x%*%beta,inverse variance = taue)
    # beta[j]<-delta[j]*alpha[j]
    # int~flat
    # taue~gamma(a1,b1)
    # alpha[j]~N(0,inv.var=prec.beta)
    # delta[j]~Bern(inprob)
    ########################################################:

    n<-length(y)
    p<-ncol(x)

    #initial values:
    int<-mean(y)
    beta<-alpha<-inout<-rep(0,p)
    taue<-1/var(y)

    #keep track of stuff:
    keep.beta<-matrix(0,runs,p)
    keep.int<-keep.taue<-rep(0,runs)

    txx<-t(x)%*%x
    txy<-t(x)%*%y
    tx1<-t(x)%*%(0*y+1)

    #LET'S ROLL:
    for(i in 1:runs){

      taue<-rgamma(1,n/2+a1,sum((y-int-x%*%beta)^2)/2+b1)
      int<-rnorm(1,mean(y-x%*%beta),1/sqrt(n*taue))
   
      #update alpha
      z<-x%*%diag(inout)
      V<-solve(taue*t(z)%*%z+prec.beta*diag(p))
      M<-taue*t(z)%*%(y-int)
      alpha<-V%*%M+t(chol(V))%*%rnorm(p)
      beta<-alpha*inout
 
      #update inclusion indicators: 
      r<-y-int-x%*%beta
      for(j in 1:p){
        r<-r+x[,j]*beta[j]
        log.p.in<-log(inprob)-0.5*taue*sum((r-x[,j]*alpha[j])^2)
        log.p.out<-log(1-inprob)-0.5*taue*sum(r^2)
        diff<-log.p.in-log.p.out
        if(diff>10){diff<-10}
        p.in<-exp(diff)/(1+exp(diff))
        inout[j]<-rbinom(1,1,p.in)
        beta[j]<-inout[j]*alpha[j]
        r<-r-x[,j]*beta[j]
      }

      keep.beta[i,]<-beta
      keep.int[i]<-int
      keep.taue[i]<-taue

      if(i%%update==0){
        plot(beta,main=paste("Iteration",i))
        abline(0,0)
      }
    }
   
list(beta=keep.beta[burn:runs,],
     int=keep.int[burn:runs],
     taue=keep.taue[burn:runs])}

###############        Output        ##################:
# beta,int,taue,lambda = posterior samples
#######################################################:






#A simple example:
library(mlbench)
data(BostonHousing)
?BostonHousing
dat<-BostonHousing

y<-dat$medv
x<-cbind(dat$crim,dat$zn,dat$indus,dat$chas,dat$nox,dat$rm,
         dat$age,dat$dis,dat$rad,dat$tax,dat$ptratio,dat$b,
         dat$lstat)
colnames(x)<-c("crim","zn","indus","chas","nox","rm","age",    
 "dis","rad","tax","ptratio","b","lstat")

for(j in 1:ncol(x)){
  x[,j]<-(x[,j]-mean(x[,j]))/sd(x[,j])
}


fit<-SSVS(y,x)

boxplot(fit$beta~col(fit$beta),main="Posterior dist of beta",outline=F)
hist(fit$beta[,1],breaks=100,main="Posterior of beta1",xlab="beta1")

print("inclusion prob")
print(round(colMeans(fit$beta!=0),2))