If I have two random variables: N follows Poisson distribution with mean as $\lambda$, and X follows binomial distribution with success rate of p, the output of N is the number of trials for X, according to this link A binomial multiplied by a poisson, XN should follow Poisson distribution with mean $p\lambda$
However, I was reading "Probability for Statistics and Machine Learning" for DasGupta and I found this example
I do not know how those two cases differ, but I think I am missing a fundamental point.
Best Answer
A hierarchical model is not the same as the product of two random variables.
Here is a simple example of a product: Let $$\Pr[X = x] = \theta^x (1-\theta)^{1-x}, \quad x \in \{0,1\}, \\ \Pr[Y = y] = p(1-p)^{y-1}, \quad y \in \{1, 2, 3, \ldots\}.$$
So $X$ is Bernoulli with parameter $\theta$ and $Y$ is geometric with parameter $p$. Now $XY$ is a new random variable: clearly, $$\Pr[XY = 0] = \Pr[X = 0] = 1 - \theta,$$ since $\Pr[Y = 0] = 0$. And for each positive integer $y$, $$\Pr[XY = y] = \theta p(1-p)^{y-1},$$ since $XY = y$ only if $X = 1$ and $Y = y$.
Now here is an example of a hierarchical model. Let $$\Pr[N = n] = e^{-\lambda} \frac{\lambda^n}{n!}, \quad n \in \{0, 1, 2, \ldots\}, \\ \Pr[X = x \mid N = n] = \binom{n}{x} p^x (1-p)^{n-x}, \quad x \in \{0, 1, 2, \ldots, n\}.$$ Note that the second distribution is conditional on the first: in other words, $N$ is Poisson with rate $\lambda$, but $X$ given $N = n$ is binomial with $n$ trials with success probability $p$. The product of $N$ and $X$ is not the same random variable as the marginal (or unconditional) distribution $X$, which is $$\begin{align} \Pr[X = x] &= \sum_{n=x}^\infty \Pr[X = x \mid N = n] \Pr[N = n] \\ &= \sum_{n=x}^\infty \binom{n}{x} p^x (1-p)^{n-x} e^{-\lambda} \frac{\lambda^n}{n!} \\ &= \sum_{n=x}^\infty \frac{1}{x! (n-x)!} e^{-p\lambda} e^{-(1-p)\lambda} p^x (1-p)^{n-x} \lambda^x \lambda^{n-x} \\ &= e^{-p\lambda} \frac{(p\lambda)^x}{x!} \sum_{n=x}^\infty e^{-(1-p)\lambda} \frac{((1-p)\lambda)^{n-x}}{(n-x)!} \\ &= e^{-p \lambda} \frac{(p\lambda)^x}{x!}. \end{align}$$ by the law of total probability. Therefore, the unconditional distribution of $X$ is Poisson with rate $p\lambda$. This would imply, for instance, $\Pr[X = 1] = p\lambda e^{-p\lambda}$. However, as I have stated, $NX$ is not the same thing: $$\begin{align} \Pr[NX = 1] &= \Pr[(N = 1) \cap (X = 1)] \\ &= \lambda e^{-\lambda} p. \end{align}$$
So when you write "$XN$ should follow a Poisson distribution with mean $p\lambda$," you are incorrect. The question you linked to contains a similar misunderstanding, hence your confusion. That question talks about multiplying two probability mass functions, not random variables. And then the answer that is given assumes that this multiplication comes in the form of a hierarchical model, which is not necessarily true.
As for the question that you specifically posted, regarding the eggs, this is clearly a binomial-Poisson hierarchical model; i.e., $N$ gives the number of laid eggs, and $X \mid N$ gives the number of eggs that fertilized. Then the covariance between $X$ and $N$ is calculated by the formula $$\operatorname{Cov}[X,N] = \operatorname{E}[XN] - \operatorname{E}[X]\operatorname{E}[N],$$ and in this formula, $XN$ truly means the product of the random variables $X$ and $N$. And now, by the law of total expectation, we have $$\operatorname{E}[XN] = \operatorname{E}[\operatorname{E}[XN \mid N]] = \operatorname{E}[N \operatorname{E}[X \mid N]] = \operatorname{E}[N (Np)] = p \operatorname{E}[N^2] = p (\lambda^2 + \lambda).$$
The first two equalities are simply conditioning on the outcome of $N$; i.e., it is a shorthand way to write $$\operatorname{E}[XN] = \sum_{n=0}^\infty \operatorname{E}[nX \mid N = n]\Pr[N = n] = \sum_{n=0}^\infty \left( n \operatorname{E}[X \mid N] \right)\Pr[N = n].$$ The third equality arises from the fact that, conditioned on $N$, $X \mid N$ is binomial and has expectation $Np$. The fourth equality simply factors the constant parameter $p$ out of the expectation, leaving the final equality which is the second moment of a Poisson random variable.