# Unbiased Estimator of Largest Mean of Two Normal Distributions

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})$$?

• The bias would depend on $\mid \mu_X-\mu_Y\mid$, being very small if that is large. Finding an exactly unbiased estimator will not be easy, practically go for bootstrap or jackknife ... Commented May 4, 2020 at 11:50
• How would you practically go about making an estimate using bootstrap or jackknife methods? Commented May 4, 2020 at 14:57
• Since $\max(a,b)=(a+b)/2+|a-b|/2$, it would be enough to get an unbiased estimator for $|\mu_x-\mu_y|/2$, and add that to the average of the estimators of the means.
– user225256
Commented Apr 30, 2021 at 19:33

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

• Nice answer. However, max_hat appears twice in the code (the first globally and the second locally). The first instance is not needed since it is repeated within 'myboot'. Commented May 8, 2020 at 16:52
• If you use rstan or other Bayesian modeling software based on Monte Carlo simulation to draw posterior distribution samples, you can forget about attempting to get "unbiased" estimates and quite simply compute the maximum of each pair of posterior draws for the two unknown parameters. This gives rise to the posterior distribution of the maximum of the two means, from which all your inference may be easily obtained. Commented Jul 25, 2021 at 11:31

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)$$.