# 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 ... May 4 '20 at 11:50
• How would you practically go about making an estimate using bootstrap or jackknife methods? May 4 '20 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. Apr 30 at 19:33

## 1 Answer

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'. May 8 '20 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. Jul 25 at 11:31