PS: I actually just tried the example above (well I had to fixed a few typos), with the lastest 5.12.0 pymc version, and it does not work (shape mismatch). Looking at the latest doc, and the shapes in the error message, it is clear that it cannot work because pmx.distributions.histogram_approximation operates on the first dimension only → it returns univariate histograms for each obs. It does not account for correlations.
The issue is probably best described with copulas. However, the linked notebook does not seem to address the issue directly, but works to transform the data to a normal distribution first. It does not look like a general solution, but I may be mistaken.