Conditioned Random Poisson Variables in Python

Question

In Python I would like to generate a set $X_1,..,X_n$ of random variables whose distribution is that of $n$ i.i.d. Poisson random variables with mean $1$ conditioned on

$$\sum_{i=1}^nX_i = n-1 $$ and $$1 + \sum_{j=1}^{t} \left( X_j -1\right) > 0 \text{ for all } t \in\{1,..,n-1 \}$$ At the moment I keep generating a set, check the conditions and throw it away until I find one that does satisfy these conditions. Is there any way to immediately generate a set that satisfies (at least one of) these conditions?

Answer 1

The following code only ensures the constraint set (1) that requires $1 + \sum_{j=1}^{t} \left( X_j -1\right) > 0$, for all $t \in\{1,..,n-1 \}$,

i.e., $\sum_{j=1}^{t} X_j > t - 1 \text{ for all } t \in\{1,..,n-1 \}$.

import numpy as np
import matplotlib.pylab as plt

np.random.seed(111)

n = 1000
s = 0
X = []
for i in range(1, n + 1):
    X_i = np.random.poisson(1, 1)[0]
    while s + X_i <= i - 1: # constraint set (1)
        X_i = np.random.poisson(1, 1)[0]
    s += X_i
    X.append(X_i)
plt.hist(X)
plt.show()

Now, in order to ensure the constraint (2) $\sum_{i=1}^nX_i = n-1$ additionally (along with (1) ), we can try the following:

np.random.seed(1)

n = 1000
X = np.random.poisson(1, n)
while (not all(np.cumsum(X[:-1]) > range(n-1))) or (not np.sum(X) == n - 1):
    X = np.random.poisson(1,n)

all(np.cumsum(X[:-1]) > range(n - 1)) # constraint set (1)
# True
np.sum(X) == n - 1 # constraint (2)
# True

plt.hist(X)
plt.show()

On average, the second code fragment will be much slower that the first one, since it's hard to satisfy the last constraint along with the first constraint set.