CDFs for discrete variables and PMFs only allow integer values

[ethan-alt]

Description

Functions such as bernoulli_cdf and bernoulli_pmf only accept integer values. This is mathematically imprecise, as CDFs and PMFs should be defined for all real numbers (e.g., pbinom and dbinom in R accept non-integer values).

This is a problem, e.g., if you have a matrix Y of values where each column has a different distribution (e.g., first column is Gaussian, second column is Bernoulli). Such constructs are necessary to create a general Gaussian copula function.

Example

binomial_cdf(0.5, 0.75);

Expected Output

0.25

Currently, no code will compile if a non-integer value is input as the first argument to binomial_cdf(0.5, 0.75).

Current Version:

v4.0.1

[bob-carpenter]

Thanks for reporting.

This is mathematically imprecise, as CDFs and PMFs should be defined for all real numbers (e.g., pbinom and dbinom in R accept non-integer values).

Stan's definition is decidedly non-standard. I'll take the blame with my obsession for type checking. I wasn't expecting anyone to check the values of a cdf outside the support of the pmf (or pdf), but this isn't the first time it's come up.

We could relax this with some work on the parser, code generator, math library, and documentation. We'd either have to make the cdf arguments non-differentiable or define derivatives for all of them (constant between support, undefined at support). That technically breaks the assumptions behind HMC, which requires smooth functions, but it can still be workable. And given that the derivatives are constant, they never play a role in guiding parameters.

Could you explain more about the use case where you expect the cdfs to be defined outside of the support of the pmf? An example you'd like to write in Stan that doesn't work because of this constraint would be most compelling. Do you need derivatives (they'll all be zero and they're not defined at the support of the pmf because of the discontinuity?

[ethan-alt]

The specific use case is that maybe we have a matrix, Y, containing Gaussian, Binary, and Poisson values, so Y has to be a real matrix. For the discrete variables, one must evaluate the CDFs at y and y-1 for a Gaussian copula.

It's possible to work around it by having two input matrices, but sadly the solution for tying it all together becomes rather unclean.

Derivatives would not be necessary.

[bob-carpenter]

Thanks for the elaboration. I'm still trying to imagine what the Stan code would look like. The obstacle to encoding binary, count, and real-valued entries is that Stan doesn't have heterogeneous data structures like this. All matrices are real-valued in Stan and there's no way to convert real values to integers without losing the gradient information that drives Stan's HMC sampler.

Sorry to keep asking, but do you have a sketch of what the code would look like if Stan did extend the range of the CDFs to real values?

Depending on the distribution, a workaround might be to directly implement the cdf as a function in Stan. It's unlikely we'd be able to do this quickly even if the decision is made to change the support of cdfs.

After many years, I think you're right---we should've had the cdfs defined on all of R. It was my fault for defining them the way they are defined (and my excuse is that I was trained as a computer scientist, not a statistician). That still wouldn't solve problems with derivatives, but at least you'd be able to get the values.

[ethan-alt]

I think a function that allows for the evaluation of CDFs over the real line would work fine.

For background, I have attached some code that will not work, but what I was attempting to implement.

Basically, the user would input, say, a N x 3 matrix of normal, Bernoulli, and Poisson dependent variables (along with covariates, if applicable) and have an indicator of distribution type (e.g., 1 = normal, 2 = Bernoulli, 3 = Poisson).

For the purposes of this specific application, I've developed a workaround where the user must input separate arrays for each distribution type, but the code is not as clean. Still, I feel the option to at least evaluate the CDF over the real line would be a good implementation.

copula_glm.zip

[spinkney]

@bob-carpenter this PR to my helpful stan functions repo has the code and a reproducible R example. Same as attached by Ethan but possibly a bit cleaner.

spinkney/helpful_stan_functions@e603fb9

[betanalpha]

This is mathematically imprecise, as CDFs and PMFs should be defined for all real numbers (e.g., pbinom and dbinom in R accept non-integer values).

This is incorrect. Cumulative distribution functions and probability mass functions for distributions defined on discrete spaces are defined on only the discrete elements of the set. Even if one embeds those discrete spaces into the real numbers then the discrete probability mass function and cumulative distribution functions do not specify real-valued functions. Indeed there are an infinite number of real-valued functions that interpolate between any PMF/CDF at the discrete values for which they are defined. Many software implementations of common PMF/CDF families, such as those in R, exploit analytic interpolations of discrete functions, such as the gamma function interpolation of the discrete factorial function. Admitting real-valued inputs is an artifact of the implementation, not a mathematical consequence of discrete probability theory. A program that wants to overload discrete PMF/CDF families can just use the analytic real-valued interpolations directly — for example gamma_p and gamma_q are both exposed in the Stan modeling language.

[spinkney]

To @betanalpha's point this is pretty straightforward. Here's the overloaded binomial and poisson cdfs in R

> library(gsl)
> cpoisson_cdf <- function(x, lambda) {
+   gamma_inc_Q(x + 1, lambda) 
+ }
> 
> cpoisson_cdf(2, 5)
[1] 0.124652
> ppois(2, 5)
[1] 0.124652
> 
> cbinomial_cdf <- function(x, N, p) {
+   1 - beta_inc(x + 1, N - x, p)
+ }
> 
> cbinomial_cdf(2, 10, 0.3)
[1] 0.3827828
> pbinom(2, 10, 0.3)
[1] 0.3827828

[ethan-alt]

This is mathematically imprecise, as CDFs and PMFs should be defined for all real numbers (e.g., pbinom and dbinom in R accept non-integer values).
This is incorrect. Cumulative distribution functions and probability mass functions for distributions defined on discrete spaces are defined on only the discrete elements of the set. Even if one embeds those discrete spaces into the real numbers then the discrete probability mass function and cumulative distribution functions do not specify real-valued functions. Indeed there are an infinite number of real-valued functions that interpolate between any PMF/CDF at the discrete values for which they are defined. Many software implementations of common PMF/CDF families, such as those in R, exploit analytic interpolations of discrete functions, such as the gamma function interpolation of the discrete factorial function. Admitting real-valued inputs is an artifact of the implementation, not a mathematical consequence of discrete probability theory. A program that wants to overload discrete PMF/CDF families can just use the analytic real-valued interpolations directly — for example gamma_p and gamma_q are both exposed in the Stan modeling language.

It is not incorrect. The definition of a mass function for a random variable X evaluated at value x is "the probability that X is equal to x, i.e., p(x) = P(X = x). For example if X is a binomial random variable with probability 0.25, p(0) = 0.75, p(1) = 0.25, p(0.5) = 0 because P(X = 0.5) = 0 since X can only be 0 or 1.

For CDFs, F(x) = P(X <= x) is similarly always well defined. In the above example, F(-1) = 0 because X has to be at least 0, F(0) = 0.75 because P(X <= 0) = P(X = 0) = 0.75 and F(1.5) = P(X <= 1.5) = P(X <= 1) = 1 since X cannot be larger than 1.

[bob-carpenter]

The Wikipedia entry for probability mass functions has two definitions in the first two paragraphs. Their simple definition matches @ethan-alt's (although restricts to real numbers, not any old x) and their more measure-theoretic one matches @betanalpha's. Their definition of the Bernoulli PMF has exactly the domain {0, 1} and isn't extended to real numbers.

Setting all the mathematical generalizations aside, the question for our API isn't whether it's possible to extend the definition of a discrete pmf (or cdf) to all real numbers, but whether it's more helpful to restrict the domain or to catch potential errors. As a type-sensitive language designer, I went with the catch-errors-early approach. I'm not trying to argue it's the only way things could have been designed or even that it's the best design. I'm just explaining why I defined things that way to begin with.

Stan's underlying typing also requires an integer argument for a Bernoulli pmf. It won't compile if you try to give it a real number and will throw a runtime exception if it tries to evaluate the probability of an integer other than 0 or 1. I made that decision before Stan 1.0 to err on the side of catching errors as quickly as possible. Nevertheless, it's possible to define even a mixed continuous/discrete variable in Stan using mixtures (e.g., spike and slab prior with proper delta function spike).

[bob-carpenter]

(from @ethan-alt) Basically, the user would input, say, a N x 3 matrix of normal, Bernoulli, and Poisson dependent variables (along with covariates, if applicable) and have an indicator of distribution type (e.g., 1 = normal, 2 = Bernoulli, 3 = Poisson).

This is how R works and I can see where you're trying to go with this. I decided to make Stan strongly typed and to distinguish integer and real values like in C++ and Java (more like Java in that you can't assign real values to integers). You're not the only one who wished we'd typed Stan like R.

My motivation was that it's generally considered bad form in numerical programming to use floating point for discrete problems or to even try to compare floating point numbers with == (for example, there are two zero values in standard IEEE floating point, which are == but do not have the same behavior everywhere).

[spinkney]

@bob-carpenter can't we get around that by defining new functions for the continuous/relaxed version of discrete cdfs? So we could have a "relaxed_poisson_cdf" or "continuous_approx_poisson_cdf" or whatever name people like the most.

[nhuurre]

@ethan-alt: The definition of a mass function for a random variable X evaluated at value x is "the probability that X is equal to x, i.e., p(x) = P(X = x).

Technically true, but @betanalpha's point was that there are alternative "interpolations" that are commonly used. See for instance @spinkney's post right above yours that begins "this is pretty straightforward" and then gives the wrong definition. That sort of confusion can be avoided only by leaving the functions undefined outside of their natural domain.

The Stan language has enough combinatorial functions that you can build your own discontinuous CDFs if you need them.

real real_binomial_cdf(real x, int N, real p) {
  return 1 - inc_beta(floor(x)+1, N-floor(x), p);
}
real real_poisson_cdf(real x, real lambda) {
  return 1 - gamma_p(floor(x), lambda);
}

[spinkney]

It's not "wrong" it's the "continuous" or "interpolated" version that doesn't use floor...

[ethan-alt]

@bob-carpenter I agree whether PMFs should be extended along the real line is a matter of some debate (although it does mention in the article you posted that the Radon-Nikodym derivative is typically extended to the real line--so that even the measure-theoretic definition might include all real numbers), but it is not a debate for CDFs. CDFs are the only thing that discrete and continuous random variables have in common. It's fundamentally important that CDFs are defined on the entire real line. I will also say that the C code that powers dbinom and pbinom in R is set to be a double type. Nevertheless, I understand the discontinuity (no pun intended), and I'm not saying that what is currently done in Stan is "wrong," just that it "could be more accurate."

I certainly can understand that you don't want to break existing code, which is why I think @spinkney 's solution is a good one--perhaps have a separate function like binomial_cdf_real or something.

@nhuurre Of course, one can always build their own functions. But that's true of anything. Stan is not intended to be a general purpose language. It is intended to be utilized as a software package for people in applied statistics. As such, it makes sense that, at a minimum, the precise mathematical definitions of statistical functions are available.

I'm not on the dev team of course, so perhaps my vote counts the least in this discussion. But as an avid user and fan of Stan, who dabbles in statistical package development, I hope it counts for something.

[bob-carpenter]

(@spinkey): It's not "wrong" it's the "continuous" or "interpolated" version that doesn't use floor...

Is it intended as a pdf rather than pmf, and if so, is it normalized? I've never quite understood these tactics, though there's a really neat paper by Will Grathwohl et al. that leans on this idea to add gradient-like information to discrete distributions that looks very promising for things like count variables: https://arxiv.org/abs/2102.04509

(@ethan-alt) CDFs are the only thing that discrete and continuous random variables have in common. It's fundamentally important that CDFs are defined on the entire real line.

(@ethan-alt) I will also say that the C code that powers dbinom and pbinom in R is set to be a double type.

dbinom() returns NaN if the count parameter is not an integer and throws a warning and returns 0 if the outcome is not an integer. Arguably this is all just to get around that R isn't big on exceptions (aka conditions in R) and tries its hardest to not "fail" (which I consider a huge drawback to writing robust code in R, but that's another topic). pbinom() behaves like a usual CDF, and thus also produces NaN and prints a warning if the count parameter isn't an integer.

(@ethan-alt) I certainly can understand that you don't want to break existing code,

It'd break existing unit tests and doc, but I doubt changing either the pmfs (which I'm opposed to) or cdfs (which I'm OK with) would have much effect on user code. Behavior wouldn't change for existing programs that worked other than that rather than an exception out of domain you'd get a 0 or 1 depending on which direction you left the support. I don't think I've ever seen a user program that relies on failure out of domain to drive anything. The reason I put it there in the first place is that it's usually a coding mistake.

(@ethan-alt) I'm not on the dev team of course, so perhaps my vote counts the least in this discussion.

We try to listen to everyone. The devs get votes on changes like this, but I'm not even sure who counts as a dev any more. Our governing body is talking about extending that to users somehow, but I don't think that's in place yet.

(@ethan-alt). Of course, one can always build their own functions. But that's true of anything.

Our main concern is not adding so many functions that it becomes impossible to find the one you need or decide among the ones that are there. And that it doesn't become a maintenance nightmare. We have a lot of distributions, so changes like this can be a programming nightmare and maintenance burden long term. I even wrote a blog post about maintenance costs to try to explain maintenance costs to statisticians. Normally we could introduce new functions and deprecate old ones, but I don't think we're going to change the names of our cdf functions at this point. So it'd have to be new function names or modification/overload of existing ones.