BayesSemiPar<-function(y,x,
               iters=10000,burn=1000,update=1000,
               a1=0.01,b1=0.01,a2=0.01,b2=0.01,nu=1){

    ################            model        ###############:
    # y~N(int+x%*%beta,inverse variance = taue)
    # beta[j]~t_nu(0,1/lambda) <=>
    # beta[j]~N(0,taub[j]) & taub[j]~G(nu/2,nu*lambda/2) 
    # int~flat
    # taue~gamma(a1,b1)
    # lambda~gamma(a2,b2)
    ########################################################:

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

    #initial values:
    init<-lm(y~x)$coef
    int<-init[1]
    beta<-init[-1]
    taub<-rep(100,p)
    lambda<-1
    taue<-10

    #keep track of stuff:
    keep.beta<-matrix(0,iters,p)
    txx<-t(x)%*%x
 
    F1<-F2<-rep(0,n)
    #START THE MCMC:
    for(i in 1:iters){

      taue<-rgamma(1,n/2+a1,sum((y-int-x%*%beta)^2)/2+b1)

      taub<-rgamma(p,1/2+nu/2,beta*beta/2+nu*lambda/2)

      lambda<-rgamma(1,p*nu/2+a2,sum(taub)*nu/2+b2)

      int<-rnorm(1,mean(y-x%*%beta),1/sqrt(n*taue))
      
      #update beta
      VVV<-solve(taue*txx+diag(taub))
      MMM<-taue*t(x)%*%(y-int)
      beta<-VVV%*%MMM+t(chol(VVV))%*%rnorm(p)
      beta<-as.vector(beta)

      keep.beta[i,]<-beta

      F<-int+x%*%beta
      if(i%%update==0){
        plot(y,main=i)
        lines(F)
      }

      if(i>burn){
         F1<-F1+F/(iters-burn)
         F2<-F2+F*F/(iters-burn)
      }
    }
    

list(beta=keep.beta[burn:iters,],
     fitted.mn=F1,
     fitted.var=F2-F1^2)}

###############        Output        ##################:
# beta = samples for beta
# fitted = posterior mean of int + x%*%beta
#######################################################:


#Load the motorcylce data
library(splines)
library(MASS)
y<-mcycle$accel
t<-mcycle$times
y<-(y-mean(y))/sd(y)
t<-t/max(t)


#Plot the basis functions:
B<-ns(t,10)
plot(t,t,col=0,ylim=range(B),xlab="t",ylab="Bj(t)")
for(j in 1:ncol(B)){lines(t,B[,j])}


#Fit with 10 or 50 knots and Gaussian (nu=100) or Cauchy (nu=1) coefficients:
fitG10<-BayesSemiPar(y,ns(t,10),nu=100)
fitT10<-BayesSemiPar(y,ns(t,10),nu=1)
fitG50<-BayesSemiPar(y,ns(t,50),nu=100)
fitT50<-BayesSemiPar(y,ns(t,50),nu=1)


#plot the fitted values:
plot(t,y)
lines(t,fitG10$fitted.mn,lwd=1,lty=2)
lines(t,fitT10$fitted.mn,lwd=2,lty=2)
lines(t,fitG50$fitted.mn,lwd=1,lty=1)
lines(t,fitT50$fitted.mn,lwd=2,lty=1)
legend("topleft",c("J=10, Gauss","J=10, Cauchy","J=50, Gauss","J=50, Cauchy"),
  lwd=c(1,2,1,2),lty=c(2,2,1,1),inset=0.05)

X11()
par(mfrow=c(2,2),mar=c(2,2,2,2))
boxplot(fitG10$beta,outline=F,ylim=c(-6,3),main="Gaussian, J=10")
boxplot(fitT10$beta,outline=F,ylim=c(-6,3),main="Cauchy, J=10")
boxplot(fitG50$beta,outline=F,ylim=c(-6,3),main="Gaussian, J=50")
boxplot(fitT50$beta,outline=F,ylim=c(-6,3),main="Cauchy, J=50")