Open sethaxen opened 2 years ago
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.
OrderedBijector
is only appropriate for bounded real numbers. It would be convenient to also have aPositiveOrderedBijector
(similar to Stan'spositive_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
andOrderedBijector
in reverse. As a bonus, ifordered
is passed a distribution with a defined positive bijector, then it could instead returnPositiveOrderedBijector
to avoid incorrect posterior draws.