Feature: Implementation Plan for Multivariate BART Support

[miaoqingyu2]

Following up on the general request in issue #2 and my conversation with Giacomo, I'm posting a more detailed implementation plan here to help coordinate the effort for adding multivariate response support to bartz.

Here’s the basic plan I have in mind:

0. Core Structure:

• The model will support a y with shape (n_samples, d).
• All d responses will share the same underlying forest structure.

1. Priors:

• Switch the leaf node priors from a univariate Normal to a Multivariate Normal distribution.
• Swap out the Inverse-Gamma prior for sigma2 with an Inverse-Wishart prior to handle the error covariance matrix. It would also be great to add support for more flexible priors later if possible.

2. Outputs:

• The predict() method will return multivariate predictions, for instance, an array of shape (n_mcmc, n_samples, d).
• The model will return the posterior trace of the full error covariance matrix, which will allow for calculating the correlation matrix between the response variables.

I'm open to any feedback or suggestions on this plan. Looking forward to coordinating with everyone interested.