Add generalized normal distribution

[shifdz]

Is it possible to use distributions from scipy stat or tensorflow probability in PYMC3? I need to use generalized normal distribution, but it is not avialable in PYMC3. Any help is really appreciated.

[ricardoV94]

You cannot use scipy/tensorflow directly, you need a function defined in theano.

Looking at the density funcition it should be pretty straightforward to define your own: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gennorm.html

I will mark this issue as a feature request. In the meantime if you want to go ahead and add your own distribution for your project I suggest you open an topic in our discourse: https://discourse.pymc.io/

Someone will probably help you out there

[shifdz]

I am new to programming and haven't worked in Theano, so I am not sure how to define a new distribution. Can someone please help me with implementing it? I opened a topic on Discourse, but didn't get any response.

[twiecki]

@shifdz There is a guide here: #4686 (comment)

[shifdz]

@twiecki Will the generalized normal distribution be added in pymc3 any time soon?

[twiecki]

@shifdz Whenever someone creates a pull request :) It's not that hard and we're happy to help if you want to give it a go.

[mvds314]

I recently made an implementation for the generalized normal. It's basically a copy of the Laplace and other distributions you already coded up. I don't really have time to go through the process of making a PR and such, but feel free to add it so that it can be properly maintained!

import numpy as np

import pymc as pm
import pytensor.tensor as pt
from pymc.distributions.dist_math import (
    check_icdf_parameters,
    check_icdf_value,
    check_parameters,
)
from pymc.distributions.distribution import Continuous
from pymc.distributions.shape_utils import rv_size_is_none
from pytensor.tensor.random.basic import ScipyRandomVariable

import scipy.special as sc
import scipy.stats as sps

class GenNormRV(ScipyRandomVariable):
        r"""A generalized normal continuous random variable.

        The probability density function for `gennorm` in terms of its shape
        parameter :math:`\beta`, location parameter :math:`\mu` and scale
        parameter :math:`\sigma` is:

        .. math::

            f(x; \beta, \mu, \sigma) = \frac{\beta}{2 \sigma} \exp\left(-\left(\frac{|x-\mu|}{\sigma}\right)^\beta\right)

        with :math:`\beta > 0`.
        """

        name: str = "gennorm"
        signature = "(),(),()->()"
        dtype: str = "floatX"
        _print_name: Tuple[str, str] = ("GenNorm", "\\operatorname{GenNorm}")

        def __call__(self, beta, loc=0.0, scale=1.0, size=None, **kwargs):
            r"""Draw samples from a generalized normal distribution.

            Signature
            ---------

            `(), (), () -> ()`

            Parameters
            ----------
            beta
                Shape parameter :math:`\beta` of the distribution.
            loc
                Location parameter :math:`\mu` of the distribution.
            scale
                Scale parameter :math:`\sigma` of the distribution. Must be
                positive.
            size
                Sample shape. If the given size is `(m, n, k)`, then `m * n * k`
                independent, identically distributed samples are returned. Default
                is `None` in which case a single sample is returned.

            """
            return super().__call__(beta, loc, scale, size=size, **kwargs)

        @classmethod
        def rng_fn_scipy(cls, rng, beta, loc, scale, size):
            return sps.gennorm.rvs(beta, loc=loc, scale=scale, random_state=rng, size=size)

gennorm = GenNormRV()

class GenNorm(Continuous):
    rv_op = gennorm

    @classmethod
    def dist(cls, beta, loc, scale, *args, **kwargs):
        beta = pt.as_tensor_variable(beta)
        loc = pt.as_tensor_variable(loc)
        scale = pt.as_tensor_variable(scale)
        return super().dist([beta, loc, scale], *args, **kwargs)

    def support_point(rv, size, beta, loc, scale):
        mu, _, _ = pt.broadcast_arrays(loc, scale, beta)
        if not rv_size_is_none(size):
            mu = pt.full(size, mu)
        return mu

    def moment(rv, size, beta, loc, scale):
        moment, _ = pt.broadcast_arrays(loc, scale, beta)
        if not pm.distributions.shape_utils.rv_size_is_none(size):
            moment = pt.full(size, moment)
        return moment

    def logp(value, beta, loc, scale):
        return check_parameters(
            pt.log(beta / (2 * scale))
            - pt.gammaln(1.0 / beta)
            - (pt.abs(value - loc) / scale) ** beta,
            beta > 0,
            scale > 0,
            msg="beta and scale must be positive",
        )

    def cdf(value, beta, loc, scale):
        b = value - loc
        res = pt.switch(
            np.asanyarray(b > 0),
            # Note the abs is required, because both switches are always calculated
            1 - 0.5 * pt.gammaincc(1.0 / beta, pt.abs(b / scale) ** beta),
            0.5 * pt.gammaincc(1.0 / beta, pt.abs(-b / scale) ** beta),
        )
        return check_parameters(
            res, beta > 0, scale > 0, msg="beta and scale must be positive"
        )

    def logcdf(value, beta, loc, scale):
        return pt.log(GenNorm.cdf(value, beta, loc, scale))

    def icdf(value, beta, loc, scale):
        res = pt.switch(
            np.asanyarray(value > 0.5),
            loc + scale * (pt.gammainccinv(1.0 / beta, 2 - 2 * value) ** (1.0 / beta)),
            loc - scale * (pt.gammainccinv(1.0 / beta, 2 * value) ** (1.0 / beta)),
        )
        res = check_icdf_value(res, value)
        return check_icdf_parameters(
            res, beta > 0, scale > 0, msg="beta and scale must be positive"
        )

[ricardoV94]

Thanks a ton @shifdz

If anybody wants to do it, please add it to https://github.com/pymc-devs/pymc-experimental