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.

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")

