ENH: stats: add multivariate_laplace

[jiawei-zhang-a]

Is your feature request related to a problem? Please describe.

Proposed new feature or change:

I am interested in submitting a new feature proposal as an issue. Specifically, I have noticed that scipy.stats.laplace only allows for the specification of location (mean) and scale (decay) parameters, and does not have an option for Multivariate Laplace. Conversely, while scipy.stats.multivariate_normal is available for normal distribution, there is currently no equivalent for multivariate Laplace.

mean = [0, 0]
cov = [[1, 0], [0, 100]]
x, y = scipy.stats.multivariate_normal(mean, cov).T
loc, scale = 0., 1.
s = scipy.stats.laplace(loc, scale, 1000)

Describe the solution you'd like.

I want to implement a new function multivariate_laplace
scipy.stats.multivariate_laplace
A multivariate laplace random variable.

The mean keyword specifies the mean. The cov keyword specifies the covariance matrix.

Parameters:
mean: array_like
Mean of the distribution.

cov: array_like
Symmetric positive (semi)definite covariance matrix of the distribution.

allow_singularbool, default: False
Whether to allow a singular covariance matrix. This is ignored if cov is a [Covariance]

seed{None, int, np.random.RandomState, np.random.Generator}, optional
Used for drawing random variates. If seed is None, the RandomState singleton is used. If seed is an int, a new RandomState instance is used, seeded with seed. If seed is already a RandomState or Generator instance, then that object is used. Default is None.

Describe alternatives you've considered.

No response

Additional context (e.g. screenshots, GIFs)

No response

[jiawei-zhang-a]

Thanks! @rkern @j-bowhay

[mansanlab]

Here is a function that can perform your operation while you wait for the scipy fix.
multivariate-Laplace-extension-for-SciPy

[mdhaber]

@jiawei-zhang-a Please suggest this on the forum if you are still interested in this feature.

[mdhaber]

Which version of the multivariate Laplace distribution did you want? There are a few.

Here is some code for the version on Wikipedia, adapted from the multivariate normal:

import numpy as np
import matplotlib.pyplot as plt

from scipy.stats._multivariate import multi_rv_generic, mvn_docdict_params, _LOG_2PI
from scipy import stats, special, integrate
from scipy._lib import doccer
rng = np.random.default_rng(3409235845834753)

class multivariate_laplace_gen(multi_rv_generic):

    def __init__(self, seed=None):
        super().__init__(seed)
        self.__doc__ = doccer.docformat(self.__doc__, mvn_docdict_params)

    def _process_parameters(self, mean, cov):
        """
        Infer dimensionality from mean or covariance matrix, ensure that
        mean and covariance are full vector resp. matrix.
        """
        if not isinstance(cov, stats.Covariance):
            L = np.linalg.cholesky(cov)
            cov = stats.Covariance.from_cholesky(L)

        dim = cov.shape[-1]
        mean = np.array([0.]) if mean is None else mean
        message = (f"`cov` represents a covariance matrix in {dim} dimensions,"
                   f"and so `mean` must be broadcastable to shape {(dim,)}")
        try:
            mean = np.broadcast_to(mean, dim)
        except ValueError as e:
            raise ValueError(message) from e
        return dim, mean, cov

    def logpdf(self, x, mean=None, cov=1):
        dim, mean, cov = self._process_parameters(mean, cov)
        x = np.asarray(x)
        return self._logpdf(x, mean, cov)

    def pdf(self, x, mean=None, cov=1):
        return np.exp(self.logpdf(x, mean, cov))

    def _logpdf(self, x, mean, cov):
        log_det_cov, k = cov.log_pdet, cov.rank
        dev = x - mean
        if dev.ndim > 1:
            log_det_cov = log_det_cov[..., np.newaxis]
            k = k[..., np.newaxis]
        maha = np.sum(np.square(cov.whiten(dev)), axis=-1)
        v = 1 - k/2
        return (np.log(2*special.kv(v, np.sqrt(2*maha)))
                - (k*_LOG_2PI + log_det_cov - v * np.log(maha/2))/2)

    def rvs(self, mean=None, cov=1, size=1, random_state=None):
        # https://www.sciencedirect.com/science/article/pii/S0047259X12000516
        # Something doesn't seem right - in the plot, the random variates 
        # seem too disperse compared to the PDF.
        dim, mean, cov = self._process_parameters(mean, cov)
        random_state = self._get_random_state(random_state)
        # `gamma(1, ...)` same as `standard_exponential`
        Z = random_state.gamma(1, size=size)[..., np.newaxis]
        X = stats.multivariate_normal.rvs(mean=0, cov=cov, size=size, 
                                          random_state=random_state)
        return mean*Z + np.sqrt(Z)*X

multivariate_laplace = multivariate_laplace_gen()

def pdf(x, mean, cov):
    # reference PDF
    k = cov.shape[-1]
    v = (2 - k) / 2
    x = x - mean
    maha = x @ (np.linalg.inv(cov) @ x)
    det = np.linalg.det(cov)
    den = (2*np.pi)**(k/2) * det**0.5
    t1 = (maha/2)**(v/2)
    t2 = special.kv(v, np.sqrt(2*maha))
    return 2 / den * t1 * t2

n = 2
x = rng.random(n)
cov = rng.random((n, n))
cov = cov + cov.T + n*np.eye(n)
mean = rng.random(n)

# check PDF against simple reference implementation
multivariate_laplace.pdf(x, mean, cov), pdf(x, mean, cov)

# Ensure that PDF integrates to 1
# integrate.dblquad(lambda x, y: pdf([x, y], mean, cov), -np.inf, np.inf, -np.inf, np.inf)

# compare random variates against PDF
# seems too disperse
x = multivariate_laplace.rvs(mean, cov, size=100)
plt.plot(x[:, 0], x[:, 1], '.')
x = np.linspace(-10, 10)
y = np.linspace(-10, 10)
x, y = np.meshgrid(x, y, indexing='ij')
xy = np.moveaxis(np.stack((x, y)), 0, -1)
pdf = multivariate_laplace.pdf(xy, mean, cov)
plt.contour(x, y, pdf)

@jiawei-zhang-a Please suggest this on the forum if you are still interested in this feature. Otherwise, I think we can close this.