Say we have a multivariate normal with m dimensions (let's say that m=1,000), the mean vector $\mu$ is known only for the first 100 elements, but unknown for the remaining 900.

We have a single observation $X$ from this distribution. We are estimating the $\mu_i$ for each of the m dimensions using $X_i$, and we want to build a confidence interval for this estimation for each of the i elements.

The main diagonal of the covariance matrix can be assumed to have the same value $\sigma^2$ for all m elements, but like $\mu$, it is also unknown. The off-diagonal covariance is unknown, but we have some estimation of it $\hat \sigma_{i,j}$.

Now comes the trick - since we have the true $\mu_i$ values for the first E elements (let's say 100). If we had known that the off-diagonal covariance matrix is all 0, we could have just estimated $\sigma^2$ using the observations for which we have the true values using something like $\hat \sigma^2 = \sum_{i=1}^{100} \frac{(X_i-\mu_i)^2}{100}$, and apply that estimator for building the confidence interval for any $\mu_i$ when i>100.

But we have two issues here:

  1. our first 100 elements are not i.i.d (they are correlated) so it's not clear to me how to estimate $\sigma^2$
  2. when wanting to construct the confidence interval (CI) to some i (say, i=101), we know that it has a correlation with the other 100 observations we've seen - so it's not clear to me how to use that knowledge to adjust the CI.

What I describe sounds like some classic problem (feels related to James–Stein estimator), but I couldn't find any paper on it, and am not sure how it's called.

Any suggestions on how to solve it, or relevant references - will be great!

  • $\begingroup$ What I would do, depending the dimension of $\mu_{i}$, is maximum likelihood. Plug in the known $\mu$s, and estimate the correlation structure via the EM algorithm. Sure $m$ is big, but hopefully $n$ is bigger. Alternately, you can use a ridge type estimator for the means, and make the penalty for the "known" coefficients arbitrarily large, then estimate the covariance matrix as a nuissance parameter. Mean estimates are inefficient but not biased. $\endgroup$
    – AdamO
    Mar 12, 2021 at 16:04
  • $\begingroup$ Could you indicate in what way this is not a time series forecasting problem? $\endgroup$
    – whuber
    Mar 12, 2021 at 17:44

1 Answer 1


This only works if all the correlations are known and the observations from the data with known means are independent of the observations from the data with unknown means. Outside of that scenario, everything becomes more complicated and you would have to post the details of the specific scenario for any answers about that. I'm assuming you want two-sided symmetric 95% confidence intervals.

Suppose the covariance matrix for the first $100$ is $\sigma^2R$ where $R$ a known positive definite matrix with $1$'s on the diagonal and suppose the means, $\mu_i$, are known.

First, find a matrix $B$ such that $BRB'=I$ where $B'$ denotes the transpose of $B$ and $I$ denotes the identity matrix. $$Y=B\left((X_1,...,X_{100})'-(\mu_1,...,\mu_{100})\right)$$ is multivariate normal with mean $(0,...,0)$ and variance $\sigma^2I$.

Now, $\hat{\sigma}^2=\frac{1}{100}\sum_{i=1}^{100}Y_i^2$ is an unbiased estimate of $\sigma^2$. $\frac{1}{\sigma^2}\sum_{i=1}^{100}Y_i^2$ has a chi-square distribution with $100$ degrees of freedom.

For the remaining $900$ means $\mu_{101},...,\mu_{1000}$, you can construct the 900 confidence intervals in different ways depending on what statement you want to make about the confidence.

The confidence intervals $x_i \pm t_{100}(0.025) \hat{\sigma}$ each contain the true mean with probability 0.95. So, on average 45 of the 900 will not contain the true mean and the remaining 855 will contain the true mean. The per comparison error rate is 0.05.

Suppose you want to control the familywise error rate. All $900$ confidence intervals $x_i \pm q \hat{\sigma}$ must simultaneously contain the true mean with probability $0.95$. So, there will be a 5% chance that at least one of them does not contain the true mean. Let $\sigma^2 S$ denote the covariance matrix for the $900$ random variables where $S$ is a symmetric matrix with $1$'s on the diagonal. You want $$-q<\frac{X_{i}-\mu_{i}}{\hat{\sigma}}<q$$ to hold simultaneously for all $i=101,...,1000$ with probability 0.95. Note $Z=(\frac{X_{101}-\mu_{101}}{\sigma},...,\frac{X_{1000}-\mu_{1000}}{\sigma})'$ and $u=\frac{1}{\sigma^2}\sum_{i=1}^{100}Y_i^2=100\frac{\hat{\sigma}^2}{\sigma^2}$ are independent with $Z\sim N(0,S)$ and $u\sim \chi_{100}^2$. Also, $$\frac{Z}{\sqrt{u/100}}=(\frac{X_{101}-\mu_{101}}{\hat{\sigma}},...,\frac{X_{1000}-\mu_{1000}}{\hat{\sigma}})'$$ Therefore, $$(\frac{X_{101}-\mu_{101}}{\hat{\sigma}},...,\frac{X_{1000}-\mu_{1000}}{\hat{\sigma}})'$$ has a multivariate t-distribution. You can then use pmvt package in R for example to find $q$.
Here is an example assuming $S$ is the identity matrix. Note: it takes a long time to run using $m=1000$, I only checked that it works for smaller $m$. For the case where $S$ is the identity, there is an easier way to calculate $q$.


pmvt(lower=rep(-q,900),upper=rep(q,900),df=100) #check that it works

  • $\begingroup$ Thanks! The situation is that I do have many samples (X_1, X_2, etc...). So I can (theoretically) use these to estimate the correlation between them. The followup question that interests me is when I actually don't see X, but actually only see X+B (B is some bias vector.). But I'll need to think of how to frame it, and will then ask a followup question. $\endgroup$
    – Tal Galili
    Mar 12, 2021 at 17:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.