Estimating means of multiple Gaussians

Question

Let's say we have two Gaussian distributions $\mathcal{N}(\mu_1, \sigma^2I_d)$ and $\mathcal{N}(\mu_2, \sigma^2I_d)$. We are trying to get estimators $\hat \mu_1, \hat \mu_2$ to minimize the following quantity, $$\mathbb{E}_{\hat \mu_1, \hat \mu_2}\left [\frac{\|\mu_1 - \hat \mu_1\|_2^2 + \|\mu_2 - \hat \mu_2\|_2^2}{2}\right],$$ i.e., the average $L_2$ distance to the actual means. If our estimators are unbiased, this is also the average variance of our estimators. Now assume we have $N$ samples from both distributions, i.e., total $2N$ samples. I am trying to understand the lowest we can get this error using these samples. For instance if I knew that $\mu_1\approx \mu_2$, I would perhaps use all the $2N$ samples and estimate $\frac{\mu_1 + \mu_2}{2}$. On the other hand, if I knew that the means are very far apart, the best thing to do (at least for unbiased estimators) is to use the respective MLE estimates for $\mu_1$ and $\mu_2$. I am wondering if there is a general answer to this question, which depends on the following quantity, $$\zeta_\star^2 := \left\|\frac{\mu_1 - \mu_2}{2}\right\|_2^2,$$, i.e., a theorem of alternative of sorts. More generally, when there are more than two distributions, say $M$, is there a standard way to prove variance lower bounds which depend on the following quantity, $$\zeta_\star^2 := \frac{1}{M}\sum_{m=1}^{M}\|\mu_m - \bar \mu\|_2^2,$$ where $\bar \mu = \frac{1}{M}\sum_{m=1}^{M}\mu_m$.

I looked into the standard technique for proving the lower bounds for Gaussian mean estimation (Le Cam's method and Cramer-Rao lower bounds), but I could not extend them to the case when we have more than one mean. I asked ChatGPT what to do, and it gave several wrong answers but at least suggested that Bhatacharya Coefficient might be relevant to answering this question.

I would appreciate it if you have any pointers. This is not a HW question. Also, this is my first post, so I apologize for missing any guidelines.

Edit 1

After reading a bit on James-Stein Estimator and following the comments on the original question, it makes sense to clarify my problem and ask a more rigorous question. We have samples of size $N$ from $M$ different Gaussian distributions, each with mean $\mu_m$ and a known co-variance $\sigma^2I_d$. We want to get estimators $\hat\mu_1, \dots, \hat \mu_M$ to minimize the average mean squared error (AMSE): $$AMSE(\hat\mu_1, \dots, \hat \mu_M) := \mathbb{E}\left[\frac{1}{M}\sum_{m\in[M]}\|\hat \mu_m - \mu_m\|_2^2\right],$$ where the expectation is w.r.t.,

1. sampling from these distributions, i.e., obtaining $X_{m,1}, \dots, X_{m,N}\sim \mathcal{N}(\mu, \sigma^2I_d)$ for all $m\in[M]$; and
2. potential randomness in our estimators $\hat\mu_1, \dots, \hat \mu_M$.

We know that the simple sample mean estimators $\hat\mu_m^{MLE} := \frac{1}{N}\sum_{i\in[N]}X_{m,i},\ \forall m\in[M]$ are inadmissible (i.e., are dominated by other estimators for some value of the means) when $M\geq 3$. This is due to the famous Stein's example. One can calculate the MSE of this estimator as follows, $$AMSE(\hat\mu_1^{MLE}, \dots, \hat\mu_M^{MLE}) = \frac{\sigma^2 d}{N}.$$ While inadmissible, this is optimal for unbiased estimators due to the Cramer-Rao Lower bound. Now, I am wondering if there are biased estimators that obtain better AMSE. In particular, consider the estimator which averages the samples across all the distributions, i.e., $$\hat\mu_1 = \dots = \hat \mu_M = \hat \mu^{ALL} = \frac{1}{MN}\sum_{m\in[M], i\in[N]} X_{m,i}.$$ This estimator has the following AMSE, $$AMSE(\hat \mu^{ALL}, \dots, \hat \mu^{ALL}) = \frac{\sigma^2d}{MN} + \frac{1}{M}\sum_{m=1}^{M}\|\mu_m - \bar \mu\|_2^2,$$ where $\bar \mu := \frac{1}{M}\sum_{m=1}^{M}\mu_m$. Now let's assume $\|\mu_m\|_2 = B$ for all $m\in[M]$, so we even know the scale of the true parameter (which is the key quantity we vary to check if an estimator is admissible). My question is whether there is an approach to study the smallest the AMSE can be as a function of $$\zeta_\star^2 := \frac{1}{M}\sum_{m=1}^{M}\|\mu_m - \bar \mu\|_2^2,$$ i.e., the "heterogeneity" of the distributions. One can see that when $\zeta_\star = 0$, the ideal thing is to use the estimator $\hat\mu^{ALL}$, while if $\zeta_\star$ is large, it "intuitively" seems like the best estimator is some admissible estimator for the problem such as Jason-Stein applied for the concatenated vector of means of dimension $dM$. Is it possible to say something intelligent in between these regimes? If I squint a little (a lot, perhaps), I can see a trade-off here akin to the bias-variance trade-off. That's what I am trying to understand.

Thanks for any pointers in advance!

Answer 1

If these random variables take values in $\mathbb R^n$ for some $n>1,$ then the thing you're calling $\sigma^2$ would be an $n\times n$ nonnegative-definite matrix, and usually that's not denoted in that way, so I will assume here that these are scalar-valued random variables and $\|\mu-\hat\mu\|_2^2$ is synonymous with $(\mu-\hat\mu)^2.$

If $n\ge3,$ then there are some biased estimators with uniformly smaller mean square errors than the MLE – a counterintuitive and even paradoxical result discovered in the 1950s. Google "James–Stein."

You say both distributions have the same variance $\sigma^2,$ and that makes things easier. If they didn't have equal variances then there would be a question of whether the variances are "known," so that the values of the estimators are allowed to depend upon them. And here I think that if only the ratio of variances were known, that would be as much as you would want to use. Then one would use a weighted average with weights proportional to the reciprocals of the variances.

Suppose one of the two samples is $X_{1,1},\ldots,X_{1,N}$ and the other is $X_{2,1},\ldots X_{2,N}.$

The conditional probability distribution of $(X_{1,1},\ldots,X_{1,N},X_{2,1},\ldots X_{2,N})$ given the pair of sums $ \left( X_{1,1} + \cdots + X_{1,N}, \, X_{2,1} + \cdots + X_{2,N} \right)$ does not depend on the pair $(\mu_1,\mu_2).$ That is what it means to say that that pair of sums is a sufficient statistic for this family of probability distributions indexed by $(\mu_1,\mu_2).$ "Intuitvely" (if you can pardon such a vulgar word) the pair of sums has all the information in the sample that is relevant to estimation of the two parameters.

There is no function $g(X_{1,1}+\cdots+X_{1,N},\, X_{2,1}+ \cdots+ X_{2,N})$ (forbidden to depend on $(\mu_1,\mu_2,\sigma^2)$) whose expected value is zero regardless of the value of the pair $(\mu_1,\mu_2)$ except that is equal to zero with probability $1.$ I.e. there is no nontrivial unbiased estimator of zero based on this sufficient statistic. That fact is experessed by saying that this pair of sums is a complete statistic. This is essentially the same as the one-to-one nature of the two-sided Laplace transform.

Dividing both sums by $N$ we get a pair of sample means, which is complete, sufficient, and unbiased.

A corollary to the Lehmann–Scheffé theorem is that a statistic that is complete, sufficient, and unbiased is the unique estimator having the smallest mean squared error among all unbiased estimators.

Whether one can do better with a biased estimator is more than I am prepared to address right now. In some cases it can be done. For example $\frac1{N+1}\sum_{i=1}^n \left( X_i - (X_{1,1} + \cdots + X_{1,N})/N\right)^2$ is a biased estimator of $\sigma^2$ that has a smaller mean squared error than does $\frac1{N-1}\sum_{i=1}^n \left( X_i - (X_{1,1} + \cdots + X_{1,N})/N\right)^2,$ which is the best unbiased estimator based on the first sample $(X_{1,1},\ldots,X_{1,N}).$