How can regression trees be fit in WinBUGS/OpenBUGS/JAGS? There is an R package called BayesTree which can fit regression trees in Bayesian environment. However, this way only simple regression is possible. I would like to use regression trees as a part of a bigger hierarchical model instead of the simple GLM formulas (CARTs are proven to give better results in certain applications, see e.g. Hu et al 2011).
Is it possible (and how) to fit regression trees as a part of the model in WinBUGS/OpenBUGS/JAGS? Are there any such packages for these pieces of software?
 A: Here is a suggestion of how the Bayesian CART model from Hu et al (2011) might be implemented in WinBUGS.  This method requires you to fix the tree depth in advance.  Given a complete binary tree, any node at depth 1 (i.e. children of the root) can be identified by an indicator $b_1$. Any node at depth 2 can be identified by a pair of indicators $(b_1,b_2)$.  And so on.  Let $x[i,v]$ be the $v$th predictor variable in the $i$th instance.  We can encode the path through the tree by thresholding $x[i,v]$ according to the parameters of the chosen node at each depth.  This code works for trees of depth 3:
for(i in 1:n) {
  b1[i] <- 1 + step(x[i,variable1] - split1)
  b2[i] <- 1 + step(x[i,variable2[b1[i]]] - split2[b1[i]])
  b3[i] <- 1 + step(x[i,variable3[b1[i],b2[i]]] - split3[b1[i],b2[i]])
  y[i] ~ dnorm(mean[b1[i],b2[i],b3[i]], prec[b1[i],b2[i],b3[i]])
}

You just need to add appropriate priors for $variable$ (the splitting variable at each node), $split$ (the splitting threshold at each node), and $(mean, prec)$ of the leaves.  There are 8 possible leaves, but they need not all be used.  For example, if one of the splits is outside the range of the data then it effectively prunes away part of the tree.  You could encourage such pruning with an appropriate prior on splits.
