Bring back distribution moments

[ricardoV94]

With #4983 and #5087 finished, it's a good time to bring back moments for our distributions for more stable starting points.

With this we can then do the switch in #5009

We should also update the distribution developer guide to mention the implementation of moments: https://github.com/pymc-devs/pymc/blob/main/docs/source/developer_guide_implementing_distribution.md

---

How to help?

1. This PR should give a template on how to implement and test new moments for distributions: https://github.com/pymc-devs/pymc/pull/5087/files

2. In most cases we should be able to copy the moments we were using in the V3 branch. For example here is what we were doing for the Beta
   

   pymc/pymc3/distributions/continuous.py

   Line 1235 in efbacce

   self.mean = self.alpha / (self.alpha + self.beta)

   
   2.1 We used to have multiple moments for some distributions such as mean, median, mode. We only support one moment now, and probably the "higher-order" one is the most useful (that is mean > median > mode)... You might need to truncate the moment if you are dealing with a discrete distribution.
   2.2 We left some of these moments commented out inside the distribution dist classmethod. Make sure they are removed when you implement them!

   pymc/pymc/distributions/continuous.py

   Lines 538 to 539 in 8f3636d

   	# mean = median = mode = mu = at.as_tensor_variable(floatX(mu))
   	# variance = 1.0 / self.tau

3. We have to be careful with size != None and broadcasting properly when when some parameters that are not used in the moment may nevertheless inform about the shape of the distribution. E.g. pm.Normal.dist(mu=0, sigma=np.arange(1, 6)) returns a moment of [mu, mu, mu, mu, mu]. Again Add tests for distributions moments #5087 should give some template to think about this

   pymc/pymc/distributions/continuous.py

   Lines 546 to 550 in 8f3636d

   	def get_moment(rv, size, mu, sigma):
   	mu, _ = at.broadcast_arrays(mu, sigma)
   	if not rv_size_is_none(size):
   	mu = at.full(size, mu)
   	return mu

   
   3.1 In the case where you have to manually broadcast the parameters with each other it's important to add test conditions that would fail if you were not to do that. A straightforward way to do this is to make the used parameter a scalar, the unused one(s) a vector (one at a time) and size None

4. Just to keep things uniformish, please add the get_moment immediately below the dist classmethod and before logp or logcdf methods.

5. New tests have to be added in test_distributions_moments.py. Make sure to test different combinations of size and broadcasting to cover the cases mentioned in point 3.

6. Don't hesitate to ask any questions. You can grab as many distributions to implement moments as you want. Just make sure to write in this issue so that we can keep track of it.

7. Profit with your new open source KARMA!

The following distributions don't have a moment method implemented:

• pymc.distributions.continuous.Beta Add Beta moment and tests for #5078 #5145
• pymc.distributions.continuous.Kumaraswamy Add back some missing moments #5147
• pymc.distributions.continuous.Exponential Add back some missing moments #5147
• pymc.distributions.continuous.Laplace Add back some missing moments #5147
• pymc.distributions.continuous.StudentT Add back some missing moments #5147
• pymc.distributions.continuous.Cauchy Add back some missing moments #5147
• pymc.distributions.continuous.HalfCauchy HalfCauchy, Gamma, Weibull and LogNormal moments #5148
• pymc.distributions.continuous.Gamma HalfCauchy, Gamma, Weibull and LogNormal moments #5148
• pymc.distributions.continuous.Weibull HalfCauchy, Gamma, Weibull and LogNormal moments #5148
• pymc.distributions.continuous.LogNormal HalfCauchy, Gamma, Weibull and LogNormal moments #5148
• pymc.distributions.continuous.HalfStudentT Add HalfStudentT moment #5152
• pymc.distributions.continuous.ChiSquared Add ChiSquared moment #5154
• pymc.distributions.continuous.Wald Add Pareto and Wald moment #5161
• pymc.distributions.continuous.Pareto Add Pareto and Wald moment #5161
• pymc.distributions.continuous.InverseGamma Adds moments for inverse gamma distribution #5199
• pymc.distributions.continuous.ExGaussian Add ExGaussian moment #5165
• pymc.distributions.continuous.VonMises Vonmises moments #5232
• pymc.distributions.discrete.Binomial Binomial and Poisson Moment #5150
• pymc.distributions.discrete.BetaBinomial adds beta-binomial mean and test cases #5175
• pymc.distributions.discrete.Poisson Binomial and Poisson Moment #5150
• pymc.distributions.discrete.NegativeBinomial NegativeBinomial, ZeroInflatedPoisson, ZeroInflatedBinomial moments #5163
• pymc.distributions.discrete.Constant Add Constant Moment #5156
• pymc.distributions.discrete.ZeroInflatedPoisson NegativeBinomial, ZeroInflatedPoisson, ZeroInflatedBinomial moments #5163
• pymc.distributions.discrete.ZeroInflatedNegativeBinomial Corrected ZeroInflatedBinomial moment and addeed ZeroInflatedNegativeBinomial moment #5206
• pymc.distributions.discrete.ZeroInflatedBinomial NegativeBinomial, ZeroInflatedPoisson, ZeroInflatedBinomial moments #5163
• pymc.distributions.discrete.DiscreteUniform Discrete uniform and hyper geometric moment #5167
• pymc.distributions.discrete.Geometric Added moment for geometric distribution along with test for issue #5078 #5158
• pymc.distributions.discrete.HyperGeometric Discrete uniform and hyper geometric moment #5167
• pymc.distributions.discrete.Categorical add categorical moment #5176
• pymc.distributions.distribution.DensityDist Add DensityDist moment #5159
• pymc.distributions.multivariate.MvNormal Add MvNormal moment #5171
• pymc.distributions.multivariate.MatrixNormal Add MvStudentT and MatrixNormal moment #5173
• pymc.distributions.multivariate.KroneckerNormal Add moments for KroneckerNormalDistribution #5235
• pymc.distributions.multivariate.MvStudentT Add MvStudentT and MatrixNormal moment #5173
• pymc.distributions.multivariate.Dirichlet Adding Dirichlet moment and tests #5174
• pymc.distributions.multivariate.Multinomial Add Multinomial moments #5201
  • Note: This distribution had tests for the mode in test_distributions.py which should be moved/refactored to test_distributions_moments.py

  • pymc/pymc/tests/test_distributions.py

    Line 2151 in bdd4d19

    def test_multinomial_mode(self, p, n):

  • pymc/pymc/tests/test_distributions.py

    Line 2189 in bdd4d19

    def test_multinomial_mode_with_shape(self):

• pymc.distributions.multivariate.DirichletMultinomial Adds tests and mode for dirichlet multinomial distribution #5225
• pymc.distributions.continuous.AsymmetricLaplace Adding moments for AsymmetricLaplace and SkewNormal #5188
• pymc.distributions.continuous.SkewNormal Adding moments for AsymmetricLaplace and SkewNormal #5188
• pymc.distributions.continuous.Triangular Added moments for gumbel, triangular and logitnormal distributions for issue #5078 #5180
• pymc.distributions.discrete.DiscreteWeibull Add DiscreteWeibull moment #5407
• pymc.distributions.continuous.Gumbel Added moments for gumbel, triangular and logitnormal distributions for issue #5078 #5180
• pymc.distributions.continuous.Logistic Add Logistic Moment #5157
• pymc.distributions.continuous.LogitNormal Added moments for gumbel, triangular and logitnormal distributions for issue #5078 #5180
• pymc.distributions.continuous.Interpolated Adding Interpolated moment #5222
• pymc.distributions.continuous.Rice Adding moment for Rice distribution and associated test #5190
• pymc.distributions.continuous.Moyal Adding moment for Moyal distribution and corresponding tests #5179
• pymc.distributions.simulator.Simulator Add moment for Simulator distribution #5208
• pymc.distributions.multivariate.CAR Add moments for CAR distribution #5220
• pymc.distributions.continuous.PolyaGamma Adding polya gamma moment and corresponding tests #5193
• pymc.bart.BART add moment for BART distribution #5211

[ricardoV94]

@twiecki I updated #5087 which includes a possible template for adding and testing moments. Once we merge that, we can use it as an example

[ricardoV94]

Updated the top comment with a list of missing moments and some instructions

[michaeloriordan]

I'll add a PR for pymc.distributions.continuous.Beta as a warmup #5145. Will work on some more if I got this one right! 😄

[michaeloriordan]

I'll work on the next 5 distributions on the list. PR here #5147

• pymc.distributions.continuous.Kumaraswamy
• pymc.distributions.continuous.Exponential
• pymc.distributions.continuous.Laplace
• pymc.distributions.continuous.StudentT
• pymc.distributions.continuous.Cauchy

[lucianopaz]

I'll do:

• HalfCauchy
• Gamma
• Weibull
• LogNormal

[farhanreynaldo]

I would like to take on:

• Binomial
• Poisson

[twiecki]

Thanks @farhanreynaldo!

[azihna]

I'd like to work on HalfStudentT

[ricardoV94]

I'd like to work on HalfStudentT

Sure go ahead :)

[patel-zeel]

I have added ChiSquared moment.

[zoj613]

Yeah, as far as I understand. The CAR is modeled using a multivariate normal distribution with a specific covariance matrix. So it should be identical to the MvNormal if both accept the same kind of parameters mu and cov.

Sounds good. Do you want to open a PR for it? I'll mark it above as taken.

Sure, i'll open one.

[ricardoV94]

All distributions have been taken. I've just unchecked those that don't yet have a PR open to facilitate tracking in the final stage.

[ricardoV94]

I would like to try VonMises

@Domenico89 Any progress on this?

[ricardoV94]

I will add KroneckerNormal, Wishart, DiscreteWeibull.

@anirudhacharya Any progress on this?

[Domenico89]

I would like to try VonMises

@Domenico89 Any progress on this?

@ricardoV94 I should have an implementation for it, but I am getting an error for some values of the test parameters for which I haven't managed to figure out the reason yet

[ricardoV94]

I would like to try VonMises

@Domenico89 Any progress on this?

@ricardoV94 I should have an implementation for it, but I am getting an error for some values of the test parameters for which I haven't managed to figure out the reason yet

Feel free to open a PR with what you have so far and we will help you there

[anirudhacharya]

@anirudhacharya Any progress on this?

I will finish up the tests and raise a PR by tomorrow.

[anirudhacharya]

@ricardoV94 Wishart Distribution documentation says - This distribution is unusable in a PyMC3 model should I be adding tests and moments for it? Is there an example usage I can look at for writing the tests?

[twiecki]

@anirudhacharya No, we don't need that for Wishart.

[anirudhacharya]

@ricardoV94 I am having trouble writing tests for the DiscreteWeibull distribution.

Is there some example usage of this distribution that I can check out? The documentation did not provide any - https://docs.pymc.io/en/v3/api/distributions/discrete.html?highlight=Weibull#pymc3.distributions.discrete.DiscreteWeibull

[ricardoV94]

@ricardoV94 I am having trouble writing tests for the DiscreteWeibull distribution.

Is there some example usage of this distribution that I check out? The documentation did not provide any - https://docs.pymc.io/en/v3/api/distributions/discrete.html?highlight=Weibull#pymc3.distributions.discrete.DiscreteWeibull

Feel free to open a PR with what you have so far, it will be easier to help you then.

I don't know about any examples using the distribution, but the check_logp test gives an indication of some values and parameters that are being tested successfully:

pymc/pymc/tests/test_distributions.py

Lines 1740 to 1746 in a50b386

	def test_discrete_weibull(self):
	self.check_logp(
	DiscreteWeibull,
	Nat,
	{"q": Unit, "beta": NatSmall},
	discrete_weibull_logpmf,
	)

[ricardoV94]

@anirudhacharya are you still interested in the DiscreteWeibull or should we put it up for grabs again?

[anirudhacharya]

@anirudhacharya are you still interested in the DiscreteWeibull or should we put it up for grabs again?

Please put it up for grabs.

[zoj613]

@ricardoV94 The mean looks like an infinite sum. This is a job for scan (seeing that the computation is meant to be symbolic) with a terminal condition based on tolerance, right? If so I can take a stab at it.

[ricardoV94]

@zoj613 we are not interested in the exact moments, but just good enough moments to start sampling. It seems the mean is quite a costly computation, so we should pick something else, perhaps the median or the mode? All we care is that the initial point always has a non-zero probability and, provided that, it is in a place with relevant prior mass

[zoj613]

Okay I see. If my math is correct it looks like the median can be computed from the CDF expression as:

median = at.power(at.log(0.5) / at.log(q), 1 / beta) - 1

where CDF = 1 - q ** ((x + 1) ** beta)

[ricardoV94]

Okay I see. If my math is correct it looks like the median can be computed from the CDF expression as:

median = at.power(at.log(0.5) / at.log(q), 1 / beta) - 1

where CDF = 1 - q ** ((x + 1) ** beta)

Great! Feel free to open a PR for it