Chapter 9 Algorithm to fit HBART (stump)

---

Algorithm type: Metropolis within GIBBS for a hierachical BART model

Reason: We have closed posteriors for most parameters, but not for the tree structure

Data: Target variable $y$, groups $j = 1,\dots, J$

Result: Posterior distributions for all parameters

---

Initialisation;

Hyper-parameters values for $\alpha, \beta, k_2$;

Number of groups $J$;

Number of trees P;

Total number of nodes $N_p$

Number of observations $N =\sum_{j = 1}^{J} n_j$;

Number of iterations I;

• define $\mu_{\mu} = 0$, $\mu^{0}$, $\tau^{0}$, and $\mu_j^{0}, j = 1,\dots, J$, $k_1^0$, as the initial parameter values

• for i from 1 to I do:

  • for p from 1 to P do:

  • sample $\mu_p$ from the posterior $N(\frac{\mathbf{1}^{T} \Psi^{-1} R_p }{\mathbf{1}^{T} \Psi^{-1} \mathbf{1} + (k_2/P)^{-1}}, \tau^{-1} (\mathbf{1}^{T} \Psi^{-1} \mathbf{1} + (k_2/P)^{-1}))$

    • for j in 1:J do:
      • sample $\mu_{jp}$ from the posterior $MVN(\frac{P \mu /k_1 + \bar R_{pj} n_j}{(n_j + P/k_1)}, \tau^{-1} (n_j + P/k_1))$
    • end

  • set $R_{ijp} = Y_{ij} - \sum_{t = 1}^{p} \mu_{jt}$, or, in other words, the residuals for tree $p$ are the vector $y$ minus the sum of the $\mu_{jp}$ values up to tree p

  • end

• Define $\hat f_{ij}$ as the current overall prediction for observation $R_{ij}$, which will be $\sum_{p \neq pp}^{P} \mu_{jp}$, where pp represents the current tree index.

• sample $\tau$ from the posterior $Ga(\frac{N + J N_p + N_p}{2} + \alpha, \frac{\sum_{i= 1}^{N}(y_i - \hat f_i)^2}{2} + \frac{P \sum_{j, l, p}(\mu_{j, l, p} - \mu_{l, p})^2}{2 k_1} + \frac{P \sum_{l, p}\mu_{l, p}^2}{2 k_2} + \beta)$

• sample $k_1$ from a Uniform(0, u)

  • calculate p_old as the conditional distribution for the current $k_1$

  • calculate p_new as the conditional distribution for the newly sampled $k_1$

  • sample U from a Unif(0, 1)

  • if (a > u) accept the new $k_1$; else keep the current $k_1$

• end