Add OrderStatistic and JointOrderStatistics distributions

[sethaxen]

This PR will add two new distributions:

• OrderStatistic: a univariate distribution representing the distribution of the rank ith draw in an IID sample of length n from some (continuous or discrete) univariate distribution dist.
• JointOrderStatistics: a continuous multivariate distribution representing the joint distribution of the rank r draws in the same sample.

These have a number of uses. For example, they can represent transformations that have been performed to data (e.g. outlier removal and calculation of summary statistics). They also can be used to draw confidence intervals, e.g. of ECDF plots. For JointOrderStatistics, in the extreme case where r=1:n, i.e. when the distribution is over all ranks, it can be used in a PPL like Turing to restrict an IID array of parameters to the subset of arrays that are ordered.

Relates #1643 #1655, closes #1284

[codecov-commenter]

Codecov Report

Patch coverage: 98.11% and project coverage change: -0.03 ⚠️

Comparison is base (9cf6a74) 85.92% compared to head (2804662) 85.89%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1668      +/-   ##
==========================================
- Coverage   85.92%   85.89%   -0.03%     
==========================================
  Files         139      142       +3     
  Lines        8376     8560     +184     
==========================================
+ Hits         7197     7353     +156     
- Misses       1179     1207      +28     

Impacted Files	Coverage Δ
src/Distributions.jl	100.00% <ø> (ø)
src/multivariates.jl	44.82% <ø> (ø)
src/univariate/orderstatistic.jl	96.66% <96.66%> (ø)
src/multivariate/jointorderstatistics.jl	98.66% <98.66%> (ø)
src/univariates.jl	75.45% <100.00%> (+0.90%)	⬆️

... and 4 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[sethaxen]

I've only seen the expression of the density of joint order statistics 1:n or [i,j], so for completeness, here's the derivation for the more general case used in JointOrderStatistics:

Derivation

Given the density function $f_X$ of the sample distribution, with CDF $F_X$, the joint density for all order statistics is
$$f_{1\ldots n:n}(x_1, \ldots, x_n) = n! \prod_{k=1}^n f_X(x_k)$$

Suppose now we marginalize all $x_k$ for $i < k < j$, where $i$ and $j$ are integers in $1\ldots n$. The resulting density is
$$f_{1\ldots i,j\ldots n:n}(x_1, \ldots, x_i, x_j, \ldots x_n) = n! \left(\prod_{k=1}^i f_X(x_k)\right) \left(\prod_{k=j}^n f_X(x_k)\right) C_{ij}(x_i, x_j),$$
where
$$C_{ij}(x_i, x_j) = \int_{x_i}^{x_j} \cdots \int_{x_i}^{x_{i+3}}\int_{x_i}^{x_{i+2}} f_X(x_{i+1}) f_X(x_{i+2}) \cdots f_X(x_{j-1}) \mathrm{d}x_{i+1} \mathrm{d}x_{i+2} \cdots \mathrm{d}x_{j-1}.$$

Let $u_k = F_X(x_k) - F_X(x_i)$. Then $\mathrm{d}u_k = f_X(x_k) \mathrm{d}x_k$, and
$$C_{ij}(x_i, x_j) = \int_{0}^{u_j} \cdots \int_{0}^{u_{i+3}}\int_{0}^{u_{i+2}} \mathrm{d}u_{i+1} \mathrm{d}u_{i+2} \cdots \mathrm{d}u_{j-1} = \frac{u_j^{j-i-1}}{(j-i-1)!} = \frac{(F_X(x_j) - F_X(x_i))^{j-i-1}}{(j-i-1)!}$$

We can repeat this to get the joint density of any subset of order statistics.

[sethaxen]

I haven't added them to the docs yet. Univariate distributions are currently separated into Continuous and Discrete categories, but OrderStatistic can be either. Plus, JointOrderStatistics and OrderStatistic are so closely related, it might make sense to give them their own docs page. Could even be a page devoted to distributions of statistics of samples, in case in the future someone adds asymptotic distributions of statistics like quantiles, etc.

[sethaxen]

Probably maybe worth creating samplers for these two cases.

[sethaxen]

For the continuous case, we can write down the PDF and not much else (except the mean and covariance for a handful of dists, for a future PR).

For the discrete case, the PDF is much more complicated and probably only tractable when length(r) <= 2. For a discrete dist whose support is a finite set, we can define a distribution of the order statistics for a sample without replacement, at least in the bivariate case. The PMF, mean, and covariance are known. I wonder if that should be a special case of this distribution.

[sethaxen]

I wonder if that should be a special case of this distribution.

Nah, seems complicated and not all that useful.

[sethaxen]

Bump @devmotion

[devmotion]

Sorry for the embarassingly long delay, I had completely forgotten this PR (again).

I made a few comments. I also wonder if we could add a comparison with existing R implementations, to exploit additionallu the existing test infrastructure and run standardized tests for evaluations and sampling?

Plus, JointOrderStatistics and OrderStatistic are so closely related, it might make sense to give them their own docs page

I agree, I think a separate page in the docs could be useful.

[sethaxen]

Sorry for the embarassingly long delay, I had completely forgotten this PR (again).

No problem!

I also wonder if we could add a comparison with existing R implementations, to exploit additionallu the existing test infrastructure and run standardized tests for evaluations and sampling?

Unfortunately I don't think there are any standard R implementations of distributions of order statistics against which we can compare.

I agree, I think a separate page in the docs could be useful.

Done!

[devmotion]

Hmm, no I'm not aware of such an example. Might be good enough for now.

[devmotion]

Great PR, looks good to me @sethaxen!

I spotted only a few final things that, IMO, should be addressed before merging the PR. But I already approve now to indicate that it seems fine to me otherwise.

[sethaxen]

Great PR, looks good to me @sethaxen!

Great, thanks for the review!

[sethaxen]

I realized that when we evaluated pdf(::JointOrderStatistics, x) for unsorted x, an error would be raised by logdiffcdf, so I fixed that. This should now be ready to merge.

[devmotion]

Great! I'll merge it within the next days, if there are no objections.