MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top


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]$


$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.


share|cite|improve this question

You can randomly sample from the full simplex (classically, one can do this by choosing iid exponential random variables with mean one and normalizing).

If you then order the resulting elements of the simplex, you get a (uniformly chosen) element of the 'non-decreasing' set that you're interested in. So these two things let you do 'perfect simulation', and thus approximate any function.

If you're only interested in the volume, the above suggests that the answer should just be the volume of the simplex over n!

This, or I'm misreading the question... even if I am (which seems likely), this is certainly the type of question for which MCMC methods are often useful. This should let you get rid of the worst of the rejection sampling problems. If you feel like being fancy, randomness recyclers also work for questions of this type... but unlike MCMC there don't seem to be any easy/general writeups of that method.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.