creating a multivariate skew normal distribution python

Question

How can I create a multivariate skew normal function, where then by inputting x and y points we can create a surface diagram in 3d (x,y and z coordinates)

Answer 1

I wrote a blog post about this, but here is complete working code:

from   matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
from   scipy.stats import (multivariate_normal as mvn,
                           norm)


class multivariate_skewnorm:
    
    def __init__(self, a, cov=None):
        self.dim  = len(a)
        self.a    = np.asarray(a)
        self.mean = np.zeros(self.dim)
        self.cov  = np.eye(self.dim) if cov is None else np.asarray(cov)

    def pdf(self, x):
        return np.exp(self.logpdf(x))
        
    def logpdf(self, x):
        x    = mvn._process_quantiles(x, self.dim)
        pdf  = mvn(self.mean, self.cov).logpdf(x)
        cdf  = norm(0, 1).logcdf(np.dot(x, self.a))
        return np.log(2) + pdf + cdf


xx   = np.linspace(-2, 2, 100)
yy   = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(xx, yy)
pos  = np.dstack((X, Y))

fig  = plt.figure(figsize=(10, 10), dpi=150)
axes = [
    fig.add_subplot(1, 3, 1, projection='3d'),
    fig.add_subplot(1, 3, 2, projection='3d'),
    fig.add_subplot(1, 3, 3, projection='3d')
]

for a, ax in zip([[0, 0], [5, 1], [1, 5]], axes):
    Z = multivariate_skewnorm(a=a).pdf(pos)
    ax.plot_surface(X, Y, Z, cmap=cm.viridis)
    ax.set_title(r'$\alpha$ = %s, cov = $\mathbf{I}$' % str(a), fontsize=18)

That code will generate this figure: