# Sampling ordered values inside a 2-simplex

**URL:** <https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233>\
**Category:** v5\
**Tags:** modeling\
**Created:** [April 12, 2024, 11:17am UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233 "2024-04-12T11:17:10Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![sergiopasra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sergiopasra/32/563_2.png) [@sergiopasra](https://discourse.pymc.io/u/sergiopasra)\
**Post date:** [April 12, 2024, 11:17am UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233/1 "2024-04-12T11:17:10Z")

</div>

I’m trying to sample values x1 and x2 defined as:

0 \<= x1 \<= x2 \<= 1  
x1 + x2 \<= 1

The second condition can be fulfilled by doing:

```auto
x_3d = pm.Dirichlet("x_3d", [1, 1, 1])
x_2d = x_3d[:2]

```

The values in `x_2d` add up to less than or equal to 1 and are both positive. But how do I sort the values?  
Is it possible to apply the `ordered` transform to `x_2d`?

---

<div class="post-metadata">

**Author:** ![Genesis](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/genesis/32/6017_2.png) [@Genesis](https://discourse.pymc.io/u/Genesis)\
**Post date:** [April 12, 2024, 11:45am UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233/2 "2024-04-12T11:45:54Z")

</div>

Coincidentally, I am also trying to work this out. My understanding that right now this is not possible because of the bug in the ordered transform ([Ordered transform incompatible with constrained space transforms (ZeroSum, Simplex, etc.) · Issue #6975 · pymc-devs/pymc · GitHub](https://github.com/pymc-devs/pymc/issues/6975))

I am playing around with doing it “manually” but I am _very_ early in the process. All I got so far is that the following transformation can give us an ordered simplex from an unconstrained space:

```auto
import pytensor.tensor as tt

def ordered_simplex(vals):
  n = vals.shape[-1]
  v_exp = tt.exp(vals)
  weights = tt.arange(start=n, stop=0, step=-1)
  return tt.cumsum(v_exp)/tt.dot(v_exp, weights)

```

But I don’t think on its own this is enough… I am not sure right now:

1. Whether to define a PyMC Transform subclass using this (which would require also defining and inverse and a log determinant of a Jacobian) and then passing this to Dirichlet in order to override the default Simplex transform that the distribution uses
  - Also if I go this way the transform would only work for MCMC; it’s not going to produce ordered simplices using forward sampling and I am not sure yet how to do **that**

2. Whether to eschew Dirichlet entirely and just sample some, say, Normals and just pipe them directly through `ordered_simplex`. But (a) I am not sure how “advisable” that is and (b) I think that wouldn’t allow it to be used for an observed random variable, for which you’d end up needing to definte a custom distribution. And at that stange I am not sure whether it’s any different/better than overwriting the default transform for Dirichlet (but again, not sure how to deal with defining the `rv` operation for it)

Maybe @ricardoV94 has some thoughts?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [April 12, 2024, 4:44pm UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233/3 "2024-04-12T16:44:14Z")

</div>

Just to clarify, it’s not a bug, those transforms are not built to be composable. Namely order transform is a distorting transform. An ordered variable will in general not correspond to any simple distribution.

Also observed variables can’t be transformed, so I’m abut confused by your text.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [April 12, 2024, 4:46pm UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233/4 "2024-04-12T16:46:04Z")

</div>

> [@sergiopasra](#):
>
> The values in `x_2d` add up to less than or equal to 1 and are both positive. But how do I sort the values?  
> Is it possible to apply the `ordered` transform to `x_2d`?

You can apply a ChainedTransform (ordered + simplex) to x\_3d, then x\_3d[:2] will follow both properties

---

<div class="post-metadata">

**Author:** ![sergiopasra](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sergiopasra/32/563_2.png) [@sergiopasra](https://discourse.pymc.io/u/sergiopasra)\
**Post date:** [April 16, 2024, 5:05pm UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233/5 "2024-04-16T17:05:08Z")

</div>

But if you order a 3-vector and then you discard one component, the result is different to discard one component and then sort the two remaining values.

```auto
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

rng = np.random.default_rng()
n = 1000

# Sample, drop one coordinate, then sort
sx3 = scipy.stats.dirichlet.rvs([1,1,1], size=n, random_state=rng)
sx_sorted = np.sort(sx3[:,:2], axis=1)

# Sample, sort, drop one coordinate
sy3 = scipy.stats.dirichlet.rvs([1,1,1], size=n, random_state=rng)
sy_sorted = np.sort(sy3, axis=1)[:,:2]

plt.scatter(sx_sorted[:,0], sx_sorted[:,1], marker='+')
plt.scatter(sy_sorted[:,0], sy_sorted[:,1], marker='.');
plt.show()

```

![two_ways](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/6/611055bfc35ef32c76a69b77a60a1ca7e58a67a6.png)

The points marked with ‘+’ fullfill both conditions. Orange points are biased

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [April 16, 2024, 5:10pm UTC](https://discourse.pymc.io/t/sampling-ordered-values-inside-a-2-simplex/14233/6 "2024-04-16T17:10:53Z")

</div>

Yes, that makes sense. Doesn’t the ChainTransform followed by drop give you the less biased case?