Hi,
I am trying to compute a distribution by integrating over all non-increasing categorical distributions of a given size $n$.
For instance, for $n=2$, each categorical distribution must follow $u_0 \geq u_1 \geq 0$ and $u_0 + u_1 = 1$. For $n=3$, it must follow $u_0 \geq u_1 \geq u_2 \geq 0$ and $u_0 + u_1 + u_2 = 1$.
For a given $n$, I want to integrate over all of the categorical distributions that are consistent with these constraints.
For example, if $n=2$, then I want to compute $p = \int_{U:~u_0 \geq u_1, u_0 + u_1 = 1} U \partial U$.
I was able to simplify this a bit by transforming:
$u_0 \in [\frac{1}{n}, 1]$
$u_1 \in [\frac{1-u_0}{n}, u_1]$
$u_2 \in [\frac{1-u_0-u_1}{n}, u_2]$
$\vdots$
$u_{n-1} = 1 - u_0 - u_1 - \cdots - u_{n-2}$
Using this, I can compute what I want using nested integrals:
$p_0 = \int_{1/2}^1 u_0 \partial u_0 = 3/8$
$p_1 = \int_{1/2}^1 1-u_0 \partial u_0 = 1/8$
If $n=3$, you'd compute
$p_0 = \int_{1/3}^1 \int_{(1-u_0)/2}^{u_0} u_0 \partial u_1 \partial u_0$
$p_1 = \int_{1/3}^1 \int_{(1-u_0)/2}^{u_0} u_1 \partial u_1 \partial u_0$
$p_2 = \int_{1/3}^1 \int_{(1-u_0)/2}^{u_0} 1-u_0-u_1 \partial u_1 \partial u_0$
I want to compute a handful of the $p_k$ $p_{k_1}, p{k_2}, \cdots $ for large $n$ (on the order of a few thousand). To be exact, all I really need are values proportional to these $p_k$.
I've managed to make an algorithm that can do this integration and compute a single $p_k$ in $O(n^3)$, but the constant isn't very fast, and it gets quite slow for large $n$.
This really feels like it must be a classic problem of some sort, but I just can't seem to see or find a way to get the closed form solution or a faster algorithm.
Also, I only need to do this numerically, so sampling is an option. Unfortunately, I could only work out a rejection sampling strategy (you throw away the entire configuration U as soon as one variable violates an inequality), so the number of iterations to generate one sample becomes roughly exponential in $n$. If you know of a way that I can sample nondecreasing uniform values without rejection, I'd love to know (I also tried generating $n$ uniform values and sorting them, but I don't know how to force them to sum to 1 and still keep the same distribution; scaling didn't seem to give me the right answer).
If this looks familiar to anything you've seen, I would really, really appreciate the help. I've been stuck on this for a few days, and I'm sure you can appreciate how that makes your brain itch like there's a little splinter you can't quite get out.
Oliver