Return to Answer

Brown, L. D., Cai, T. T., & DasGupta, A. (2001). Interval estimation for a binomial proportion. Statistical science, 101-117.

You can also approach this problem differently if you notice that under uniform prior $p$ follows $\mathrm{Beta}(n/2+1, n/2+1)$ distribution (Pires and Amado, 2008, see also here). Knowing this, you can find such $n$ that maximizes probability of $p$ being within $\pm \varepsilon$ of $0.5$, i.e. that $ \Pr(X \le m) \geq p-\varepsilon $ and $ \Pr(X \le m) \leq p+\varepsilon $.

betapr <- function(n, eps) diff(pbeta(c(0.5-eps, 0.5+eps), n/2+1, n/2+1))
betapr <- Vectorize(betapr, "n")

That for $\varepsilon = 0.05$ returns

> betapr(n, eps = 0.05)
 [1] 0.1271114 0.2020283 0.2662452 0.3162281 0.3942332 0.5281417 0.6875108 0.8881346 0.9752351 0.9984892

So with small sample probability that $p$ is anywhere close to $0.5$ is quite small. Beta distribution can be also used for obtaining exact confidence intervals.

Brown, L. D., Cai, T. T., & DasGupta, A. (2001). Interval estimation for a binomial proportion. Statistical science, 101-117.

Pires, A. M., & Amado, C. (2008). Interval estimators for a binomial proportion: Comparison of twenty methods. REVSTAT‒Statistical Journal, 6(2), 165-197.

Another thing that you can do is to find sample size that enables you to estimate median with enough precision. If you think of this problem the other way around, we want $np$ out of $n$ values of $X$ be smaller or equal than $m$. If we define new random variable in terms of counts $Y = \sum_{i=1}^n [x_i \le m ]$, then it follows binomial distribution with parameters $n$ and $p$. Wald confidence interval for $p$ can be easily calculated as

set.seed(123)
 
R <- 1e3
N <- 1e6
n <- c(1, 5, 10, 15, 25, 50, 100, 250, 500, 1000)
 
res <- matrix(NA, R, length(n))

for (i in 1:R) {
  for (j in seq_along(n)) {
    m <- median(rnorm(n[j]))
    res[i, j] <- pnorm(m)
  }
}

wald <- function(n, p=0.5) matrix(c(-1, 1), ncol = 1) %*% (1.96 * sqrt(p*(1-p)/n)) + p

> apply(res, 2, quantile, c(0.025, 0.975))
       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
2.5%  0.022 0.135 0.233 0.264 0.308 0.368 0.406 0.436 0.453  0.47
97.5% 0.980 0.861 0.771 0.728 0.682 0.626 0.594 0.561 0.540  0.53
> wald(n)
      [,1]   [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,] -0.48 0.0617 0.19 0.247 0.304 0.361 0.402 0.438 0.456 0.469
[2,]  1.48 0.9383 0.81 0.753 0.696 0.639 0.598 0.562 0.544 0.531

Regarding your question about sample size for estimating median, it can be easily calculated. If you think of this problem the other way around, we want $np$ out of $n$ values of $X$ be smaller or equal than $m$. If we define new random variable in terms of counts $Y = \sum_{i=1}^n [x_i \le m ]$, then it follows binomial distribution with parameters $n$ and $p$. Wald confidence interval for $p$ can be easily calculated as

set.seed(123)
R <- 1e5
n <- c(1, 5, 10, 15, 25, 50, 100, 250, 500, 1000)
res <- matrix(NA, R, length(n))

for (i in 1:R) {
  for (j in seq_along(n)) {
    m <- median(rnorm(n[j])) # sample median
    res[i, j] <- pnorm(m)    # population Pr(X <= m)
  }
}

wald <- function(n, p=0.5) matrix(c(-1, 1), ncol = 1) %*% (1.96 * sqrt(p*(1-p)/n)) + p

> wald(n)
      [,1]   [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,] -0.48 0.0617 0.19 0.247 0.304 0.361 0.402 0.438 0.456 0.469
[2,]  1.48 0.9383 0.81 0.753 0.696 0.639 0.598 0.562 0.544 0.531
> apply(res, 2, quantile, c(0.025, 0.975))
        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
2.5%  0.0258 0.146 0.233 0.266 0.313 0.366 0.403 0.438 0.456 0.469
97.5% 0.9745 0.853 0.769 0.733 0.687 0.635 0.596 0.562 0.544 0.531

As you can see, if $n=1$ then sample median could be possibly any value of $X$, as sample grows, sample median approaches population median, where it is for you to decide if it is precise enough.

set.seed(123)

R <- 1e3
N <- 1e6
n <- c(1, 5, 10, 50, 100, 250, 1000)
k <- length(n)
res <- matrix(NA, R, n[k])

# WARNING: very slow
for (i in 1:R) {
  print(i/R)
  X <- rnorm(N)
  x_small <- sample(X, n[k])
  for (j in 1:k) {
    m <- median(x_small[1:n[j]])
    res[i, j] <- sum(X <= m)
  }
}

res <- res/N

wald <- function(n, p=0.5) matrix(c(-1, 1), ncol = 1) %*% (1.96 * sqrt(p*(1-p)/n)) + p

> wald(n)
      [,1]   [,2] [,3]  [,4]  [,5]  [,6]  [,7]
[1,] -0.48 0.0617 0.19 0.361 0.402 0.438 0.469
[2,]  1.48 0.9383 0.81 0.639 0.598 0.562 0.531
> apply(res, 2, quantile, c(0.025, 0.975))
        [,1]  [,2] [,3]  [,4]  [,5]  [,6]  [,7]
2.5%  0.0234 0.152 0.25 0.366 0.403 0.437 0.467
97.5% 0.9817 0.831 0.76 0.631 0.597 0.564 0.529

set.seed(123)

R <- 1e3
N <- 1e6
n <- c(1, 5, 10, 15, 25, 50, 100, 250, 500, 1000)

res <- matrix(NA, R, length(n))

for (i in 1:R) {
  for (j in seq_along(n)) {
    m <- median(rnorm(n[j]))
    res[i, j] <- pnorm(m)
  }
}

wald <- function(n, p=0.5) matrix(c(-1, 1), ncol = 1) %*% (1.96 * sqrt(p*(1-p)/n)) + p

> apply(res, 2, quantile, c(0.025, 0.975))
       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
2.5%  0.022 0.135 0.233 0.264 0.308 0.368 0.406 0.436 0.453  0.47
97.5% 0.980 0.861 0.771 0.728 0.682 0.626 0.594 0.561 0.540  0.53
> wald(n)
      [,1]   [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,] -0.48 0.0617 0.19 0.247 0.304 0.361 0.402 0.438 0.456 0.469
[2,]  1.48 0.9383 0.81 0.753 0.696 0.639 0.598 0.562 0.544 0.531