Use PDMat variates for LKJ, Wishart and InverseWishart

[devmotion]

We have LKJCholesky as a distribution of Cholesky factors such that the corresponding unfactorized matrix follows the LKJ distribution of the same parameters. This is convenient since (copied from the docstring)

Sampling `LKJCholesky` is faster than sampling `LKJ`, and often having the correlation
matrix in factorized form makes subsequent computations cheaper as well.

The same would be true for e.g. Wishart and InverseWishart (ref https://arxiv.org/abs/2310.15884).

However, one drawback of LKJCholesky is that many statistics of interest such as mean etc. could only be defined if it is identified with the corresponding LKJ distribution: #1938 (review)

If you're actually interested in the LKJ (or Wishart or InverseWishart) distribution AND would like variates to be readily (and efficiently!) factorized, IMO a better approach would be to use PDMat variates, i.e., to directly sample the Cholesky factors (as in LKJCholesky) but return the unfactorized matrix (with readily available factorization) as a PDMat.

[sethaxen]

I like this plan. Since changing the variate-type of LKJ would be breaking, a less satisfactory but still non-breaking way would be to 1) add PDMatVariate 2) add LKJPDMat, and 3) deprecate LKJCholesky. In the next breaking release, LKJPDMat could be renamed to LKJ, and LKJCholesky could be removed.

[sethaxen]

This proposal would supersede #1787