Unbiased Estimator of Largest Mean of Two Normal Distributions

Question

Given samples from two normal distributions:

$X_i \stackrel{iid}{\sim} \mathcal{N}(\mu_X, \sigma_X^2)$ for $i = 1,...,n$

$Y_i \stackrel{iid}{\sim} \mathcal{N}(\mu_Y, \sigma_Y^2)$ for $i = 1,...,n$

How can I form an unbiased estimator of $\max(\mu_X, \mu_Y)$?

Clearly, $\overline{X}$ and $\overline{Y}$ are unbiased estimators of $\mu_X$ and $\mu_Y$ respectively, but $\mathbb{E}[\max(\overline{X},\overline{Y})] > \max(\mu_X, \mu_Y)$.

If an unbiased estimator does not exist, is there at least some way to reduce the bias of $\max(\overline{X},\overline{Y})$?

Answer 1

Finding an exactly unbiased estimator is probably impossible, so a practical solution is bootstrapping. I will here show nonparametric bootstrap, but small modification give a parametric bootstrap. So we assume the data is a sample from the distribution $F$, and the interest parameter is a function of $F$, $t(F)$. In your case $\theta=t(F) =\max(\mu_X,\mu_Y)$. This is estimated by the plug-in estimate $\hat{\theta}=\max(\bar{X}, \bar{Y})$. The bias of this estimator is $$ \DeclareMathOperator{\b}{bias} \DeclareMathOperator{\bh}{\hat{bias}} \DeclareMathOperator{\E}{\mathbb{E}} \b_F(\hat{\theta},\theta)=\E_F t(X,Y) - t(F) $$ which we can estimate under bootstrapping by $$ \bh_{\hat{F}} =\E_{\hat{F}} t(X^*,Y^*) -t(\hat{F}) =\frac1B \sum_i^B t(X_i^*,Y_i^*) -\hat{\theta} $$ where $B$ is the number of bootstrap resamples, and superscript $^*$ signifies a bootstrap resample. We can do this in R:

sigma <-3
mu1 <- 0
mu2 <- 0.67
n1 <- n2 <- 20
set.seed(7*11*13) # My public seed

# Simulate some observed data:

x1 <- rnorm(n1, mu1, sigma)
x2 <- rnorm(n2, mu2, sigma)

mu1_hat <- mean(x1)
mu2_hat <- mean(x2)
max_hat <- max(mu1_hat, mu2_hat)

### then for the bootstrapping

B <- 2000

myboot <- function(x1, x2, B) {
    max_hat <- max(mean(x1), mean(x2))
    # Then B bootstrap samples
    boots <- numeric(length=B)
    n1 <- length(x1) ; n2 <- length(x2)
    for (i in seq_along(boots)) {
        boots[i] <- max(mean(sample(x1, n1, replace=TRUE)),
                        mean(sample(x2, n2, replace=TRUE)))
    }
    bias_hat <- mean(boots) - max_hat
    return(list(mu1_hat=mean(x1), mu2_hat=mean(x2),
                max_hat=max_hat, bias_hat = bias_hat, boots=boots))
}

res <- myboot(x1, x2, B)

res[1:4]
$mu1_hat
[1] -0.007525858

$mu2_hat
[1] 0.8717599

$max_hat
[1] 0.8717599

$bias_hat
[1] 0.1455065

We can observe the skewness in the bootstrap distribution:

Here are some papers that treats a generalized version of the problem: this and this one, both papers uses multi-stage sampling.

---

In a comment I mentioned the jackknife as a possible solution. That does actually not work in this example, at least not the standard jackknife using leave-one-out, it would need some adapted jackknife leaving out more points, and is probably not worth the trouble. But it is interesting to think about why it does not work! As a help to that, I offer this example, continuing the example above:

xdata <- cbind(c(x1, x2), rep(1:2, c(n1, n2)))

theta <- function(x, xdata) {
    max(tapply(xdata[x, 1], xdata[x, 2], mean))
}

jackknife(1:(n1+n2), theta, xdata)
$jack.se
[1] 0.8787011

$jack.bias
[1] 0

$jack.values
 [1] 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599
 [8] 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599
[15] 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599 0.8717599 0.8921795
[22] 0.8671260 0.6355972 0.7516054 0.7971055 0.9763174 1.1825116 0.8518249
[29] 0.9742421 1.1263595 0.8067769 0.8461670 0.9041567 0.5017359 1.1859780
[36] 0.6198090 0.5474274 0.9121232 1.2174494 0.8387057

Answer 2

There is no unbiased estimator, for reasons essentially as in this answer or this answer.

Suppose we had the identity $$\max(\mu_x,\mu_y)=E[f(X,Y)]$$ and, in particular, had the right-hand side always defined. Then for a standard normal $Y$ we would have,

$$\max(\mu_x,0)=\!\iint \frac{f(X,Y)}{(2\pi)^n}\exp\left(-\frac{\sum(X-\mu_x)^2}{2\sigma_x^2}-\frac{\sum Y^2}{2}\right) dX\,dY$$

So the derivative of the left-hand side with respect to $\mu_x$ would be discontinuous: 0 when $\mu_x$ is negative, 1 when $\mu_x$ is positive, and undefined when $\mu_x$ is zero. But the derivative of the right-hand side would be a continuous function of $\mu_x$.

Since this is impossible, the original identity is also impossible, and there is no unbiased estimator for $\max(\mu_X,\mu_Y)$.