6
$\begingroup$

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.

$\endgroup$
7
$\begingroup$

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):

enter image description here

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

$\endgroup$
  • $\begingroup$ How could you modify these approach if you wanted to condition on S greater than a value, i.e., generate X_1,\ldots,X_n conditional on S > 10, for example? Would you need to first inverse transform sample from S|S >10? $\endgroup$ – user21359 Apr 25 '18 at 13:48
  • $\begingroup$ @user21359: when conditioning on a positive probability event such as $S>10$, it is not possible to use a change of variable but one must constrain the original space for this event to happen. Check the keyword "cross-entropy" for efficient simulation in such cases. $\endgroup$ – Xi'an Apr 25 '18 at 14:07
  • $\begingroup$ Thanks. Do you have any examples of something similar? A search of cross entropy most yielded information about estimating the probability of rare events rather than sampling from distributions. $\endgroup$ – user21359 Apr 25 '18 at 16:40
  • $\begingroup$ I mean is there an example where someone uses cross entropy to estimate a similar problem kind of problem? I think I’m too new to this concept to construct it myself. $\endgroup$ – user21359 Apr 25 '18 at 18:36
0
$\begingroup$

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$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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