This is an interesting problem. It seems like finding an analytical solution for the confidence intervals will be very difficult, so I would instead try to look for reasonable estimates.
I would start by treating the $\hat p_i$ as i.i.d. samples from a population with unknown mean and unknown variance. This would give you confidence intervals of
$$\hat\mu\pm z\frac{\hat\sigma}{\sqrt n},$$
where $\hat\mu$ and $\hat\sigma$ are the sample mean standard deviation, $n$ is the sample size (your $M$) and $z$ is the $z$-score (i.e. for two-sided 95% confidence, $z_{0.05}=1.96$).
The accuracy of this confidence interval will depend on the validity of the implicit normality assumption. In this case, there are a few points on your side:
- The sample size is large
- The values of $\hat p_i$ are bounded by $[0,1]$
- $p_i=0$ more than half the time.
Combining these three facts justifies the use of the normality assumption for the lower confidence interval at least. This is because, the fact that $p_i=0$ more than half of the time implies that the sample will contain a large number of values with $\hat p_i=0$, and thus the sample will span the entire lower tail of the distribution. The normality assumption generally breaks down when there are very long tails in the population distribution which won't show up in a sample (warning, that is not a rigorous statement). If you could also guarantee that, say $p_i>0.5$ at least $5\%$ of the time, then you would be safe on both ends.
However, the upper confidence interval could be problematic in situations where the sample mean $\hat\mu$ is very close to $0$. Consider, for example, the extreme case where all of the $p_i$ are equal to $0$ except for one: $p_1=1$. Depending on the probability $q_1$ associated with this ball, we could easily have a case where, with probability, say $0.5$ this ball is never drawn during the sample, in which case we would have $\hat \mu=0$ but the confidence interval would be invalid.
To address this specific case, we can explicitly construct a confidence interval when the sample mean is $0$. If the total weight assigned to non-zero balls is $q$, we can see that the probability of never sampling such a ball is, $(1-q)^M$, and in order for the probability to be less than, say $0.05$, we must have $q\geq 1-0.05^{1/M}$.
If follows that, if $\hat\mu=0$, then we can say with $95\%$ confidence that the probability that a single round of the experiment returns a non-zero value is at most $1-0.05^{1/M}$, which means that, with $95\%$ confidence, the true mean lies in the interval
$$[0,1-0.05^{1/M}].$$
For your specific constants, this would yield an interval of $[0,0.0003]$.
Of course, there are intermediate scenarios. As a general rule of thumb,if the sample mean $\hat\mu$ is bounded away from $0$, then the standard confidence intervals should be pretty accurate, and the closer to $0$ that $\hat\mu$ is, the more you have to worry. I'd be tempted to suggest something like
$$[\hat\mu- z_\alpha\frac{\hat\sigma}{\sqrt n},
\hat\mu+ z_\alpha\frac{\hat\sigma}{\sqrt n} + 1-\alpha^{1/M}]
$$
as an "out-of-the-box" $1-\alpha$ confidence interval, but without closer examination, I'm not fully confident that this will work well in situations where $\hat\mu$ is small but non-zero.