BayesMixReg<-function(y,X,runs=5000,N=5,D=1){   
  #y:       vector of observations
  #X:       matrix of predictors
  #N:       Number of mixture components
  library(mvtnorm)
  library(MCMCpack)

  n<-length(y)
  p<-ncol(X)
  tXXinv<-solve(t(X)%*%X)

  #Generate initial values
  beta<-rep(0,p)
  theta<-rep(0,N)
  tau1<-tau2<-1
  probs<-rep(1/N,N)
  g<-sample(1:N,n,replace=T,prob=probs)

  #keep track of stuff
  keepbeta<-keepgamma<-matrix(0,runs,p)
  x.grid<-seq(-10,10,0.1)
  density<-matrix(0,runs,length(x.grid))

  #START THE MCMC!
  for(i in 1:runs){


     #update beta:
     VAR<-tXXinv/tau1
     MN<-tau1*VAR%*%t(X)%*%(y-theta[g])
     beta<-as.vector(rmvnorm(1,MN,VAR))
     resids<-y-X%*%beta
   
     #update the variances
     tau1<-rgamma(1,n/2+0.1,sum((resids-theta[g])^2)/2+0.1)
     tau2<-rgamma(1,N/2+0.1,sum(theta^2)/2+0.1)

     #update theta
     for(j in 1:N){
       VAR<-1/(tau1*sum(g==j)+tau2)
       MN<-VAR*sum(tau1*resids[g==j])
       theta[j]<-rnorm(1,MN,sqrt(VAR))
     }

     #update g
     for(s in 1:n){
       g[s]<-sample(1:N,1,prob=probs*dnorm(resids[s],theta,1/sqrt(tau1)))
     }

     #new v
     P<-rep(D,N)
     for(j in 1:N){P[j]<-P[j]+sum(g==j)}     
     probs<-rdirichlet(1,P)


     #keep track of stuff
     keepbeta[i,]<-beta
     for(j in 1:N){
        density[i,]<-density[i,]+probs[j]*dnorm(x.grid,theta[j],1/sqrt(tau1))
     }

     if(i%%100==0){
       par(mfrow=c(1,1))
       h<-hist(y-X%*%beta,breaks=25,main=i)
       lines(x.grid,max(h$count)*density[i,]/max(density[i,]),lwd=2,col=4)
     }
  }

list(beta=keepbeta,grid=x.grid,density=density)}



#A simulated example:

set.seed(0820)
n<-200
p<-5
x<-matrix(rnorm(n*p),n,p)
x[,1]<-1

g<-rep(1,n)
g<-rbinom(n,1,.8)
y<-x[,1]+2*x[,2]+g*rnorm(n)+(1-g)*rnorm(n,3,3)
print(summary(lm(y~x[,-1])))

fit<-BayesMixReg(y,x)

print("Posterior summary of the betas")
mn<-apply(fit$beta[1000:5000,],2,mean)
sd<-apply(fit$beta[1000:5000,],2,sd)
low<-apply(fit$beta[1000:5000,],2,quantile,0.025)
high<-apply(fit$beta[1000:5000,],2,quantile,0.975)
results<-cbind(mn,sd,low,high)
rownames(results)<-paste("x",0:4,sep="")
print(round(results,3))

#plot the estimated density:
d<-apply(fit$density[1000:5000,],2,mean)
plot(fit$grid,d,ylim=range(d),type="l",
ylab="post mean dense",xlab="y")