# How to sample ordered variables with different distributions in turing

**URL:** <https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455>\
**Category:** Probabilistic Programming\
**Tags:** question, package, turing\
**Created:** [January 7, 2024, 4:55am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455 "2024-01-07T04:55:58Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![gonzalo](https://avatars.discourse-cdn.com/v4/letter/g/4bbf92/32.png) [@gonzalo](https://discourse.julialang.org/u/gonzalo)\
**Post date:** [January 7, 2024, 4:55am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/1 "2024-01-07T04:55:58Z")

</div>

Hi julia bros,

I want to use turing to sample ordered variables each with different distributions. For example,

> X~Normal(0,1) Y~Exp(1), Z~Uniform(0,1)

and require

> X\>Y\>Z

I know this can be done by rejection sampling but in case there are lots of ordered variables the efficiency is very low. Is there a way in turing or other packages in julia to help?

Gonzalo

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 7, 2024, 12:27pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/2 "2024-01-07T12:27:15Z")

</div>

Would try something like

```julia
@model function demo()
    Z ~ Uniform(0, 1)
    Y ~ truncated(Exponential(1); lower = Z)
    X ~ truncated(Normal(); lower = Y)
    [X, Y, Z]
end

chain = sample(demo(), NUTS(), 1_000; progress = true)

```

Not sure if it works correctly though, as it returns some samples violating the constraints … maybe a bug?

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 7, 2024, 1:33pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/3 "2024-01-07T13:33:10Z")

</div>

```julia
chain = sample(demo(), SMC(), 1_000; progress = true);

```

Seems to return non-violating samples.

---

<div class="post-metadata">

**Author:** ![gonzalo](https://avatars.discourse-cdn.com/v4/letter/g/4bbf92/32.png) [@gonzalo](https://discourse.julialang.org/u/gonzalo)\
**Post date:** [January 8, 2024, 2:25am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/4 "2024-01-08T02:25:47Z")

</div>

Thank you. In your code Z ~ Uniform(0,1) and the others are conditioned on its value but that’s not accurate because Z’s distribution is not affected. A rejection sampling can show that X,Y,Z’s distributions are not what they are specified at the beginning.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 8, 2024, 6:56pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/5 "2024-01-08T18:56:58Z")

</div>

Well, yes, that was my interpretation of your definition.  
Seems you want a different density though. In this case, you will probably need to define a transformation from \mathbb{R}^3 \to (x, y, z) where x \in [0, 1], y \> x, z \> y and then use `Turing.@addlogprob!` to explicitly define the desired density – including the log-Jacobian of the transformation – or define a custom `Bijector` and use this to implicitly transform the desired density. Have never done this in Turing though.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 8, 2024, 7:46pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/6 "2024-01-08T19:46:17Z")

</div>

Following up, does this do what you want?

```julia
struct CustomDistribution <: ContinuousMultivariateDistribution
end

Base.length(c::CustomDistribution) = 3

Bijectors.bijector(c::CustomDistribution) = Bijectors.Stacked([Bijectors.Logit(0, 1), identity, identity]) ∘ inverse(Bijectors.OrderedBijector())

Distributions.logpdf(c::CustomDistribution, v::AbstractVector) = sum(logpdf.([Uniform(0, 1), Exponential(1), Normal(0, 1)], v))

@model function demo()
    xyz ~ CustomDistribution()
end

chain = sample(demo(), NUTS(), 1_000; progress = true)

```

Note: This is just hacked together as MWE. Please check the docs of `Distributions` and `Bijectors` in more detail on how to define and use custom distributions!

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 8, 2024, 10:15pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/7 "2024-01-08T22:15:14Z")

</div>

> [@bertschi](#):
>
> `Bijectors.Logit(0, 1)`

Is there a bijector which simply wraps the reals mod 1 ? It isn’t 1-1 but it is continuous and a local bijection.  
Perhaps I should go and look at the Bijectors code, but here I am asking this question.

(as for the relevance, using such a uniform bijector will generate different sampling distribution, which may somehow correspond better to some desired distribution).

UPDATE: There is a problem though with using such a bijector, as it doesn’t conserve ordering (it is non-monotonic) and wouldn’t work in combination with the second bijector as suggested.

---

<div class="post-metadata">

**Author:** ![gonzalo](https://avatars.discourse-cdn.com/v4/letter/g/4bbf92/32.png) [@gonzalo](https://discourse.julialang.org/u/gonzalo)\
**Post date:** [January 9, 2024, 4:46am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/8 "2024-01-09T04:46:05Z")

</div>

Thank you so much! That appears nicely solve my questions. I am not familiar with Bijectors but what you present gives the same result to rejection sampling.

However if the distributions for X,Y,Z are changed, for example X,Y,Z iid ~ Exponential(1), then I use your code and got mean(X) = 0.28 which however should be 1/3 (can be shown analytically or with rejection sampling).

> function rejection\_sampling\_exp3(n\_samples::Int)  
> samples\_X = Float64  
> samples\_Y = Float64  
> samples\_Z = Float64
> 
> ```
> while length(samples_X) < n_samples
> x = rand(Exponential(1))
> y = rand(Exponential(1))
> z = rand(Exponential(1))
>     
> if x < y < z
> push!(samples_X, x)
> push!(samples_Y, y)
> push!(samples_Z, z)
> end
> end
> 
> return samples_X, samples_Y, samples_Z
> 
> ```
> 
> end
> 
> n\_samples = 100000
> 
> samples\_X, samples\_Y, samples\_Z = rejection\_sampling\_exp3(n\_samples)
> 
> mean\_X = mean(samples\_X)  
> mean\_Y = mean(samples\_Y)  
> mean\_Z = mean(samples\_Z)
> 
> println(“Mean X: $mean\_X”)  
> println(“Mean Y: $mean\_Y”)  
> println(“Mean Z: $mean\_Z”)

I guess this is related to transforming exponentially distributed X~Exp(1) by the Bijectors.

Is there anything that I got wrong? Thanks so much!

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 9, 2024, 6:01am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/9 "2024-01-09T06:01:17Z")

</div>

This following does rejection sampling, but does it with Turing:

```julia
using Turing

@model function demo()
    Z ~ Uniform(0, 1)
    Y ~ Exponential(1)
    X ~ Normal()
    if !(Z < Y < X)
        DynamicPPL.@addlogprob! -Inf
        return
    end
    [X, Y, Z]
end

chain = sample(demo(), NUTS(), 1_000; progress = true)

```

The average of Z in my test came out to be 0.338

---

<div class="post-metadata">

**Author:** ![gonzalo](https://avatars.discourse-cdn.com/v4/letter/g/4bbf92/32.png) [@gonzalo](https://discourse.julialang.org/u/gonzalo)\
**Post date:** [January 9, 2024, 9:04am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/11 "2024-01-09T09:04:41Z")

</div>

Thanks you much! working pretty well in my exampke. That’s very super clever way.

---

<div class="post-metadata">

**Author:** ![gonzalo](https://avatars.discourse-cdn.com/v4/letter/g/4bbf92/32.png) [@gonzalo](https://discourse.julialang.org/u/gonzalo)\
**Post date:** [January 9, 2024, 10:34am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/12 "2024-01-09T10:34:17Z")

</div>

By the way, do you know how with VI this can be performed?

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 9, 2024, 11:40am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/13 "2024-01-09T11:40:30Z")

</div>

> [@Dan](#):
>
> Is there a bijector which simply wraps the reals mod 1 ? It isn’t 1-1 but it is continuous and a local bijection.

As you mention yourself, it would not be 1-1, i.e., not a bijector. This will be a problem when trying to compute the Jacobian correction it has on the density – which is needed in your case as you want to define a density on the constrained, i.e., transformed, space.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 9, 2024, 12:04pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/14 "2024-01-09T12:04:37Z")

</div>

> [@gonzalo](#):
>
> Thank you so much! That appears nicely solve my questions. I am not familiar with Bijectors but what you present gives the same result to rejection sampling.
> 
> However if the distributions for X,Y,Z are changed, for example X,Y,Z iid ~ Exponential(1), then I use your code and got mean(X) = 0.28 which however should be 1/3 (can be shown analytically or with rejection sampling).

In this case, you would need to change the bijector to match the support of your distribution:

- `Logit(0, 1)`: Transforms from a constrained space x \in [0, 1], i.e., the support of uniform, to an unconstrained value y \in \mathbb{R}:

```julia-repl
julia> using Bijectors
with_logabsdet_jacobian(Bijectors.Logit(0, 1), 0.1) # (unconstrained value, Jacobian correction)

```

Turing uses this to define a density p(y) on the unconstrained space which nevertheless has the desired density p(x) when mapped onto the constrained space. By change of variable, this requires to adjust the density by the log-abs-Jacobian as p(y) = p(x) |dx/dy| where dx/dy denotes the Jacobian of the inverse transformation. The [Stan manual](https://mc-stan.org/docs/stan-users-guide/changes-of-variables.html) has a nice explanation on this, but also the source code of `Bijectors` is rather clear, i.e., just try `@less logpdf(transformed(Uniform(0, 1), Bijectors.Logit(0, 1)), 0.2)`

Long story short, in case of the `Exponential` distribution, you need a bijector from \mathbb{R}^+ \to \mathbb{R} instead, i.e., just replace `Bijectors.Logit(0, 1)`  
by `bijector(Exponential(1))` – should probably have used `bijector(Uniform(0, 1))` in the first place.

Regarding VB: This should work out-of-the-box for the model using `CustomDistribution`.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 9, 2024, 11:42pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/15 "2024-01-09T23:42:44Z")

</div>

> [@bertschi](#):
>
> This will be a problem when trying to compute the Jacobian correction it has on the density

This is a bit of a tangent to the thread. But I think even though it is not a bijector, using this transformation is valid, and might be useful. The two caveats needed are:  
(i) Continuity of density/transform on the stiched domains of the covering space - without which some algos might get confused (although this may not be a strict requirement)  
(ii) A ‘uniform’ covering. For example the mod1 cover of [0,1) with R has the required symmetry (the covering can be partitioned into branches)  
With these conditions, the stable distribution on the original space should be maintained (and the sampling will work as intended).

Such more general transformation would be a nice addition (and maybe Stan doesn’t have them).

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 10, 2024, 7:34pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/16 "2024-01-10T19:34:06Z")

</div>

Still not sure that something like `mod 1` is a good idea: First, the posterior on the unconstrained space will be multi-modal as many points map onto the same constrained parameter and accordingly also have the same likelihood. Secondly, and even worse, this means that the posterior will not be proper as it cannot be integrated – without any additional density defined on the unconstrained space directly.  
Why do you think such an approach would be better than a 1-1 mapping, e.g., via inverse logit?

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 10, 2024, 7:51pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/17 "2024-01-10T19:51:04Z")

</div>

> [@bertschi](#):
>
> Why do you think such an approach would be better than a 1-1 mapping, e.g., via inverse logit?

I have an intuition it might converge faster.  
Speed is not so much of an issue here, but perhaps a different problem can amplify the speed difference.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 10, 2024, 7:55pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/18 "2024-01-10T19:55:08Z")

</div>

Would be surprised by that – and in any case, it would only work as most practical samplers are not able to explore the true posterior which is improper and multi-modal.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 10, 2024, 8:01pm UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/19 "2024-01-10T20:01:14Z")

</div>

> [@bertschi](#):
>
> improper and multi-modal

I don’t understand this comment fully:  
Samplers are supposed to sample the distribution which is often multi-modal.  
Samplers are not integrating, just sampling. They don’t really know about a global property of the distribution, like whether it is proper.  
But I’m not really up to date with all the modern sampling tricks.  
This is off-thread topic. So maybe we should let this idea slip through the cracks unless someone picks it up (maybe I will).

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [May 29, 2024, 8:06am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/20 "2024-05-29T08:06:07Z")

</div>

> [@Dan](#):
>
> I don’t understand this comment fully:  
> Samplers are supposed to sample the distribution which is often multi-modal.  
> Samplers are not integrating, just sampling. They don’t really know about a global property of the distribution, like whether it is proper.

This is strictly not correct. MCMC algorithms _are_ numerical integration methods. If certain conditions are satisfied, the resulting samples can be used to approximate certain expectations numerically. Moreover, when those conditions are satisfied, there’s a central limit theorem that bounds the error in the approximation of expectations. MCMC methods _cannot_ sample from an improper distribution. “Improper distributions” are not distributions at all. One could never draw exact samples from an “improper distribution,” so it is also not possible to draw samples using MCMC methods. Multi-modal distributions are in principle no problem for MCMC, but in practice most widely used MCMC algorithms would require infinite time to jump between widely-separated modes. So when possible, it is a good idea to parametrize a model in such a way that it is unimodal.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [May 29, 2024, 8:26am UTC](https://discourse.julialang.org/t/how-to-sample-ordered-variables-with-different-distributions-in-turing/108455/21 "2024-05-29T08:26:40Z")

</div>

> [@sethaxen](#):
>
> This is strictly not correct.

Few things are.

According to [Markov chain Monte Carlo - Wikipedia](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) :  
" In statistics, **Markov chain Monte Carlo** ( **MCMC** ) is a class of algorithms used to draw samples from a probability distribution. " but it is true that:  
" MCMC methods are primarily used for calculating numerical approximations ".

> [@sethaxen](#):
>
> One could never draw exact samples from an “improper distribution,”

You can’t even draw from the Normal Distribution on a computer (everything is finite…).

> [@sethaxen](#):
>
> but in practice most widely used MCMC algorithms would require infinite

Also inaccurate, and depends mostly on the width of the chain’s “bottlenecks” (see Mixing Time literature).