Sample random variables conditional on their sum

Question

Let $(X_1, \dots, X_n)$ be an iid sample of random variables with a known continuous distribution. I would like to simulate such a sample, conditional on the value of its sum, that is: $$ X_1, \dots, X_n | \sum_{k=1}^n X_k = s $$

If I was conditioning on, say $\sum_{k=1}^n X_k \in [s - \varepsilon, s+\varepsilon]$, then a solution would be to simulate an iid sample and discard the cases where the sum does not fall in $[s - \varepsilon, s+\varepsilon]$. It wouldn't be very efficient if $\varepsilon$ was small, but it would work. The problem here is that I'm trying to condition on a zero probability event, so that logic is not applicable.

I know there are cases (normal distribution, at least), where everything can be calculated analytically, but I seek a more general solution here.

Answer 1

If you seek the conditional density of $(X_1,...,X_{n-1})$ given $$S=\sum_{k=1}^n X_k$$ a change of variable from $$(X_1,...,X_{n})\sim\prod_{i=1}^n f(x_i)$$ to $$\left(X_1,...,X_{n-1},S\right)\sim\prod_{i=1}^{n-1}f(x_i)\times f(s-x_1-\cdots-x_{n-1})$$ [with Jacobian equal to 1] shows that this conditional density is proportional to$$f(x_1)\cdots f(x_{n-1})\,f(s-x_1-\cdots-x_{n-1})$$

Therefore there exists a closed form expression for the conditional density and one can thus call a generic simulation method to simulate from it, like accept-reject, Gibbs sampling, or a Metropolis-Hastings algorithm.

The resolution even extends to independent variables that are not identically distributed.

Note: A similar question was asked a while ago, but none of the answers mentions this generic solution.

For instance, if $f$ is the N$(0,1)$ density and $n=4$, a Metropolis-within-Gibbs sampler for this problem would be of the form

T=1e3 #Gibbs steps
n=3 #n-1
s=3.1415 #imposed sum
x=matrix(rnorm(n),T,n)
for (t in 2:T){
  x[t,]=x[t-1,]
  for (i in 1:n){
   prop=rnorm(1,x[t-1,i],3)
   if (runif(1)<dnorm(prop)*
    dnorm(s-sum(x[t,-i])-prop)/
    dnorm(x[t-1,i])/dnorm(s-sum(x[t,]))) 
     x[t,i]=prop}}

Here is the outcome of the simulation of the three (first) components $x_1$ (brown), $x_2$ (red), and $x_3$ (yellow):

[reproduced from my blog] I recently came upon an unexpected property shown by Lindqvist and Taraldsen (Biometrika, 2005) that to simulate a sample ${\bf y}$ conditional on the realisation of a sufficient statistic, $T({\bf y})=t⁰$, it is sufficient (!!!) to simulate the components of ${\bf y}$ as ${\bf y}=G({\bf u},θ)$, with ${\bf u}$ a random variable with fixed distribution, e.g., a $U(0,1)$, and to solve in $θ$ the fixed point equation $$T({\bf y})=T\circ G({\bf u},θ)=t⁰$$ assuming there exists a single solution to this equation.

To borrow a simple example from the authors, take an exponential sample ${\bf y}$ to be simulated given the sum statistic being fixed. As it is well-known, the conditional distribution of ${\bf y}$ is then a (rescaled) Beta and the proposed algorithm ends up being a standard Beta generator. For the method to work in general, $T({\bf y})$ must factorise through a function of the ${\bf u}$’s, a so-called pivotal condition. If this condition does not hold, it gets more complicated: the authors introduce a pseudo-prior distribution on the parameter $θ$ to make it independent from the ${\bf u}$’s conditional on $T({\bf y})=t⁰$. While the setting is necessarily one of exponential families and of sufficient conditioning statistics, I find it amazing that this property is not more well-known.

Answer 2

I think it's worth noting the normal case even though the questioner seeks a more general answer. Following the logic shown here How to generate two groups of $n$ random numbers in $U(0,1)$ such that sum of these two groups equal? leads to the result that for the normal case the conditional distributions are all normal.

So, for example, if we want to generate realizations of $X_1, X_2, \ldots X_n$ that sum to $z$ where the $X_i$ are all identically normal, here is the resulting approach:

(1) Generate $X_1 \sim N \left({z \over n},{\frac{n-1}{n}} \sigma^2 \right)$

(2) Generate $X_2 \sim N \left({z-x_1 \over n-1},{\frac{n-2}{n-1}} \sigma^2 \right)$

For general $i,$

(3) Generate $X_i \sim N \left({z-\sum_{j=1}^{i-1}x_j \over {n-i+1}},{\frac{n-i}{n-i+1}} \sigma^2 \right)$

(4) Generate $X_{n-1} \sim \left({z-\sum_{j=1}^{n-2}x_j \over {2}},{\frac{\sigma^2}{2}} \right)$

(5) Let $X_n = z - \sum_{j=1}^{n-1} x_j$