Why would scipy/numpy wrapped in @as_op cause much faster sampling than using pytensor operations?

The LKJCholskyCov implicitly assumes the sd_dist is a distribution on the positive real line. It applies a LogTransform that is only really valid/useful for these, when doing the relevant transformation / computing the jacobian (my details here are hazy, but @aseyboldt can probably clarify).

A TruncatedNormal would suffer the same fate (divergences if the sampler gets close to the boundaries) because that would require an IntervalTransform to make it properly unbounded, but the LKJ implies/hardcodes a LogTransform. The LogTransform may actually make things worse because of the backward exponentiation when evaluating the logp.

I did try to make this explicit when I wrote sd_dist: A positive scalar distribution..., but should have been explicit about being unbounded on the right.

Would be good to update the docs to reflect this. Would you be interested in doing it?