In general, if $p_1, p_2, \ldots, p_m \in (0,1)$ are IID realizations from some probability distribution with mean $\mu_p$ and standard deviation $\sigma_p$, and $n_1, n_2, \ldots, n_m \in \mathbb Z^+$ are IID realizations of from another probability distribution with mean $\mu_n$ and standard deviation $\sigma_n$, and for each $i = 1, 2, \ldots, m$, we have random variables $$X_i \sim \operatorname{Binomial}(n_i, p_i),$$ and we are interested in the distribution of $S = \sum_{i=1}^m X_i$, then we have by linearity of expectation $$\operatorname{E}[S] = \sum_{i=1}^m \operatorname{E}[X_i].$$ In turn, for each $X_i$, we have by the law of total expectation $$\operatorname{E}[X_i] = \operatorname{E}[\operatorname{E}[X_i \mid (n_i \cap p_i)]] = \operatorname{E}[n_i p_i] = \operatorname{E}[n_i]\operatorname{E}[p_i] = \mu_n \mu_p;$$ thus $$\operatorname{E}[S] = m\mu_n \mu_p.$$ This assumes that $n_i$ and $p_i$ are independent for each $i$ (from which it follows that each $X_i$ is independent). The variance calculation is done in a similar fashion; $$\operatorname{Var}[S] \overset{\text{ind}}{=} \sum_{i=1}^m \operatorname{Var}[X_i],$$ whence by the law of total variance
$$\begin{align*}
\operatorname{Var}[X_i]
&= \operatorname{Var}[\operatorname{E}[X_i \mid (n_i \cap p_i)]] + \operatorname{E}[\operatorname{Var}[X_i \mid (n_i \cap p_i)]] \\
&= \operatorname{Var}[n_i p_i] + \operatorname{E}[n_i p_i (1-p_i)] \\
&= (\sigma_n^2 \sigma_p^2 + \sigma_n^2 \mu_p^2 + \sigma_p^2 \mu_n^2) + \mu_n \operatorname{E}[p_i(1-p_i)] \\
&= (\sigma_n^2 \sigma_p^2 + \sigma_n^2 \mu_p^2 + \sigma_p^2 \mu_n^2) + \mu_n (\mu_p - (\sigma_p^2 + \mu_p^2)).
\end{align*}$$
To understand the variance of $n_i p_i$, note that for two independent random variables $A$, $B$, with means and standard deviations $\mu_A, \sigma_A, \mu_B, \sigma_B$, respectively,
$$\begin{align*}\operatorname{Var}[AB]
&= \operatorname{E}[(AB)^2] - \operatorname{E}[AB]^2 \\
&= \operatorname{E}[A^2 B^2] - \operatorname{E}[A]^2 \operatorname{E}[B]^2 \\
&= \operatorname{E}[A^2]\operatorname{E}[B^2] - \mu_A^2 \mu_B^2 \\
&= (\operatorname{Var}[A] + \operatorname{E}[A]^2)(\operatorname{Var}[B] + \operatorname{E}[B]^2) - \mu_A^2 \mu_B^2 \\
&= (\sigma_A^2 + \mu_A^2)(\sigma_B^2 + \mu_B^2) - \mu_A^2 \mu_B^2 \\
&= \sigma_A^2 \sigma_B^2 + \sigma_A^2 \mu_B^2 + \sigma_B^2 \mu_A^2. \end{align*}$$
Note that my computation of the variance differs from yours. I have substantiated my results by simulating $m = 10^6$ observations from $X_i$ where $n_i \sim \operatorname{Poisson}(\lambda)$ and $p_i \sim \operatorname{Beta}(a,b)$, for $\lambda = 11$ and $(a,b) = (7,3)$. This should result in $\operatorname{Var}[X_i] = 1001/100$; your results do not match. I should also point out that the reason that your computation does not work is because the total variance of each $X_i$ is not merely due to the expectation of the conditional variance of $X_i$ given $n_i$ and $p_i$; the other term in the law of total variance must also be included, which captures the variability of the conditional expectation of $X_i$. In other words, there is variation in $X_i$ coming from the binomial variance even when $n_i$ and $p_i$ are fixed, but there is also additional variation in the location of $X_i$ arising from the fact that $n_i$ and $p_i$ are not fixed.
I'll just treat your example. I'm not sure how difficult it would be to express this in closed form; I haven't tried. We must distinguish between the cases where the last throw, the one that fulfills all the conditions, is a $1,2,$ or $3$. Suppose it is a $1$. Then we know that we threw $k\geq5$ $2$'s and $j\geq6$ $3$'s and that in the $k+j+3$ rolls prior to the last we rolled exactly $3$ $1$'s. We can make similar analyses when the last roll was a $2$ or a $3$. The expected number of rolls is
$$\sum_{k=5}^\infty\sum_{j=6}^\infty(k+j+4){k+j+3\choose3,k,j}3^{-(k+j+4)}+\\
\sum_{i=4}^\infty\sum_{j=6}^\infty(i+j+5){i+j+4\choose4,i,j}3^{-(i+j+5)}+\\
\sum_{i=4}^\infty\sum_{k=5}^\infty(i+k+6){i+k+5\choose5,i,k}3^{-(i+k+6)}
$$ where, of course, the first sum deals with the case where a $1$ is rolled last, the second where $2$ is last, and the third where $3$ is last.
Best Answer
Preliminary (TL;DR)
Background
In his 1991 publication, Norman C. Beaulieu answered your question w/ what he dubbed, the generalized multinomial distribution (GMD). My explanation will focus on the GMD's utility.
Notation
Pmf of GMD
$$P\left[X = x\right] = \sum_{\mathfrak{s} \in \mathfrak{S}_{([c], m)}} \left\{\prod_{k = 1}^t \left\{p_{k,\mathfrak{s}_k}\right\}\right\}$$
So far, I've identified it as being the superclass of 7 distributions! Namely...
Examples
Games
Questions
Answers w/o knowledge of GMD
$\Longrightarrow P\left[X = x\right] = 3/5$.
$\Longrightarrow P\left[X = x\right] = 1/6$.
$\Longrightarrow P\left[X = x\right] = 2/15$.
$\Longrightarrow P\left[X = x\right] = 20412/78125$.
$\Longrightarrow P\left[X = x\right] = 224/140625$.
$= \frac{1}{8}\sum_{i = 0}^7 \left\{\exp\left(\frac{-j5\pi i}{4}\right) \prod_{k = 1}^7 \left\{\left(\frac{0.5k + 14.5}{30}\right)\left(\exp\left(\frac{j\pi i}{4}\right) - 1\right) + 1\right\}\right\}$
$\Longrightarrow P\left[X_2 = 5\right] = 308327/1440000$.
Answers w/ Knowledge of GMD
Final Words
I know my answer was very long (& went far beyond what OP asked for) but this had been flying around inside my head for quite some time & this q seemed like the most suitable landing strip.
I performed the last 6 calculations using the function gmdPmf (which I defined in Mathematica)...
Please edit, if you know of any ways to make it shorter/faster. Congrats, if you made it to the end! (: