Bijectors.ordered(d) is incorrect if bijector(d) is not Identity

[sethaxen]

OrderedBijector is only appropriate for bounded real numbers. It would be convenient to also have a PositiveOrderedBijector (similar to Stan's positive_ordered vector type), which wraps a distribution required to have support ℝ_{>0}^D.

Since the exponential is monotonic, the ordering can happen in the latent space, so the bijector is just a composition of Exp and OrderedBijector in reverse. As a bonus, if ordered is passed a distribution with a defined positive bijector, then it could instead return PositiveOrderedBijector to avoid incorrect posterior draws.

[sethaxen]

More generally, bijector(d) for a given continuous univariate distribution d is guaranteed to be strictly monotone, so we have 2 cases: the bijector is increasing or decreasing.

Suppose we have an increasing distribution Gamma(1, 1) and a decreasing distribution -1 * Gamma(1, 1):

julia> dinc = Product(fill(Gamma(1, 1), 10));

julia> ddec = Product(fill(-1 * Gamma(1, 1), 10));

julia> inv(bijector(dinc))(-2:1:2)
5-element Vector{Float64}:
 0.1353352832366127
 0.36787944117144233
 1.0
 2.718281828459045
 7.38905609893065

julia> inv(bijector(ddec))(-2:1:2)
5-element Vector{Float64}:
 -0.1353352832366127
 -0.36787944117144233
 -1.0
 -2.718281828459045
 -7.38905609893065

In both cases we can construct an ordered vector of IID elements from these distributions using OrderedBijector:

julia> binc = inv(Bijectors.OrderedBijector()) ∘ bijector(dinc)
Composed{Tuple{Bijectors.TruncatedBijector{1, Float64, Float64}, Inverse{Bijectors.OrderedBijector, 1}}, 1}((Bijectors.TruncatedBijector{1, Float64, Float64}(0.0, Inf), Inverse{Bijectors.OrderedBijector, 1}(Bijectors.OrderedBijector())))

julia> binc(sort(rand(dinc)))
10-element Vector{Float64}:
 -2.022615883470717
 -1.1840408057206144
 -1.1534482444313512
 -2.213221747442945
 -1.7515607559882922
 -0.7470736409590497
 -0.08155306617081076
 -1.9816915175623477
 -1.140842898952562
 -3.1439472077666215

julia> inv(binc)(randn(10))
10-element Vector{Float64}:
      6.378410409802754
     21.967437824138944
     30.37754781370191
    131.00375891230624
    363.5112963770221
   1521.7038838401243
  16676.78769661897
  53626.234403927
  67313.59427742827
 450855.07619026705

julia> bdec = inv(Bijectors.OrderedBijector()) ∘ Bijectors.Permute(10:-1:1) ∘ bijector(ddec)
Composed{Tuple{Bijectors.TruncatedBijector{1, Vector{Float64}, Vector{Float64}}, Bijectors.Permute{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Inverse{Bijectors.OrderedBijector, 1}}, 1}((Bijectors.TruncatedBijector{1, Vector{Float64}, Vector{Float64}}([-Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Bijectors.Permute{SparseArrays.SparseMatrixCSC{Float64, Int64}}(sparse([10, 9, 8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 10, 10)), Inverse{Bijectors.OrderedBijector, 1}(Bijectors.OrderedBijector())))

julia> bdec(sort(rand(ddec)))
10-element Vector{Float64}:
 -3.6307491577258975
  0.25197946987880476
  0.665987749920476
 -1.350510882468354
 -3.7174559742545403
 -5.182266331047505
 -1.8555492211027504
 -1.2402074771788663
 -1.3204217432748224
 -0.7245043375034853

julia> inv(bdec)(randn(10))
10-element Vector{Float64}:
 -70855.84745924114
 -10452.339583021583
  -7756.019025544983
  -3688.2810163762733
   -723.4524977291115
    -12.08537333219846
    -11.016823194748373
     -6.350400824704886
     -3.1599212688906304
     -1.9292043847246376

Note that currently if we use Bijectors.ordered(dinc), the bijector of the resulting distribution evaluates the bijectors in the reverse order and so is incorrect:

julia> binc
Composed{Tuple{Bijectors.TruncatedBijector{1, Float64, Float64}, Inverse{Bijectors.OrderedBijector, 1}}, 1}((Bijectors.TruncatedBijector{1, Float64, Float64}(0.0, Inf), Inverse{Bijectors.OrderedBijector, 1}(Bijectors.OrderedBijector())))

julia> binc_wrong = bijector(Bijectors.ordered(dinc))
Composed{Tuple{Inverse{Bijectors.OrderedBijector, 1}, Bijectors.TruncatedBijector{1, Float64, Float64}}, 1}((Inverse{Bijectors.OrderedBijector, 1}(Bijectors.OrderedBijector()), Bijectors.TruncatedBijector{1, Float64, Float64}(0.0, Inf)))

julia> binc_wrong(sort(rand(dinc)))
10-element Vector{Float64}:
  -1.3132541534180986
 -Inf
 -Inf
 -Inf
 -Inf
 -Inf
 -Inf
 -Inf
  -2.7515816917132407
  -0.32023364008981076

ordered(d) should probably raise an error if bijector(d) is not an Identity. To do anything better here, we would need a distribution like Power(d::UnivariateDistribution, n::Int), for which we could specialize ordered to do the right thing.