Statistical distribution arithmetic

[Happypig375]

Companion to #1376. An approximate value (such as a variable with measurement error) should really be interval arithmetic. However, most measurements in life come with a distribution. For example as normal distribution. Interval arithmetic necessarily assumes a uniform distribution by asserting that the variable can be simply described using its minimum and maximum values.

The idea here is to figure out how we can process such distributions. For example, the summation of two normally distributed variables are still normally distributed while the summation of the square of independent standard normal variables distribute as chi-squared. Concepts like dependence also need to be exposed.

Whether we are able to process these distributions in integral form properly for all functions is part of this investigation. Where convenient, a simplified node representing a well known distribution should be used instead. Therefore, this must only start after integration work is done.

[Rafael-SOWNet]

Two notes for the investigation, nothing implemented (the Proposal stands, and the gate you state — after the integration work — is the right one, since the general case is integration).

Intervals and distributions are different claims, not a special case of each other. An interval asserts the support — the value is in [lo, hi] — and says nothing about how likely each point is; it is the distribution-free statement, sound for every distribution with that support, which is why it is what measurement error is usually propagated as when the distribution is unknown. A uniform distribution on the same interval is a much stronger claim (it fixes the density), and a sum of two intervals is again an interval whereas a sum of two uniforms is triangular. So the two belong to different layers of the type design in #1237 / #1019 item 10: an interval-valued number is a number type with an order-aware arithmetic; a distribution-valued variable is a variable whose declaration carries a law, and its arithmetic is the arithmetic of laws.

What the arithmetic of laws needs from the library, and what it already has. The closed families are tables — a linear combination of independent normals is normal, a sum of squares of independent standard normals is χ², a sum of independent Poissons is Poisson, products and quotients of log-normals are log-normal — and those are pattern rules over a Distribution node with typed parameters, cheap and exact. Outside the tables there are three mechanisms, each already partly in the library: for Y = f(X) with f monotone the density is p_X(f⁻¹(y)) |d f⁻¹/dy|, which needs inverse functions and differentiation (both here) and, for a non-monotone f, a sum over the preimages — the root set of f(x) = y, which the solver gives; for X + Y independent it is the convolution ∫ p_X(t) p_Y(y − t) dt, an integral of a product of densities, which is #718's work and where most of the difficulty sits (piecewise supports, Providedf conditions); and the moment/characteristic-function route, which turns convolution into multiplication and needs Fourier/Laplace transforms the library does not have. Dependence is the part with no shortcut: it has to be declared (a joint law, or a copula, or a functional relation Y = g(X)), because two variables with the same marginals and different joint laws give different sums — which is also the argument for the declaration living on the variable, as in #1237, rather than being inferred.

A sketch of the order, for when the time comes: (1) the node and its parameters, with the closed tables and independence assumed unless declared otherwise; (2) monotone transforms through the inverse and the derivative; (3) convolution through the integrator, with the density's support as a Providedf so the piecewise machinery carries it; (4) declared dependence. Each step is checkable numerically by sampling, which is the harness I would write first.

[Happypig375]

Thanks for clarifying on intervals vs statistical distributions. Aside from continuous distributions modelled as integrals, there is also the discrete distribution from a known population size. I don't think discrete distributions are modelled as integrals (integration deals with the continuous, right? unless we want to add the Dirac delta function which isn't really a function? maybe we can expand it to a limit instead...), but they are also not sets because sets have implicit deduplication. We need the probability mass function - what happens if a probability mass function is combined with a probability density function?