Given N independently distributed poisson variables $X_i\sim Pois(\lambda_i), i=1...N$, whats the average number of nonzero $X_i$ given $\sum X_i=M$? Given $N$ independently distributed poisson variables $X_i\sim Pois(\lambda_i), i=1...N$, whats the average number of nonzero $X_i$ (denoted as $L$) given $\sum_{i=1}^N X_i = M$?
For example, when $N=2$ and $M=2$, the result is:
$$
E(L)=\frac{\frac{1}{2!}(\lambda_1^2+\lambda_2^2)\exp(-\lambda_1-\lambda_2)+2\frac{1}{1!}\lambda_1\lambda_2\exp(-\lambda_1-\lambda_2)}{\frac{1}{2!}(\lambda_1+\lambda_2)^2\exp(-\lambda_1-\lambda_2)}=1+2\frac{\lambda_1\lambda_2}{(\lambda_1+\lambda_2)^2}
$$
Experimental results seem to show a logarithmic relationship between $M$ and $L$ ($L \approx a + b \log(M)$), but the explicit equation or even some approximation are both so hard to find by myself. Any hint on this?
 A: I think your computation for $N=2$ and $M=2$ is not correct, given that your final result does not depend on $\lambda_1, \lambda_2$.

Let $\Lambda_N := \sum_{j=1}^N \lambda_j$ and let $S_N := \sum_{j=1}^N X_j$.
\begin{align}
& E\left[\sum_{i=1}^N \mathbf{1}_{X_i > 0} \;\middle|\; \sum_{j=1}^N X_j = M\right]
\\
&= \sum_{i=1}^N P(X_i > 0 \mid S_N = M)
\\
&= N - \sum_{i=1}^N P(X_i = 0 \mid S_N = M)
\\
&= N - \sum_{i=1}^N \frac{P(X_i = 0 , S_N = M)}{P(S_N = M)}
\\
&= N - \sum_{i=1}^N \frac{P(X_i = 0) P(\sum_{j \ne i} X_j = M)}{P(S_N = M)}
\\
&= N - \frac{1}{e^{-\Lambda_N} \Lambda^M/M!} \sum_{i=1}^N e^{-\lambda_i} e^{-(\Lambda_N - \lambda_i)} \frac{(\Lambda_N - \lambda_i)^M}{M!}
\\
&= N - \sum_{i=1}^N \left(1 - \frac{\lambda_i}{\Lambda_N}\right)^M
\end{align}

Computational verification.
For $\lambda_1=\lambda_2=1$ and $M=2$, my expression above is $2 - \frac{1}{2^2} - \frac{1}{2^2} = 1.5$, and my simulations match this.
For $\lambda_1=1$, $\lambda_2=2$, $\lambda_3=3$ and $M=7$, my expression above is $3 - (5/6)^7 - (4/6)^7 - (3/6)^7 \approx 2.655$, and my simulations match this.
import numpy as np
import pandas as pd
from scipy import stats

def generate_sample(lambdas):
  return [stats.poisson.rvs(lam, size=1)[0] for lam in lambdas]

def generate_sample_conditioned(lambdas, sum_constraint):
  DEADLINE = 1000
  for _ in range(DEADLINE):
    sample = generate_sample(lambdas)
    if np.sum(sample) == sum_constraint:
      return sample

def run_experiment(lambdas, sum_constraint, num_trials):
  num_nonzero = []
  for i in range(num_trials):
    sample = generate_sample_conditioned(lambdas, sum_constraint)
    num_nonzero.append(np.sum(np.array(sample) > 0))
  return np.mean(num_nonzero)

def compute_theoretical(lambdas, sum_constraint):
  return len(lambdas) - np.sum((1 - lambdas / np.sum(lambdas)) ** sum_constraint)

lambdas = [1, 1]
M = 2
print(compute_theoretical(lambdas, M))
print(run_experiment(lambdas, M, 10000))

1.5
1.4992

lambdas = [1, 2, 3]
M = 7

print(compute_theoretical(lambdas, M))
print(run_experiment(lambdas, M, 10000))

2.6545781893004112
2.6506

A: The distribution of $n$ Poisson independent variables conditioned to the sum is a multinomial (proof). And the marginal of a multinomial is a binomial.
That is, if $Y = \sum_{i=1}^n X_i$, then $X_i | Y$ follows a Binomial $(Y, p_i)$ distribution with $p_i=\lambda_i/\sum \lambda_i$.
Letting $Z_i = 1$ if $X_i>0$ , $Z_i =0$ otherwise, we have
$$E[Z_i|Y]=1-P(Z_i=0|Y)=1-(1-p_i)^Y$$
and
$$E[L]=E[\sum Z_i | Y] = n- \sum_{i=1}^n (1-p_i)^Y$$
In particular, for iid variables, $\lambda_i=\lambda$:
$$E[L]=n \left( 1 - \left(1-\frac{1}{n}\right)^Y\right)$$
