# How to sample many ordered iid variables in turing
**URL:** https://discourse.julialang.org/t/how-to-sample-many-ordered-iid-variables-in-turing/127205
**Category:** Probabilistic Programming
**Tags:** turing, distributions, sampling
**Created:** [March 21, 2025, 7:12am UTC](https://discourse.julialang.org/t/how-to-sample-many-ordered-iid-variables-in-turing/127205 "2025-03-21T07:12:12Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@swsw](https://discourse.julialang.org/u/swsw)
#### Post date: [March 21, 2025, 7:12am UTC](https://discourse.julialang.org/t/how-to-sample-many-ordered-iid-variables-in-turing/127205/1 "2025-03-21T07:12:12Z")
Hi there,
I’d like to sample many ordered iid varialbles but i find when the number of samples is large the efficiency gets very low. For the subsequent code it works when the number of variables is 30 or so but not when it exceeds 40.
I specify that all of these variables X ~ Uniform(0,1) iid and that X[1]\ using Turing, Distributions
>
> @model function demo()
> X ~ filldist(Uniform(0, 1), n)
>
> for i in 1:(length(X)-1)
> if !(X[i] \< X[i+1])
> Turing.@addlogprob! -Inf
> return
> end
> end
>
> return X
> end
>
> n=50 # n=30 works
> n\_samples = 10000
> r = sort(rand(n)) \* 0.8
>
> @time chain = sample(demo(), NUTS(), n\_samples, burnin=n\_samples\*0.5; thinning=20, progress=true, initial\_params=r)
---
### Author:  [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [March 21, 2025, 10:29am UTC](https://discourse.julialang.org/t/how-to-sample-many-ordered-iid-variables-in-turing/127205/2 "2025-03-21T10:29:04Z")
Your `X` parameter should be a random draw from the joint distribution of all `n` order statistics from a `U(0, 1)` sample of length n. You really don’t want to sample this with `@addlogprob!`, as there are n! ways to permute your individual draws, and only 1 of them will be ordered.
The easy way to handle these is
```julia
@model function foo(n)
X ~ Bijectors.ordered(filldist(Uniform(), n))
end
```
For me this model hangs once n\>12 for some unknown reason.
In your specific case, there’s another trick that will probably be even faster. Namely, let Z\_i \sim \mathrm{Exponential}(1) for i \in 1\ldots n+1. If X\_i = \frac{\sum\_{j=1}^i Z\_i}{\sum\_{j=1}^{n+1} Z\_j}, then [X\_1, \ldots, X\_n] follows the distribution you want. `Distributions.JointOrderStatistics` uses this method for random sampling. You could use this in a model like:
```julia
@model function bar(n)
Z ~ filldist(Exponential(), n + 1)
Zsum = cumsum(Z)
X := view(Zsum, 1:n) ./ Zsum[n+1]
end
```
Here’s a demo for n=100:
```julia
julia> chns = @time sample(bar(100), NUTS(0.8), MCMCThreads(), 1_000, 4);
┌ Info: Found initial step size
└ ϵ = 1.6
┌ Info: Found initial step size
└ ϵ = 1.6
┌ Info: Found initial step size
└ ϵ = 1.6
┌ Info: Found initial step size
└ ϵ = 0.8
Sampling (4 threads) 100%|█████████████████████████████████████████████████████████████| Time: 0:00:10
10.578692 seconds (62.78 M allocations: 69.804 GiB, 40.24% gc time)
julia> x_mean = vec(mean(sort(rand(10_000, 100); dims=2); dims=1));
julia> x_mean_turing = vec(mean(stack(MCMCChains.get_params(chns; flatten=false).X); dims=(1, 2)));
julia> maximum(abs, x_mean - x_mean_turing)
julia> isapprox(x_mean, x_mean_turing; rtol=0.01)
true
```
---
### Author:  [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [March 24, 2025, 2:34pm UTC](https://discourse.julialang.org/t/how-to-sample-many-ordered-iid-variables-in-turing/127205/3 "2025-03-24T14:34:56Z")
Maybe this is what you mean by your original reference to the Dirichlet distribution, but another way you could do this is
```julia
@model function baz(n)
Z ~ Dirichlet(ones(n + 1))
X := cumsum(view(Z, 1:n))
end
```
This works because one can sample from `Dirichlet(ones(n + 1))` by taking the L1-norm of n+1 IID draws from a standard exponential (see [Dirichlet distribution - Wikipedia](https://en.wikipedia.org/wiki/Dirichlet_distribution#Random_variate_generation)). But the first n entries of the resulting cumulative sum is equivalent to what we did above in the `bar` model. Here, however, a different Bijector will be used to unconstrain the distribution, so the sampling efficiency may be better/worse.