What are the properties of the distribution of medians from samples of a set of integers? Given a set of integers $S = \{1, 2, 3, ..., n\}$ I draw uniformly at random without repetition $m$ integers a large number of times ($1 \leq m \leq n$). For each of these samples (each of size $m$) I compute the median. What is the expected distribution of these medians?
By intuition, we can say the sample space from which we can take the samples is the set of all permutations without repetition of the elements of $S$ of size $m$. Different permutations of this set will have the same median making certain median values more likely than others so the resulting distribution is not uniform for $1 \lt m \lt n$. For the extreme cases $m = 1$ we should observe a uniform distribution in the range $[1, n]$ and for $m = n$ all medians will be the median of $S$. I am unsure on what to state about the distribution of medians for all other values though.
 A: These order statistics follow a shifted beta-binomial distribution
Let $X_1,...,X_m \sim \text{SRSWOR}(n)$ denote the sampled values.  If you have an odd number of values in the sample then the median is just the middle order statistic.  I am going to generalise from this and show how you can obtain the distribution of any order statistic $X_{(k)}$ over the order values $k =1,...,m$.  (Getting the median in the case where you have an even number of values is a bit trickier, since it involves the joint distribution of two adjacent order statistics; I will leave this as an exercise.)
There are a couple of different ways you can derive the distribution for this order statistic.  The simplest method is through a direct combinatorial argument.  There are ${n \choose m}$ possible samples that can be drawn from the population (with equal probability).  To get the event $X_{(k)} = x$ we have ${x-1 \choose k-1}$ ways of choosing the previous order statistics and ${n-x \choose m-k}$ ways of choosing the subsequent order statistics.  Consequently, we get a shifted beta-binomial distribution:
$$\mathbb{P}(X_{(k)} = x)
= \frac{{x-1 \choose k-1} {n-x \choose m-k}}{{n \choose m}}
= \text{BetaBin}(x-k|n-m,k,m-k+1).$$
(To confirm that this is a valid probability mass function ---i.e., that it sums to one over the support--- note this related summation result for products of binomial coefficients.)  It is also possible to obtain this result by using the relationship of the order statistics to the hypergeometric distribution, but this latter method is much more complicated (omitted here).
With a bit of algebra, it is possible to derive the general form for the factorial moments of this distribution, and through this obtain the mean and variance.  These are given by:
$$\mathbb{E}(X_{(k)}) = \frac{n+1}{m+1} \cdot k
\quad \quad \quad \quad \quad 
\mathbb{V}(X_{(k)}) = \frac{(n+1)(n-m)}{(m+1)^2 (m+2)} \cdot k(m-k+1).$$
Further examination of this distribution can be found in O'Neill (2022).  The distribution of the scaled order statistics can be shown to have a mixture distribution that is closely related to the order statistics from IID observations from a continuous uniform distribuion on the unit interval.

Implementation in R: The beta-binomial distribution is available in various packages in R (e.g., the extraDistr and TailRank packages).  It is also quite simple to program the specific variant shown above (using the required parameterisation and shift) to give the distribution of the order statistic, which I will give below as the FPOS distribution (standing for finite-population-order-statistic).  Appealing to standard notation in these problems, I will use N for the population size and n for the sample size.  The density function is programmed using computation in log-space, in order to ensure computational stability.
#Probability mass function for FPOS distribution
dFPOS <- function(x, k, n, N, log = FALSE) {

  #Check input types
  if (!is.numeric(x))    stop('Input x must be numeric')
  if (!is.numeric(k))    stop('Order k must be a positive integer')
  if (!is.numeric(n))    stop('Sample size n must be a positive integer')
  if (!is.numeric(N))    stop('Population size N must be a positive integer')
  if (!is.logical(log))  stop('Input log must be a logical value')
  ki <- as.integer(k)
  ni <- as.integer(n)
  Ni <- as.integer(N)
  if (k != ki)           stop('Order k must be a positive integer')
  if (n != ni)           stop('Sample size n must be a positive integer')
  if (N != Ni)           stop('Population size N must be a positive integer')

  #Check input lengths and values
  if (length(k) != 1)    stop('Order k should be a single positive integer')
  if (length(n) != 1)    stop('Sample size n should be a single positive integer')
  if (length(N) != 1)    stop('Population size N should be a single positive integer')
  if (length(log) != 1)  stop('Input log should be a single logical value')
  if (k < 1)             stop('Order k must be positive')
  if (n < 1)             stop('Sample size n must be positive')
  if (N < 1)             stop('Population size N must be positive')
  if (k > n)             stop('Order k cannot be greater than sample size n')
  if (n > N)             stop('Sample size n cannot be greater than population size N')
  
  #Compute log-probabilities
  LOGPROBS <- rep(0, N)
  for (t in 1:N) {
    LOGPROBS[t] <- lchoose(t-1, k-1) + lchoose(N-t, n-k) - lchoose(N, n) }
  
  #Generate output
  L <- length(x)
  OUT <- rep(-Inf, L)
  for (i in 1:L) {
    xx <- x[i]
    if (xx %in% (1:N)) {
      OUT[i] <- LOGPROBS[xx] } }
  
  #Return output
  if (log) { OUT } else { exp(OUT) } }

Using this function we can obtain the distribution of the median (for an odd sample size) by taking $m=2k-1$ and finding the distribution of the order statistic $X_{(k)}$.  In the code below we compute and plot the mass function for the median for the case where $n=16$ and $m=9$.
#Set parameters
n <- 16
m <- 9

#Compute and plot the mass function
PROBS <- dFPOS(1:n, k=(m+1)/2, N = n, n = m)
names(PROBS) <- 1:n
barplot(PROBS, xlab = 'Value', ylab = 'Probability', col = 'blue')


A: First, consider the situation where the sample size $m=2k+1$ is odd so that the median is uniquely defined.
For $x\in S,$ the chance the median $M$ equals $x$  is the chance that at $k$ of the sample values are from the set $\{1,2,\ldots, x-1\}$ and another $k$ of the sample values are from $\{x+1,\ldots, n\},$ with the remaining value equal to $x$ itself.  Such samples are in one-to-one correspondence with ordered pairs $(I,J)$ where $I$ is a $k$-subset of $x-1$ values and $J$ is a $k$-subset of $n-x$ values.  Since there are $\binom{n}{2k+1}$ equiprobable samples, this chance therefore is

$$\Pr(M=x) = \frac{\binom{x-1}{k}\binom{n-x}{k}}{\binom{n}{m}} ;\ x \in S, m=2k+1.$$

When $m=2k$ is even, let the median be the midpoint of the interval that splits the sample into two equal halves.  When $2 \le 2x \le 2n$ is an integer, the median equals $x$ when the two closest numbers to $x$ in the sample are equally far from $x$--let's call them $i\lt x$ and (therefore) the other is $2x-i \gt x$--and the remaining $2(k-1)$ elements of the sample are split into two equal-size groups.
To restate this a little more precisely, there must be a (unique) integer $i \lt x$ for which (a) $i$ and $2x-i$ are in the sample, (b) exactly $k-1$ elements of the sample are in the set $\{1,2,\ldots, i-1\}$ (implying $i\ge k$) and (c) exactly $k-1$ elements are in the set $\{2x-i+1, 2x-i+2, \ldots, n\}$ (a set with $n - (2x-i)$ elements).  Because $i$ is unique, the events so described are disjoint: this implies their probabilities sum to the chance that $x$ is the median.  Thus

$$\Pr(M = x) = \frac{1}{\binom{n}{m}}\sum_{i=k}^{\lfloor x\rfloor} \binom{i-1}{k-1}\binom{n-2x+i}{k-1};\ m=2k, 2x\le n+1.$$

Because the distribution is symmetric about $(n+1)/2,$ use $n+1-x$ in place of $x$ in this formula when $2x \gt n+1.$
Here are some plots of the median distributions for various sample sizes from  $n=39$ values:

The alternating probabilities of the half-integral and integral medians for even sample sizes $m$ are apparent.
Simulations confirm the correctness of these formulas using a chi-squared test.
n <- 7
N <- 5e2
for (m in 1:n) {
  X.2 <- replicate(N, 2*median(sample.int(n, m)))
  p <- pmed(seq(1/2, n, by=1/2), n, m)
  i <- p > 0
  observed <- tabulate(X.2, nbins=2*n)[i]
  p <- p[i]
  Y <- rbind(observed, expected = N*p)
  colnames(Y) <- which(i)/2
  
  cat("m =", m, "; n =", n, "\n")
  print(signif(Y, 3))
  
  if (sum(i) > 1) print(chisq.test(observed, p=p))
}

Typical output tallies the simulated medians, compares them to the formulas, and displays the results of the chi-squared test.

m = 4 ; n = 7 
          2.5    3 3.5    4 4.5    5  5.5
observed 69.0 34.0 107 67.0 118 45.0 60.0
expected 57.1 42.9 114 71.4 114 42.9 57.1

  Chi-squared test for given probabilities

data:  observed
X-squared = 5.4006, df = 6, p-value = 0.4936


It performs the calculations in R using this quick and dirty implementation of the two formulas.  Using logarithms of the binomial coefficients makes the calculation feasible for very large $n$ (millions or more), whereas computing with the binomial coefficients themselves would overflow double-precision floats once $n$ exceeds $1000$ or so.
pmed <- Vectorize(function(x, n, m) {
  k <- floor(m/2)
  if (!isTRUE(m == 2*k)) {
    if (x == round(x)) {
      exp(lchoose(x-1, k) + lchoose(n-x, k) - lchoose(n, m))
    } else {
      0
    }
  } else {
    x <- pmin(x, n+1-x)
    x.0 <- floor(x-1/2)
    if(x.0 >= k) {
      i <- seq(k-1, x.0-1, by=1)
      q <- lchoose(i, k-1) + lchoose(n+1-2*x+i, k-1) 
      sum(exp(q - lchoose(n, m)))
    } else {
      0
    }
  }
}, "x")

