Implement Wishart(v, 𝗩) matrix-variate distribution

[jzeuzs]

This implements the Wishart Distribution or $\mathcal{W}_p (\nu, \mathbf{V})$, enabling matrix-variate statistical modeling.

Sampling is implemented via the Bartlett decomposition ($\mathbf{X} = \mathbf{L} \mathbf{A} \mathbf{A}^T \mathbf{L}^T$). This provides a numerically stable and efficient way to draw samples using the Cholesky factor of the scale matrix rather than using matrix multiplication.

Further, because the covariance of a Wishart distribution is a 4th-order tensor, this implementation flattens the result into a $p^2 \times p^2$ Matrix in order to allow us to stay within the existing Statistics trait's design constraints while providing full covariance/correlation information.

Results have been verified with scipy.stats.wishart.

[Axect]

Thanks for this, the Bartlett decomposition sampling and the scipy verification are much appreciated. The math (multivariate gamma, the vectorized p^2 x p^2 covariance) checks out on my side. A few points I'd like to sort out before merging.

On the O3 gating of pdf / ln_pdf.

• sample_with_rng genuinely needs O3 because cholesky has no pure-Rust path yet, so that gating is fine.
• But ln_pdf only uses det() and inv(), and both of those already have pure-Rust LU fallbacks (_ => self.lu().det() / self.lu().inv()). So the unimplemented!() on the non-O3 branch of ln_pdf is stricter than it needs to be.
• The side effect is that a user on the default feature set who calls pdf gets a runtime panic instead of a working result. Could you drop the O3 gate on ln_pdf and keep it only on the sampling path?

On input/support validation.

• Right now ln_pdf computes a finite value for inputs outside the support. For example, in even dimension x = -I has det = 1 > 0, so it returns a plausible-looking density even though the true value should be 0.
• MatDist::Wishart(df, scale) also does not check df > p - 1, so ChiSquared::new(df - i) can panic during sampling.
• For reference, the existing Dirichlet ln_pdf already asserts its support (x_i in (0,1), sum ~ 1), so matching that level of checking here would be consistent with the rest of the crate. Returning ln_pdf = -inf / pdf = 0 outside the support would be even nicer than a panic.

On tests.

The current tests focus on fixed 1x1 / 2x2 cases. It would be good to also cover:

• invalid df and non-SPD scale / x
• sample_with_rng reproducibility with a fixed seed
• empirical mean / covariance of many samples converging to mean() / cov()

Minor. In src/util/print.rs, impl Printable for Vec<Matrix> and impl Printable for &Vec<Matrix> are identical. You can forward one to the other to avoid the duplication.

Nice work overall. Once these are addressed I'm happy to merge. Thanks again.