$\newcommand{\R}{\mathbb{R}} \newcommand{\ep}{\epsilon} $ Apparently, the still best available general result from which your desired bound will follow is the now pretty old result due to Talagrand, Theorem 1.3, part (ii). Indeed, the probability in the tiltle of your question can be rewritten as \begin{equation*} P\Big(\sup_{\|u\|_\infty\le1/k}\Big(\sum_1^nf_u(X_i)-n\,Ef_u(X_1)\Big)\le t\sqrt n\Big), \end{equation*} where $f_u(x):=(1+u\cdot x)/2$, $u\cdot x:=u_{(1)}x_{(1)}+\dots+u_{(k)}x_{(k)}$ and $\|u\|_\infty:=|u_{(1)}|+\dots|u_{(k)}|$$\|u\|_\infty:=\max(|u_{(1)}|,\dots,|u_{(k)}|)$ for $u=(u_{(1)},\dots,u_{(k)}),x=(x_{(1)},\dots,x_{(k)})\in\R^k$, \begin{equation*} t:=\frac{\ep_n\sqrt n}{2kM}, \tag{1} \end{equation*} $X_i:=(X_{i1},\dots,X_{i,k})$, $X_{ij}:=1_{Y_i=j}$, $Y_1,\dots,Y_n$ are iid random variables such that $P(Y_1=j)=p_{(j)}$ for $j=1,\dots,k$, and $(p_{(1)},\dots,p_{(k)}):=p$.
Applying the mentioned theorem by Talagrand, we get your desired inequality with \begin{equation*} \delta=(Kt/\sqrt k)^k e^{-2t^2}, \end{equation*} where $K>0$ is an absolute real constant and $t$ is as in (1).