Skip to content

Commit 3add916

Browse files
ENH: Moyal distribution (#3870)
* Added Moyal distribution to continuous class. Added tests for Moyal distribution. * Release notes updated. * Added float32 exception for Moyal dist in tests.
1 parent 18a2c3b commit 3add916

File tree

5 files changed

+154
-3
lines changed

5 files changed

+154
-3
lines changed

RELEASE-NOTES.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
- `sample_posterior_predictive` can now feed on `xarray.Dataset` - e.g. from `InferenceData.posterior`. (see [#3846](https://github.com/pymc-devs/pymc3/pull/3846))
1313
- `SamplerReport` (`MultiTrace.report`) now has properties `n_tune`, `n_draws`, `t_sampling` for increased convenience (see [#3827](https://github.com/pymc-devs/pymc3/pull/3827))
1414
- `pm.sample` now has support for adapting dense mass matrix using `QuadPotentialFullAdapt` (see [#3596](https://github.com/pymc-devs/pymc3/pull/3596), [#3705](https://github.com/pymc-devs/pymc3/pull/3705) and [#3858](https://github.com/pymc-devs/pymc3/pull/3858))
15+
- `Moyal` distribution added (see [#3870](https://github.com/pymc-devs/pymc3/pull/3870)).
1516

1617
### Maintenance
1718
- Remove `sample_ppc` and `sample_ppc_w` that were deprecated in 3.6.

pymc3/distributions/__init__.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@
4848
from .continuous import LogitNormal
4949
from .continuous import Interpolated
5050
from .continuous import Rice
51+
from .continuous import Moyal
5152

5253
from .discrete import Binomial
5354
from .discrete import BetaBinomial
@@ -170,6 +171,7 @@
170171
'Interpolated',
171172
'Bound',
172173
'Rice',
174+
'Moyal',
173175
'Simulator',
174176
'fast_sample_posterior_predictive'
175177
]

pymc3/distributions/continuous.py

Lines changed: 132 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@
4040
'HalfCauchy', 'Gamma', 'Weibull', 'HalfStudentT', 'Lognormal',
4141
'ChiSquared', 'HalfNormal', 'Wald', 'Pareto', 'InverseGamma',
4242
'ExGaussian', 'VonMises', 'SkewNormal', 'Triangular', 'Gumbel',
43-
'Logistic', 'LogitNormal', 'Interpolated', 'Rice']
43+
'Logistic', 'LogitNormal', 'Interpolated', 'Rice', 'Moyal']
4444

4545

4646
class PositiveContinuous(Continuous):
@@ -1946,7 +1946,7 @@ class StudentT(Continuous):
19461946
plt.show()
19471947
19481948
======== ========================
1949-
Support :math:``x \in \mathbb{R}``
1949+
Support :math:`x \in \mathbb{R}`
19501950
======== ========================
19511951
19521952
Parameters
@@ -4381,3 +4381,133 @@ def logp(self, value):
43814381
TensorVariable
43824382
"""
43834383
return tt.log(self.interp_op(value) / self.Z)
4384+
4385+
4386+
class Moyal(Continuous):
4387+
R"""
4388+
Moyal log-likelihood.
4389+
4390+
The pdf of this distribution is
4391+
4392+
.. math::
4393+
4394+
f(x \mid \mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}\left(z + e^{-z}\right)},
4395+
4396+
where
4397+
4398+
.. math::
4399+
4400+
z = \frac{x-\mu}{\sigma}.
4401+
4402+
.. plot::
4403+
4404+
import matplotlib.pyplot as plt
4405+
import numpy as np
4406+
import scipy.stats as st
4407+
plt.style.use('seaborn-darkgrid')
4408+
x = np.linspace(-10, 20, 200)
4409+
mus = [-1., 0., 4.]
4410+
sigmas = [2., 2., 4.]
4411+
for mu, sigma in zip(mus, sigmas):
4412+
pdf = st.moyal.pdf(x, loc=mu, scale=sigma)
4413+
plt.plot(x, pdf, label=r'$\mu$ = {}, $\sigma$ = {}'.format(mu, sigma))
4414+
plt.xlabel('x', fontsize=12)
4415+
plt.ylabel('f(x)', fontsize=12)
4416+
plt.legend(loc=1)
4417+
plt.show()
4418+
4419+
======== ==============================================================
4420+
Support :math:`x \in (-\infty, \infty)`
4421+
Mean :math:`\mu + \sigma\left(\gamma + \log 2\right)`, where :math:`\gamma` is the Euler-Mascheroni constant
4422+
Variance :math:`\frac{\pi^{2}}{2}\sigma^{2}`
4423+
======== ==============================================================
4424+
4425+
Parameters
4426+
----------
4427+
mu: float
4428+
Location parameter.
4429+
sigma: float
4430+
Scale parameter (sigma > 0).
4431+
"""
4432+
4433+
def __init__(self, mu=0, sigma=1., *args, **kwargs):
4434+
self.mu = tt.as_tensor_variable(floatX(mu))
4435+
self.sigma = tt.as_tensor_variable(floatX(sigma))
4436+
4437+
assert_negative_support(sigma, 'sigma', 'Moyal')
4438+
4439+
self.mean = self.mu + self.sigma * (np.euler_gamma + tt.log(2))
4440+
self.median = self.mu - self.sigma * tt.log(2 * tt.erfcinv(1 / 2)**2)
4441+
self.mode = self.mu
4442+
self.variance = (np.pi**2 / 2.0) * self.sigma**2
4443+
4444+
super().__init__(*args, **kwargs)
4445+
4446+
def random(self, point=None, size=None):
4447+
"""
4448+
Draw random values from Moyal distribution.
4449+
4450+
Parameters
4451+
----------
4452+
point: dict, optional
4453+
Dict of variable values on which random values are to be
4454+
conditioned (uses default point if not specified).
4455+
size: int, optional
4456+
Desired size of random sample (returns one sample if not
4457+
specified).
4458+
4459+
Returns
4460+
-------
4461+
array
4462+
"""
4463+
mu, sigma = draw_values([self.mu, self.sigma], point=point, size=size)
4464+
return generate_samples(stats.moyal.rvs, loc=mu, scale=sigma,
4465+
dist_shape=self.shape,
4466+
size=size)
4467+
4468+
def logp(self, value):
4469+
"""
4470+
Calculate log-probability of Moyal distribution at specified value.
4471+
4472+
Parameters
4473+
----------
4474+
value: numeric
4475+
Value(s) for which log-probability is calculated. If the log probabilities for multiple
4476+
values are desired the values must be provided in a numpy array or theano tensor
4477+
4478+
Returns
4479+
-------
4480+
TensorVariable
4481+
"""
4482+
scaled = (value - self.mu) / self.sigma
4483+
return bound((-(1 / 2) * (scaled + tt.exp(-scaled))
4484+
- tt.log(self.sigma)
4485+
- (1 / 2) * tt.log(2 * np.pi)), self.sigma > 0)
4486+
4487+
def _repr_latex_(self, name=None, dist=None):
4488+
if dist is None:
4489+
dist = self
4490+
sigma = dist.sigma
4491+
mu = dist.mu
4492+
name = r'\text{%s}' % name
4493+
return r'${} \sim \text{{Moyal}}(\mathit{{mu}}={},~\mathit{{sigma}}={})$'.format(name,
4494+
get_variable_name(mu),
4495+
get_variable_name(sigma))
4496+
4497+
def logcdf(self, value):
4498+
"""
4499+
Compute the log of the cumulative distribution function for Moyal distribution
4500+
at the specified value.
4501+
4502+
Parameters
4503+
----------
4504+
value: numeric
4505+
Value(s) for which log CDF is calculated. If the log CDF for multiple
4506+
values are desired the values must be provided in a numpy array or theano tensor.
4507+
4508+
Returns
4509+
-------
4510+
TensorVariable
4511+
"""
4512+
scaled = (value - self.mu) / self.sigma
4513+
return tt.log(tt.erfc(tt.exp(-scaled / 2) * (2**-0.5)))

pymc3/tests/test_distributions.py

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@
3030
Bound, Uniform, Triangular, Binomial, SkewNormal, DiscreteWeibull,
3131
Gumbel, Logistic, OrderedLogistic, LogitNormal, Interpolated,
3232
ZeroInflatedBinomial, HalfFlat, AR1, KroneckerNormal, Rice,
33-
Kumaraswamy
33+
Kumaraswamy, Moyal
3434
)
3535

3636
from ..distributions import continuous
@@ -1213,6 +1213,13 @@ def test_rice(self):
12131213
self.pymc3_matches_scipy(Rice, Rplus, {'b': Rplus, 'sigma': Rplusbig},
12141214
lambda value, b, sigma: sp.rice.logpdf(value, b=b, loc=0, scale=sigma))
12151215

1216+
@pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
1217+
def test_moyal(self):
1218+
self.pymc3_matches_scipy(Moyal, R, {'mu': R, 'sigma': Rplusbig},
1219+
lambda value, mu, sigma: floatX(sp.moyal.logpdf(value, mu, sigma)))
1220+
self.check_logcdf(Moyal, R, {'mu': R, 'sigma': Rplusbig},
1221+
lambda value, mu, sigma: floatX(sp.moyal.logcdf(value, mu, sigma)))
1222+
12161223
@pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
12171224
def test_interpolated(self):
12181225
for mu in R.vals:

pymc3/tests/test_distributions_random.py

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -466,7 +466,12 @@ class TestGeometric(BaseTestCases.BaseTestCase):
466466
distribution = pm.Geometric
467467
params = {'p': 0.5}
468468

469+
470+
class TestMoyal(BaseTestCases.BaseTestCase):
471+
distribution = pm.Moyal
472+
params = {'mu': 0., 'sigma': 1.}
469473

474+
470475
class TestCategorical(BaseTestCases.BaseTestCase):
471476
distribution = pm.Categorical
472477
params = {'p': np.ones(BaseTestCases.BaseTestCase.shape)}
@@ -821,6 +826,12 @@ def ref_rand(size, mu, sigma):
821826
return expit(st.norm.rvs(loc=mu, scale=sigma, size=size))
822827
pymc3_random(pm.LogitNormal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand)
823828

829+
def test_moyal(self):
830+
def ref_rand(size, mu, sigma):
831+
return st.moyal.rvs(loc=mu, scale=sigma, size=size)
832+
pymc3_random(pm.Moyal, {'mu': R, 'sigma': Rplus}, ref_rand=ref_rand)
833+
834+
824835
@pytest.mark.xfail(condition=(theano.config.floatX == "float32"), reason="Fails on float32")
825836
def test_interpolated(self):
826837
for mu in R.vals:

0 commit comments

Comments
 (0)