Custom Distribution - Multivariate SkewNormal

The random method of CustomDist uses python code so no need to convert it to PyTensor. You can just use what the blogpost does.

There are some examples at the end of the CustomDist documentation.

For the shape error, can you provide a small reproducible snippet that triggers it? The smaller the better