Mean of bootstrapped medians I have a large set (10^6) of very small sets (5-10 entries) of numbers. For each of them, I want to compute a sensible average measure. An example might be this:
$$
x = [-2, 5, -1, 3, 100]
$$
Because of the nature of the data, I know that the 100 is not an outlier, but it should also not be taken too seriously (hard to explain better without going into boring details). Hence, I want use use neither the median (which loses too much information), nor the arithmetic mean (which gives too much weight to the 100). Using the intercept of a robust regression $x$ ~ $1$ also effectively treats the 100 as an outlier.
So, I came up with the following approach:

*

*For a set with $n$ numbers ($n=5$ in the above example), select with replacement $n$ samples.

*Compute the median of these $n$ samples.

*Repeat steps 1 and 2 $m>>1$ times, and compute the mean of the $m$ medians.

For a set with $n$ numbers, there are $n^n$ possibilities for sampling with replacement, so if $n$ is small, this can also be done deterministically by exhausting all possibilities.
The results I obtained ($7.8096$ for the above example) correspond extremely well to what I would intuitively consider a good average measure, given what I know about the data.
I'm surprised I've never heard about this measure, given how natural it feels for computing a robust average that incorporates all data entries. Does anyone know more about it and can point me to relevant literature? Or more generally to information about the distribution of bootstrapped medians of finite sets of numbers?

In principle we will have to handle ties in the data, so let's adopt a convenient notation.  Let the data have values $x_1\lt x_2 \lt \cdots \lt x_d$ where the value $x_k$ appears with multiplicity $\nu_k \ge 1.$ Write $n_k = \nu_1 + \nu_2 + \cdots + \nu_k$ for the partial sums, so that $n = \nu_d$ is the dataset size.
You describe a random variable created by sampling $m = n$ values $X_1,X_2,\ldots, X_m$ from these data (with replacement) and computing their $q = 1/2$ quantile $X.$  Given any index $k,$ what is the chance that $x_k$ is this quantile?
Let's say that a "success" at step $i$ is that $X_i \le x_k.$  This event has probability $n_k/n.$  Because $X \le x_k$ exactly when there have been $qm$  or more successes when drawing the sample of size $m,$
$$\Pr(X \le x_k) = S\left(qm; m, n_k/n\right)$$
where $S$ is the cumulative Binomial survival function (computable as an incomplete Gamma integral or, for small $m,$ directly using Binomial coefficients).  Abbreviating the right hand side as $p_k,$ we obtain $\Pr(X = x_k) = p_k - p_{k-1}$ and the expectation is
$$E[X] = \sum_{i=1}^k x_k\left(p_k - p_{k-1}\right).$$
This is (explicitly) a linear combination of the order statistics of the original $n$ data values $y_1, y_2, \ldots, y_n$ where the coefficient of $y_i$ is $(p_k-p_{k-1})/\nu_k$ when $y_i = x_k.$  The highest weights tend to be concentrated near the $q$ quantile of the data.
To illustrate, on the left is a bar chart of $n=24$ data showing $(x_k, \nu_k).$  In the middle is a plot of the distribution $(x_k, p_k - p_{k-1})$ for $q = 30\%.$  This shows how the weights are concentrated around a few data values.  At the right is a histogram of the $30^\text{th}$ percentiles of a thousand independent bootstrap samples (using the default quantile calculation in R).  The difference between the average of these thousand quantiles and the computed expectation is insignificant ($Z = -0.005$).

In an application with no ties likely and predictable dataset sizes $n$ (such as the five to ten in the question) you can precompute the $p_k$ for each possible $n,$ making the calculation extremely efficient.  If there aren't too many ties in any case you wouldn't be wrong to just sort the data and use the precomputed $p_k.$
This R code generated the figures.
#
# Make some data.
#
set.seed(17)
nu <- table(y <- round(rnorm(24, -2, 3)))
#
# Preprocessing.
#
x <- as.numeric(names(nu))
n <- cumsum(nu)
n. <- n[length(n)]
#
# Specify the subsampling.
#
m <- n.
q <- 3/10
#
# Compute the expectation of X.
#
dp <- diff(c(0, pbinom(q * m, m, n/n., lower.tail = FALSE)))
dpm <- rev(diff(c(0, pbinom((1-q) * m, m, n/n., lower.tail = FALSE))))
dp <- (dp + dpm) / 2 # Assures agreement between `pbinom` and `quantile` conventions
xq.hat <- sum(x * dp)
#
# Confirm with a simulation.
#
n.sim <- 1e3
z <- replicate(n.sim, quantile(sample(y, m, replace = TRUE), q))
#
# Test the difference between the mean simulated value and the computation.
#
delta <- mean(z) - xq.hat
se <- sd(y) / sqrt(length(y))
Z <- delta / se
(signif(c(Delta = delta, SE = se, Z = Z), 2))
#
# Show what happened in plots.
#
par(mfrow = c(1,3))
plot(x, nu, type = "h", lwd = 5, main = "Original Data", 
     ylab = expression(nu), ylim = c(0, max(nu)))
plot(x, dp, type = "b", lwd = 2, xlab = "x", ylab = "Probability", 
     main = "Distribution of X", pch = 21)
hist(z, main = paste0("Simulated ", q, " Quantiles"), xlab = "X",
     xlim = range(x))
abline(v = xq.hat, lwd = 2, col = "Red")
par(mfrow = c(1,1))

