Writing a new distribution

Questions
1

Hi
I’m going to write a distribution which is not in the pymc3’s library. How should I do that? and how should I incorporate that in the codes?? (For example, suppose that Cauchy distribution is not available and I want to write that)

2

There is a very crude guide in here: Probability Distributions in PyMC3 — PyMC3 3.10.0 documentation

You have two options, to create a standard distribution by inheriting from the right subclasses (Discrete, Continuous, PositiveContinuous, etc), or using a pm.Densitydist that takes a callable logp and optional random as arguments.

The later approach might look like:

def my_cauchy(alpha, beta):
  def logp(value):
    # NotImplemented: Should return -np.inf if beta is negative
    log_density = -tt.log(np.pi) -tt.log(beta) -tt.log(1 + ((value - alpha) / beta) ** 2)
    return log_density
return logp

with pm.Model() as m:
  alpha = pm.Normal('alpha', 0, 2)
  beta = pm.Gamma('beta', 2, 2)
  likelihood = pm.Densitydist('likelihood', logp=my_chauchy(alpha, beta), observed=data)

For the former approach the best is to look at how is done in the source code: pymc3/continuous.py at v3 · pymc-devs/pymc3 · GitHub

2 Likes
3

so helpful… thanks

4

Just out of curiosity: what distribution are you planning to create? If commonly used it could end up in the code base, I guess. Also, I’m just about to start writing a few myself (astronomy inspired: power-law distributions, broken power-laws and Schechter functions).

5

I’m exactly trying to write power-law and Schechter distributions!! :slight_smile:
They are so common in astronomy and should be available in pymc3…

6

If you want to collaborate, hit me up!

1 Like
7

Hi, by power-law distribution, do you mean fitting a power law to a power spectra, e.g. a Fourier power spectrum? Also curious as this is what I would like to do.

8

I keep meaning to add the Generalized Extreme Value distribution too - is this of interest to others? I got as far as realizing that because of theano, it is not as simply as using the SciPy implementation and so it’ll be option 2 for that one.

1 Like
9

What we meant here was to implement power-law distribution functions (P(X) ~X**alpha).

2 Likes
10

Generalized EV might be useful - I’ve written an InverseWeibull for a recent project that I intend to PR once I get a minute!

2 Likes
11

for curious readers!, this is related to the topic