fit data to a multivariate skew-normal distribution

Hi I would like to use the optimization tool bo to fit data to a multivriate skew-normal pdf. z = skewnormal(x,y) basica a 2D distribution fit. I found a explanation of makin the function in python (http://gregorygundersen.com/blog/2020/12/29/multivariate-skew-normal/) i would like to generate this as an anonymous function so that i can run a parametric fit.

Please help.