# Custom Distribution - Multivariate SkewNormal

**URL:** <https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490>\
**Category:** General\
**Created:** [December 19, 2023, 7:12pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490 "2023-12-19T19:12:01Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 19, 2023, 7:12pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/1 "2023-12-19T19:12:01Z")

</div>

I am trying to implement a multivariate skewnormal distribution as a custom distribution in PyMC. I am going based on this implemenation of it: [A Python Implementation of the Multivariate Skew Normal](https://gregorygundersen.com/blog/2020/12/29/multivariate-skew-normal/)

My question is what is the best way to implement this? I have tried to use it as a custom distribution but have run into numerous errors. Would it be better to extend the MvNormal class or use CustomDist?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 20, 2023, 11:35am UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/2 "2023-12-20T11:35:12Z")

</div>

CustomDist is the best route

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 20, 2023, 4:15pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/3 "2023-12-20T16:15:36Z")

</div>

Does this implementation make sense?

```auto

class multivariate_skewnorm(multi_rv_generic):

    def __init__ (self, mu, shape, cov=None):
        self.dim = len(shape)
        self.shape = np.asarray(shape)
        self.mean = np.asarray(mu)
        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(value, self.dim)
        pdf = mvn(self.mean, self.cov).logpdf(x)
        cdf = norm(0, 1).logcdf(np.dot(x, self.shape))
        return _squeeze_output(np.log(2) + pdf + cdf)

    def rvs(self, size=1, random_state=43):
        random_state = self._get_random_state(random_state)
        aCa = self.shape @ self.cov @ self.shape
        delta = (1 / np.sqrt(1 + aCa)) * self.cov @ self.shape
        cov_star = np.block([[np.ones(1), delta],
                            [delta[:, None], self.cov]])
        if size == None:
            size = 1
        else:
            size = size[0]
        x = random_state.multivariate_normal(np.zeros(self.dim+1), cov_star, size=size)
        x0, x1 = x[:, 0], x[:, 1:]
        inds = x0 <= 0
        x1[inds] = -1 * x1[inds]
        return x1

@as_op(itypes=[pt.dmatrix, pt.dmatrix, pt.dmatrix, pt.dvector], otypes=[pt.dscalar])
def logp(value, mu, cov, alpha):
    return multivariate_skewnorm(mu, alpha, cov).logpdf(value)

def random(mu, cov, alpha, rng=None, size=None):
    return multivariate_skewnorm(mu, alpha, cov).rvs(size=size, random_state=23)

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 20, 2023, 6:07pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/4 "2023-12-20T18:07:54Z")

</div>

~That looks like the old pymc3 API, I suggest you use PyMC. CustomDist should require much less code: [https://www.pymc.io/projects/docs/en/latest/api/distributions/generated/pymc.CustomDist.html~](https://www.pymc.io/projects/docs/en/latest/api/distributions/generated/pymc.CustomDist.html~)

Nevermind, I got confused. You’re just showing the scipy implementation.

One thing to note is that you won’t be able to use NUTS because you’re using a custom Op. It would be better if you could define the logp using PyTensor operations, in which case you would get gradients for free, but those percentiles may be tricky.

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 20, 2023, 6:22pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/5 "2023-12-20T18:22:22Z")

</div>

Gotcha.

I am probably in over my head if I don’t fully understand the math behind the MV SkewNormal then?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 20, 2023, 8:53pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/6 "2023-12-20T20:53:54Z")

</div>

Yeah could be tricky, special because it looks like it may not have a closed form solution so you need to do some approximation

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 20, 2023, 9:04pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/7 "2023-12-20T21:04:47Z")

</div>

Any resources to look at if I wanted to take a stab at it?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 20, 2023, 9:23pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/8 "2023-12-20T21:23:10Z")

</div>

Actually the function \_process\_quantiles is not doing anything fancy so you should be able to implement a `logp` function for CustomDist without too much trouble.

> <https://github.com/scipy/scipy/blob/b3ba73a1526441824935e0c3d424415892bec84a/scipy/stats/_multivariate.py#L1185>

You can get the MvNormal logp and Univariate Normal logcdf from pymc via `pm.logp(pm.MvNormal.dist(mu, cov), x)` and `pm.logcdf(pm.Normal.dist(0, 1), x)`

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 20, 2023, 10:19pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/9 "2023-12-20T22:19:05Z")

</div>

I believe this is the correct logp:

```auto
def logp(value, mu, cov, alpha):
    return pm.math.log(2) + pm.logp(pm.MvNormal.dist(mu, cov), value) + pm.logcdf(pm.Normal.dist(0, 1), pm.math.dot(value, alpha))

```

What would the random function be? Would I just convert the rvs function above to use pm.MvNormal.dist?

edit: these are my attempts at logp and random. I am running into a shape mismatch error somewhere:

```auto
def logp(value, mu, cov, alpha):
    mv_normal = pm.MvNormal.dist(mu=mu, cov=cov)
    
    # Calculate the log probability of the value under the multivariate normal distribution
    log_prob_mv_normal = pm.logp(mv_normal, value - mu)
    
    # Calculate omega inverse
    omega_inv = pt.nlinalg.matrix_inverse(cov)
    
    # Calculate the argument for the standard normal CDF
    arg_cdf = pm.math.dot(pm.math.dot(alpha.T, omega_inv), (value - mu))
    
    # Instantiate the standard normal distribution for the CDF
    std_normal = pm.Normal.dist(mu=0, sigma=1)
    
    # Calculate the log CDF of the argument under the standard normal distribution
    log_cdf_std_normal = pm.logcdf(std_normal, arg_cdf)
    
    # Return the log of the skew-normal density
    return pm.math.log(2) + log_prob_mv_normal + log_cdf_std_normal

```

```auto
def random(mu, cov, alpha, rng=None, size=None):
    aCa = alpha @ cov @ alpha
    delta = (1/ pm.math.sqrt(1 + aCa)) * cov @ alpha
    a = pt.concatenate((pt.ones(1), delta)).reshape((1, -1))
    b = pt.concatenate((delta[:, None], cov), axis=1)
    cov_star = pt.concatenate((a, b))
    x = pm.MvNormal.dist(pt.zeros(alpha.shape[0]+1), cov=cov_star, shape=(alpha.shape[0]+1, -1))
    x0 = x[:, 0]
    x1 = x[:, 1:]
    inds = x0 <= 0
    new_x1 = set_subtensor(x1[inds], -1*x1[inds])
    new_x1 = new_x1 + mu
    return new_x1.eval()

```

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 21, 2023, 12:34am UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/10 "2023-12-21T00:34:18Z")

</div>

Should the logp function return the log prob of each row or for the entire dataset? When using it in the CustomDist, it passes the entire dataset into the value argument but I am unsure if it handles the summation on the backend or if I should be summing all the log probs.

---

<div class="post-metadata">

**Author:** ![cluhmann](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/cluhmann/32/3083_2.png) [@cluhmann](https://discourse.pymc.io/u/cluhmann)\
**Post date:** [December 21, 2023, 1:01am UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/11 "2023-12-21T01:01:10Z")

</div>

I believe the summation is handled my PyMC.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 21, 2023, 6:55am UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/12 "2023-12-21T06:55:51Z")

</div>

> [@Blake\_Walters](#):
>
> What would the random function be?

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

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 21, 2023, 4:12pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/13 "2023-12-21T16:12:23Z")

</div>

I figured out the shape error. However I am now getting this warning and the sampling is very slow:

```auto
/usr/local/lib/python3.10/dist-packages/pytensor/tensor/slinalg.py:391: LinAlgWarning: Ill-conditioned matrix (rcond=1.5831e-19): result may not be accurate.

```

I tried running with numpyro on a GPU in Colab but it just ran through the sampling without accepting any steps. I am guessing there is an issue with the compatibility of custom distributions and numpyro?

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [December 21, 2023, 4:37pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/14 "2023-12-21T16:37:39Z")

</div>

The error is saying that `cov_star` is not full rank, so the matrix inverse is breaking. You could try adding some small numbers to the diagonal to stabilize it?

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 21, 2023, 4:52pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/15 "2023-12-21T16:52:19Z")

</div>

How do I go about choosing what to add? I tested it with 0.1 but it is still showing the same warning.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 21, 2023, 5:05pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/16 "2023-12-21T17:05:07Z")

</div>

To debug, have you checked if unit covariance (without having to invert) is working correctly already?

---

<div class="post-metadata">

**Author:** ![Blake\_Walters](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/blake_walters/32/3546_2.png) [@Blake\_Walters](https://discourse.pymc.io/u/Blake_Walters)\
**Post date:** [December 21, 2023, 5:42pm UTC](https://discourse.pymc.io/t/custom-distribution-multivariate-skewnormal/13490/17 "2023-12-21T17:42:39Z")

</div>

The warning does not show up there. It still samples slowly but I should note that it is running on a 24x24 covariance matrix so maybe this is to be expected? The dataset is only ~4000 rows however.

I tested it on mock data generated from my random function and it still gets the same warning and runs at the same speed.

Here is the mock model and data as well as my updated logp and random functions if anyone wants to test it.

> <https://gist.github.com/bwalters13/60bc2fce973bb7f453899dd15ad7c73e>