Rice distribution is misdefined

[nbud]

Description of your problem

The Rice distribution Rice(nu, sigma) is the distribution of sqrt(X^2+Y^2) where X~N(nu, sigma^2) and Y~N(nu, sigma^2), with X and Y independent. This parametrisation is the one notably used in the Wikipedia article on Rice distribution.

scipy.stats.rice used a different distribution: shape parameter b=nu/sigma and scale sigma.

The current Rice distribution in pymc inconsistently mixes the two parametrisations: nu is sometimes the location parameter of the normal distributions, sometimes the shape parameter.

import pymc3 as pm
from scipy import stats
import matplotlib.pyplot as plt
import numpy as np

sigma = 5
nu = 3

u = pm.distributions.Rice.dist(nu, sigma)
v = stats.rice(b=nu/sigma, scale=sigma)

np.random.seed(123)
mean_u = u.mean.eval()

n = 10000
samples_u = u.random(size=n)
samples_v = v.rvs(size=n)

plt.figure()
plt.hist(samples_u, density=True, bins=50, label="pymc3.Rice")
plt.hist(samples_v, density=True, bins=50, label="scipy Rice dist")
plt.axvline(mean_u, label="pymc3.Rice.mean", color="C2")
plt.axvline(np.mean(samples_u), label="empirical mean of pymc3.Rice", color="C3")
plt.legend()

assert np.isclose(mean_u, v.mean())  # pass

I'm writing a PR to fix this.

[twiecki]

Good catch!

[nbud]

Closed by #3289, see also #3287