Sample_ypred <- function(fit,NDraws=10){
  #extract the original data from the model (fixed and random effects)
  X <- fit$X
  M <- fit$M
  Samples <- list()
  #iterate for the number of draws required
  for (i in 1:NDraws){
    #samples the parameters from the posterior distributions of the fitted model
    CoeffFix <- apply(fit$beta.sample,2,function(x){return(sample(x,size=nrow(X)))})
    CoefRandom <- apply(fit$gamma.sample,2,function(x){return(sample(x,size=nrow(X)))})
    #Add all the valeus to get a prediction
    FixedPred <- CoeffFix*X
    RandomPred <- CoefRandom*M
    Pred <- rowSums(FixedPred)+rowSums(RandomPred)
    #Store it in a list
    Samples[[length(Samples)+1]]<-Pred
  }
  #Return the NDraws
  Preds <- do.call(cbind,Samples)
  return(Preds)
}