incorrect formula in `bilby.core.prior.conditional.DirichletElement` docstring

[Qazalbash]

I have noticed that the mentioned formula does not match the implementation,

bilby/bilby/core/prior/conditional.py

Line 332 in 768e1e0

p(x_n) \propto (x_\max - x_n)^{(N - n - 2)}

$$ p(x_n) \propto (x_\max - x_n)^{N - n - 2} $$

where $N$ is the dimension, $n$ is the order.

The implementation is,

$$ x_n\mid x_\max\sim B(1, N-n-1) $$

where, $\mathrm{B}$ is the beta distribution with the $\mathrm{supp}(B)=[0,1-x_\max]$.

Simplifying the above formula gives,

$$ \forall n < N-1,\qquad p(x_n\mid x_\max) = \frac{(N - n - 1)(1-x_\max - x_n)^{N - n - 2}}{(1-x_\max)^{N - n - 1}},\qquad x_n\in[0,1-x_\max] $$

and for $n=N-1$,

$$ x_{N-1}\mid x_\max \sim \mathrm{Unif}(0,1-x_\max) $$

---

If you find my claim true, I will be happy to make a PR!

[ColmTalbot]

I think the issue is elsewhere in the documentation, and at least a lack of clarity. I think we should say $i &lt; n - 1$ as the (n - 1)th component is not sampled as that is a delta function to ensure closure of the Dirichlet distribution. You can check this with $i=1$ pretty straightforwardly. I think you're also misremembering the definition of the beta distribution.

$$ p(x; \alpha, \beta) = \frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)} $$

if we set $\alpha = 1$, $\beta = N - n - 2$ we get

$$ p(x; 1, N - n - 2) = \frac{(1 - x)^{N - n - 2}}{B(1, N - n -2)} $$

[Qazalbash]

As per the implementation, $\beta=N-n-1$

bilby/bilby/core/prior/conditional.py

Lines 363 to 367 in 768e1e0

	super(DirichletElement, self).__init__(
	minimum=0, maximum=1, alpha=1, beta=n_dimensions - order - 1,
	name=label + str(order),
	condition_func=self.dirichlet_condition
	)

---

if we set α = 1 , β = N − n − 2 we get

p ( x ; 1 , N − n − 2 ) = ( 1 − x ) N − n − 2 B ( 1 , N − n − 2 )

I think you also made a mistake here, the correct beta distribution would be,

$$ p ( x ; 1 , N − n − 2 ) = \frac{( 1 − x ) ^ { N − n − 3} }{B ( 1 , N − n − 2 )} $$

And how can be the $\mathrm{supp}(p)=[0,x_\max]$, when the code is setting the maximum to $1-x_\max$,

bilby/bilby/core/prior/conditional.py

Lines 377 to 381 in 768e1e0

	def dirichlet_condition(self, reference_parms, **kwargs):
	remaining = 1 - sum(
	[kwargs[self.label + str(ii)] for ii in range(self.order)]
	)
	return dict(minimum=reference_parms["minimum"], maximum=remaining)

And from the implementation POV, if $\mathrm{supp}(p)=[0,1-x_\max]$, then $p(x_n\mid x_\max)\propto (1-x_\max-x_n)^{N-n-2}$

---

I am talking from the POV of implementation and the description in the docstring. I have checked my derivation multiple times. If you like, I can share here too.

I have implemented this using NumPyro; you can compare with that too.

[ColmTalbot]

Okay, I think I'm caught up now. I missed that the issue is that $x_\max \rightarrow 1 - x_\max$ in the parentheses. I agree that this should be changed. I'm also not a fan of using $x_\max$ here now that I'm looking at it again. We could use something like $X_n$ instead of $x_\max$. The N = 2, n=1 example should also change to 0 < x_1 < 1 - x_0.

---

I did have a typo, but I think not the way you changed it to. My original example should have been if we set $\alpha = 1$, $\beta = N - n - 1$ (as it is in the code) we get

$$ p(x; 1, N - n - 2) = \frac{(1 - x)^{N - n - 2}}{B(1, N - n - 1)} $$

[Qazalbash]

I have derived the exact formula; you can use that as the docstring if you like. And I agree $x_\max$ is a confusing notation; it should be changed.

[ColmTalbot]

@Qazalbash if you're willing to add this fix, please assign the issue to yourself.

[Qazalbash]

Sure, I will make a PR shortly. Do you have any suggestions on $x_\max$?

[ColmTalbot]

I would probably go for $X_i$, but that isn't a strong preference.