There is no particularly simple way to write this probability, but you can express it as a mixture using the MaxCount distribution (see this related answer) if you want. This is not necessarily a simplification, since the MaxCount distribution is itself similarly complicated and is usually computed recursively. Nevertheless, this way of expressing things makes it relatively simple to compute the probability of interest using existing probability functions that are available for this distributional family.
To facilitate analysis of this problem, let $M \equiv \max(X_2,...,X_k)$ denote the maximum frequency amongst the lower-probability sides, let $\mathscr{M}_{+} \equiv \{ X_1 > M \}$ denote the event where the biased side has the (strictly) highest frequency, and let $\mathscr{M}_0 \equiv \{ X_1 \geqslant M \}$ denote the event where the biased side has the (non-strictly) highest frequency. Using the law of total probability, the probabilities of these two events can be written in mixture form as:
$$\begin{align}
\mathbb{P}(\mathscr{M}_0)
&= \mathbb{P}(X_1 \geqslant M) \\[12pt]
&= \sum_{x=0}^n \mathbb{P}(X_1 \geqslant M | X_1 = x) \cdot \mathbb{P}(X_1 = x) \\[6pt]
&= \sum_{x=0}^n \mathbb{P}(M \leqslant x | X_1 = x) \cdot \mathbb{P}(X_1 = x) \\[6pt]
&= \sum_{x=\lceil n/k \rceil}^n \mathbb{P}(M \leqslant x | X_1 = x) \cdot \mathbb{P}(X_1 = x) \\[6pt]
&= \sum_{x=\lceil n/k \rceil}^n F_\text{MaxCount}(x|n-x,k-1) \cdot \text{Bin}(x|n,p), \\[18pt]
\mathbb{P}(\mathscr{M}_+)
&= \mathbb{P}(X_1 > M) \\[12pt]
&= \sum_{x=0}^n \mathbb{P}(X_1 > M | X_1 = x) \cdot \mathbb{P}(X_1 = x) \\[6pt]
&= \sum_{x=0}^n \mathbb{P}(M < x | X_1 = x) \cdot \mathbb{P}(X_1 = x) \\[6pt]
&= \sum_{x=\lceil (n+k-1)/k \rceil}^n \mathbb{P}(M < x | X_1 = x) \cdot \mathbb{P}(X_1 = x) \\[6pt]
&= \sum_{x=\lceil (n+k-1)/k \rceil}^n F_\text{MaxCount}(x-1|n-x,k-1) \cdot \text{Bin}(x|n,p). \\[6pt]
\end{align}$$
(Note that $\mathbb{P}(M \leqslant x | X_1 = x) = 0$ for $x < n/k$ and $\mathbb{P}(M < x | X_1 = x) = 0$ for $x < (n+k-1)/k$, which is the reason for the adjustment of the lower bound of the sum in the third step of each expansion. It is also okay to leave the lower bound at zero and include the zero terms in the sum if preferred.) This "simplification" frames the result in terms of the MaxCount distribution, which is also a complicated distribution that is computed recursively.
Computational Implementation: The above mixture form can be used for computation using the probability functions for the MaxCount distribution and binomial distribution. The maxcount probability functions are available in the occupancy package in R and the binom probability functions are in the standard packages. In the code below we create a function prob.maxfreq that computes the probability/log-probability of interest for any n, k and prob.
#Set function to compute probability of maximum frequency
prob.maxfreq <- function(n, k, prob, strict = FALSE, log = FALSE) {
#Check input n
if (!is.vector(n)) stop('Error: Input n should be a vector')
if (!is.numeric(n)) stop('Error: Input n should be a numeric vector')
if (length(n) != 1) stop('Error: Input n should be a single numeric value')
if (n != as.integer(n)) stop('Error: Input n should be an integer')
if (n < 1) stop('Error: Input n should be a positive integer')
#Check input k
if (!is.vector(k)) stop('Error: Input k should be a vector')
if (!is.numeric(k)) stop('Error: Input k should be a numeric vector')
if (length(k) != 1) stop('Error: Input k should be a single numeric value')
if (k != as.integer(k)) stop('Error: Input k should be an integer')
if (k < 1) stop('Error: Input k should be a positive integer')
#Check input prob
if (!is.vector(prob)) stop('Error: Input prob should be a vector')
if (!is.numeric(prob)) stop('Error: Input prob should be a numeric vector')
if (length(prob) != 1) stop('Error: Input prob should be a single numeric value')
if (prob > 1) stop('Error: Input prob should be a probability value')
if (prob <= 1/k) stop('Error: Input prob should be greater than 1/k')
#Check input strict
if (!is.vector(strict)) stop('Error: Input strict should be a vector')
if (!is.logical(strict)) stop('Error: Input strict should be a logical vector')
if (length(strict) != 1) stop('Error: Input strict should be a single logical value')
#Check input log
if (!is.vector(log)) stop('Error: Input log should be a vector')
if (!is.logical(log)) stop('Error: Input log should be a logical vector')
if (length(log) != 1) stop('Error: Input log should be a single logical value')
#Compute mixture parts
LOG.BINOM <- dbinom(0:n, size = n, prob = prob, log = TRUE)
LOG.ALL <- occupancy::dmaxcount.all(max.x = n, max.size = n, space = k-1, log = TRUE)
LOG.MAXCOUNT <- rep(-Inf, n+1)
if (strict) {
for (x in 1:n) {
LOG.MAXCOUNT[x+1] <- matrixStats::logSumExp(LOG.ALL[1:x, n-x+1]) }
} else {
for (x in 0:n) {
LOG.MAXCOUNT[x+1] <- matrixStats::logSumExp(LOG.ALL[1:(x+1), n-x+1]) } }
#Compute log-probability of interest
LOGPROB <- matrixStats::logSumExp(LOG.MAXCOUNT + LOG.BINOM)
#Return output
if (log) { LOGPROB } else { exp(LOGPROB) } }
We can now compute the probability of interest for some stipuled values. For example, taking $n=100$, $k=60$ and $p = 0.03$ we get the probability $\mathbb{P}(\mathscr{M}_0) = 0.1827769$. We can also estimate this probability from simulation of the sampling model to confirm that this computation is accurate.
#Set parameters
n <- 100
k <- 60
p <- 0.03
#Compute exact probability
prob.maxfreq(n = n, k = k, prob = p)
[1] 0.1827769
#Generate simulations of sampling mechanism
set.seed(1024817537)
N <- 10^6
PROB.VEC <- c(p, rep((1-p)/(k-1), k-1))
SIMS <- matrix(0, nrow = N, ncol = n)
SIMS.FREQ <- matrix(0, nrow = N, ncol = k)
MAXFREQ <- rep(FALSE, N)
for (i in 1:N) {
SIMS[i, ] <- sample.int(k, size = n, replace = TRUE, prob = PROB.VEC)
for (x in 1:k) {
SIMS.FREQ[i, x] <- sum(SIMS[i, ] == x) }
MAXFREQ[i] <- (SIMS.FREQ[i, 1] == max(SIMS.FREQ[i, ])) }
#Compute estimated probability from simulations
mean(MAXFREQ)
[1] 0.183028
