Remove degenerate cases

In some places, like the Gamma function, some effort is made to deal with degenerate
parameters like infinity. I think this is too much effort for a minute use case, as
I don't think users will insert these values to query the library to see what
degenerate distributions look like.

Having to deal with the special cases reduces clarity and complicates the API
considerably, as it forces to contend with mathematical details difficult to convey
in a numerical library. As the limit α -> inf is approached, the pointwise
limit of the pdf and the cdf is zero, so these aren't pdfs and cdfs anymore. As the
limit β -> inf is approached, the continuous distribution degenerates to a discrete
one: P(X = 0) = 1.

I propose to remove all special cases and document why they aren't handled.
This would deal with #57 and #98.

[boxtown]

I think we'd have to evaluate this on a case-by-case basis for distributions but I agree I think at least for the Gamma distribution it does not necessarily make sense to allow INF for α or β. I'll create an issue specifically for the Gamma distribution

[dhardy]

FYI, we've recently had someone interested in supporting a particular degenerate case of Exp.

The use case referenced to justify degenerate parameter values is likely erroneous, an absorbing state can be modelled by having all outgoing paths loop back to the same state. In that case any distribution would work, unless one is interested in statistics about how many times a transition, including transitions to the same state, occurred.

I disagree with the proposed solution and agree with the Wikipedia page that the case is degenerate for a reason.
To argue this, mathematical infinity is more often than not not an actual element of the probability space but a statement of what happens when a limit is taken. If a distribution is dependent on a parameter α then the function can be seen as depending on two variables p(α,x). If any IEEE infinities or other degenerate cases are to be allowed, then one would have to either

1. treat all possible permutations of taking the limit and find a way to query the library for each of them, as they're unequal in general
2. agree to take the limit in one specific order
3. some mixture of 1. and 2.

Additionally to these complications, as stated above, either cdf can exist but pdf be a non-classical function as is the case with a discrete probability or both limits exist but aren't cdf's and pdf's anymore. If degenerate cases are allowed where not all the functions exist in a classical sense there's an argument to be made that returning errors is more helpful than approximating a delta function with p(x0) = infinity and for all other x p(x) = 0

Perhaps I got lost in the mathematics on this one, I take some of that back. If there's only one parameter apart from the variable of interest then I agree, More general Exp is a reasonable extension. In the case of the gamma function however there are two parameters, so knowing in which order to take the limit doesn't have an equally reasonable default.
What also simplifies the resolution in the rand crate is that there no evaluation of pdf is sought, only a consistent extension of sampling.

[dhardy]

Thanks for the analysis and follow up comment. I agree: the exp case seems unambiguous, but in general that is not the case. It may be that allowing for degenerates in only one variable is also a sensible option (but I haven't thought enough about the gamma case).

[saona-raimundo]

Author of the More general Exp issue here.

What also simplifies the resolution in the rand crate is that there no evaluation of pdf is sought, only a consistent extension of sampling.

Exactly!! This is what I had in mind. In the case of Exp, it would be an expected behaviour to sample infinity.

In the other hand, a Gamma with both parameters infinity has no intuitive behaviour.
Having said that, one should either document the chosen behaviour (as it has been done in statrs) or not support it to start with.

I think we'd have to evaluate this on a case-by-case basis for distributions

I totally agree with that, but not only because there can be an intuitive answer for some cases but also because there can be useful (and arbitrary) implementations. Just as an example: assert_eq!(std::f64::INFINITY, 1.0 / 0.0) this behaviour simplified our code, but there is no correct answer to what 1 / 0 is.

Just as an example: assert_eq!(std::f64::INFINITY, 1.0 / 0.0) this behaviour simplified our code, but there is no correct answer to what 1 / 0 is.

As I outlined with the mathematical infinities that get confused with IEEE754 infinities, the semantics are not identical. This starts to get unintuitive because -0.0 == 0.0 which respects the mathematics but 1.0 / -0.0 != 1.0 / 0.0, so in some cases the floating point numbers are supposed to be interpreted as limits while in other contexts this isn't the case. So yes 1 / 0 does not have a sensible solution as a number satisfying 0 * a = 1, but that's not the angle to take. Floating point expressions can either be made sensible as numbers, made sensible as limits or neither of those, which leads to the at first unintuitive property of NaNs not being equal to themselves. That's akin to the situation with functions here.

[dhardy]

Yes, but you missed the bit about convergence: 1/x → +inf as x → 0 from above. "INF" is not a real mathematical value, but a model of what you observe when approaching a degenerate value.

Based on the discussion my current proposal would be

• Allow singly degenerate parameter values, limit taken first over parameter, then over argument (x in the above thread).
• Future-proofing for multidimensional arguments by also only allowing singly degenerate argument evaluations
• Return NaN if any function limits don't exist as classical functions.

I think an error is preferable to NaN on the last point, but this is at least not immediately desirable due to to being backwards-incompatible.

[saona-raimundo]

I agree with these principles and would add the following:

• If a multi-parameter case has defined behaviour, there should not be a problem in adding it.

I am not sure the consequences of following these premises, but could start with something:

1. Gamma distribution should not accept both parameters (alpha and beta) as infinity. But, it should allow one of them to be infinity.
2. Normal distribution should not accept mean infinity and std_dev infinity at the same time. But, should allow: mean infinity and std_dev finite (returning infinity always); mean finite and zero std_dev (returning a constant equal to mean always); mean finite and infinity std_dev should not be allowed, as it does not converge in a reasonable way.
3. ...

Of course, the list can be longer and implementing all these changes might not be what we want(?)
The possibility is to implement them "as it is requested".

• If a multi-parameter case has defined behaviour, there should not be a problem in adding it.

I agree, caveats below, especially since errors are returned when constructing the distributions, it's an immediate feedback on what to expect.

2. Normal distribution should not accept mean infinity and std_dev infinity at the same time. But, should allow: mean infinity and std_dev finite (returning infinity always); mean finite and zero std_dev (returning a constant equal to mean always); mean finite and infinity std_dev should not be allowed, as it does not converge in a reasonable way.

This is at odds with the solution in the rand crate. Finite mean and infinity std_dev do not make the distribution converge towards any other, but neither did the exponential distribution. This would deny extending some properties of the distribution in a meaningful way, like the mean and the std_dev, which do have sensible limits. More general Exp does not care whether the distribution converges, only if one selected property can be extended consistently. I would even argue that using the logic used to extend the exponential clock one would conclude that the limit of sampling of a normal distribution with mean 0 and infinite variance is infinity with probability 1/2 and -infinity with probability 1/2, so does have a limit.

The possibility is to implement them "as it is requested".

Strongly disagree. We should model the established mathematical practices as closely as possible. Convergence of random variables have formal methodologies for characterising random variables are degenerate cases of each other. I thus reverse my reversal on the issue and think more general Exp is a bad idea due to selectively extending properties derived from a distribution that cannot formally be shown to converge against any distribution. It breaks mathematics in subtle ways. The principled approach would be to pick one of the convergence metrics and extend to degenerate cases if such a limit in the chosen metric exists, and deny deriving any properties from inexistent distribution limits.

Edit: big rewrite for clarification purposes

[saona-raimundo]

I agree, caveats below, especially since errors are returned when constructing the distributions, it's an immediate feedback on what to expect.

Yes, indeed! and any function that can fail should return a Result, that is the idea :)

Finite mean and infinity std_dev do not make the distribution converge

Exactly, that is why I proposed to not accept this case either:

mean finite and infinity std_dev should not be allowed, as it does not converge in a reasonable way.

Convergence

More general Exp does not care whether the distribution converges, only if one selected property can be extended consistently.

Maybe it was not so explicit, but the convergence we considered to solve More general Exp was "the intuitive behaviour a user would expect", which I think is still the correct thing to do.

To translate into a mathematical convergence one needs a sequence of random variables. Sadly, there is no "good sequence" to consider. Let's say we want to claim something like "Exp(0) is the almost sure limit of a sequence of Exp(lambda) when lambda goes to zero". One can prove that the almost sure convergence depends on the sequence of lambdas we consider. Therefore, there is no way we can put a criteria like "let's accept the limit when it corresponds to the almost sure convergence". (see Appendix)

Intuitive convergence

If intuitive convergence is not almost sure convergence, what should it be?
As you point out, the convergence of distributions is not a "good" convergence either, since More general Exp would not be right in this case.

I think that, if there is a characterization of Exp(lambda) = Exp(1) / lambda, then it is easy to agree in the limit when lambda goes to zero. Basing the limit cases in characterizations like this is, of course, relative to the particular characterization and the limit might "break" under other characterizations.

Usually, such a characterization is the one used when simulating the variable, so in the context of the rand crate is easy to get an agreement. Maybe, if such a characterization is also part of the documentation in statrs, then the expected behaviour will also be "intuitive".

Practicality

The possibility is to implement them "as it is requested".

Strongly disagree. We should model the established mathematical practices as closely as possible.

This is purely from a practical point of view. Of course one would like the crate to be as complete as possible.

Appendix

Consider a sequence of variables X_n with distribution Exp(lambda_n). We will show that, depending on the sequence (lambda_n)_{n > 0}, we can have that (X_n) converges almost surely to infinity or not.

Proof:

Case 1: not convergent.
First, consider a decreasing sequence (alpha_n) that goes to zero. The idea is to take increasingly more copies of alpha_n so that we can ensure that there is no convergence to infinity. Consider alpha_1, then, take m_1 copies of an Exp(alpha_1) so that the probability that the minimum among these variables is less than one is greater than (0.9)^{1 / 2^n}. Then go to alpha_2 and do the same to define a new number m_2, and so on for each alpha_n define m_n.

Consider the sequence (lambda_n) as the decreasing sequence of all alphas considered in the previous step, with the number of repetitions considered. The idea is that the sequence (Exp(lambda_n)) has so many chances to be below one, that the liminf of the sequence is not infinity and therefore its limit either. Indeed, consider n such that lambda_{n+1} < lambda_n, i.e. just after "changing to a lower alpha". Then, the probability that in the last m_n variables the minimum is less than one is at least (0.9)^{1 / 2^n}. Therefore, the liminf is less than one with probability one.

Case 2: convergent.

Consider a sequence of Exponential variables X_n with distribution Exp(1 / n^2).
We need to prove that the probability of lim_{n \to \infty} X_n = \infty is one.
Equivalently, we can prove that, for all K > 0, we have that the probability of liminf_{n \to \infty} X_n > K is one.

Define the event E_n = {w : X_n(w) <= K}. Note that (limsup E_n)^c is the event liminf X_n > K. Let us prove that the probability of limsup E_n is zero. The probability of E_n is 1 - e^{-K / n^2}, which is less than K / n^2 (linear expansion of exponential around zero). Using Borel-Cantelli's lemma, since the sum over the probability of E_n is finite, we have that the probability of limsup E_n is zero, which means that the probability of liminf X_n > K is one. Since K is arbitrary, we conclude that the liminf of (X_n) is infinity and therefore the sequence converges to infinity almost surely.