Let $p_1, ... ,p_k$ denote the probabilities of drawing bin $1, .. ,k$, where $\sum_{i = 1}^{k} p_i= 1$. My question is if we draw $n$ times, how can I show that the probability that all bins are drawn less than $C$ times is maximized by $p_1 = p_2 = ... = p_k = 1/k$ ?
I've tried using conditional probability to solve inductively such as assuming bin $i$ is drawn $<C$ times and trying to reason about the remaining $k-1$, but no luck so far. Computing explicitly seems to be really difficult.
Here is a proof using generating functions. Let $c\geq 2$ and $k\geq 2$ be fixed.
Let $X=(X_1(n),\ldots,X_k(n))$ be the $k$-tuple of occupancy numbers at time $n$, i.e. $X_i(n)$ = number of bins of type $i$ drawn at time $n$.
Clearly $X$ has the multinomial distribution with parameters $n$ and $p=(p_1,\ldots,p_k)$.
Therefore
$$\mathbb{P}( \max X_i(n)\leq c-1)=n!\,[t^n] \prod _{i=1}^n q_c(p_it)$$
where $q_c(t):=\sum_{j=0}^{c-1}\frac{t^j}{j!}$ is the $c$th partial sum of $\exp(t)$.
The following lemma is the key.
Lemma
When $p_i,p_j$ are each replaced by their arithmetic mean $a:=(p_i+p_j)/2$
$\mathbb{P}(\max X_i(n) \leq c-1 )$ will not decrease.
If $c\leq n \leq k(c-1)$ and all other $p_\ell$ are positive $\mathbb{P}(\max X_i(n) \leq c-1 )$ will strictly increase if
$p_i\not = p_j$.
Proof Consider first the case of two factors. Let $x,y\in\mathbb{R}_+,x<y$. We show that $$ [t^m] q_c(xt)q_c(yt)\leq [t^m] \left(q_c(((x+y)/2)t\right)^2\;\;.$$It is easy to see that equality holds for $m\leq c-1$ and $m>2c-2$. Let $c\leq m \leq 2c-2$. We have $$m!\,[t^m] q_{c}(xt)q_c(yt)=(x+y)^m-\sum_{i=0}^{m-c} {m \choose i}\left(x^i y^{m-i}+y^i x^{m-i}\right) $$ For fixed sum $s=x+y$ the function $x \mapsto f(x):=\sum_{i=0}^{m-c} {m \choose i}\left(x^i y^{m-i}+y^i x^{m-i}\right)$ has the derivative $f^\prime(x)=(m-c+1){m \choose c} \left(x^{m-c+1}(s-x)^c-(s-x)^{m-c+1}x^c\right)$. Thus $x\mapsto f(x)$ is strictly decreasing (resp. increasing) on $[0,s/2]$ (resp. $[s/2,s]$), attaining its minimum at $x=s/2$. The rest is easy. End proof.
For any distribution $p=(p_1,\ldots,p_k)$ we have that $\mathbb{P}( \max X_i(n)\leq c-1)=1$ for $n\leq c-1$, and $\mathbb{P}( \max X_i(n)\leq c-1)=0$ for $n > k(c-1)$, so that only the cases $c\leq n \leq k(c-1)$ are of interest.
Corollary
For $c\leq n \leq k(c-1)$ the uniform distribution $p_1=\ldots=p_k=\tfrac{1}{k}$
uniquely maximises $\mathbb{P}(\max X_i(n) \leq c-1 )$.
Proof: Since we are maximising a continuous function (multivariate polynomial) on a compact set (simplex) the maximum is attained. Let $p$ be a maximising distribution. No $p_i$ can be $0$ (else replacing the last $0$ strictly increases the probability). But then $p$ can only be maximising if $p$ is the uniform distribution. End proof.
Remark Lemma 2 shows that the probability $\mathbb{P}(\max X_i(n) \leq c-1 )$ is a Schur-concave function on the simplex $p_i\geq 0, \sum_{i=1}^k p_i=1$. Possibly the result can be found under this heading in the literature (but a brief Google-search didn't reveal anything).
Experts may have a better answer.
This is I believe very hard if not impossible to do exactly, but if you use poissonization, i.e., model $X_i$ as independent poisson arrivals into bin $i$ with rate $np_i$ then the new process $(X_1,\ldots, X_k)$ equals the original bin arrival process $(Y_1,\ldots,Y_k)$ in expectation.
There is also a property that the probability of a monotone event can only be at most a factor of $2$ different in the bin arrival model and this model (I saw this in Mitzenmacher and Upfal's Probability and Computing textbook). Your event
$$ A=\{ \omega: \textrm{Max } X_i \leq C\} $$ is monotone. Maybe you can show strong concentration with upper and lower bounds on $A$ in the poisson model.
