MinCovDet estimation of covariance with strong bias?

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/23161

<div type='discussions-op-text'>

<sup>Originally posted by **hongfei0224** April 20, 2022</sup>
I was playing with MinCovDet (Minimum Covariance Determinant) in sklearn:
[https://scikit-learn.org/stable/modules/generated/sklearn.covariance.MinCovDet.html#sklearn.covariance.MinCovDet](url)

I see the idea of MinCovDet is to make the covariance matrix estimation robust when there is outlier in the data.
To understand the behavior, I tried the multivariate normal data **without** outliers.

What I did:
- generate multivariate_normal samples with given covariance matrix (the true covariance matrix)
- use `np.cov` to estimate empirical covariance matrix
- use `MinCovDet` to estimate robust covariance matrix
- repeat above with same "true covariance matrix"
- check the distribution of empirical covariance matrix and robust covariance matrix

I was expecting both `np.cov` and `MinCovDet` estimation should be centered around the true covariance matrix, since there is no outlier in data.

but seems not the case.

![image](https://user-images.githubusercontent.com/51205365/164061364-1daa7614-5178-46ad-b3be-0ff2e3c38a42.png)
![image](https://user-images.githubusercontent.com/51205365/164061425-00155926-f523-45a3-b872-ed83660aba5e.png)

the code to replicate results

```python
import pandas as pd
import numpy as np
from sklearn.covariance import EmpiricalCovariance, MinCovDet
import matplotlib.pyplot as plt
from tqdm.autonotebook import tqdm

res = []
res2= []
a = 2
b = 5
covar = 3 

for i in tqdm(range(1000)):
    x=np.random.multivariate_normal(mean=[0,0],cov=[[a,covar],[covar,b]], size=250).T
    res.append(np.cov(x))
    mcd = MinCovDet().fit(x.T)
    res2.append(mcd.covariance_)

pd.Series([xx[0][0] for xx in res]).hist(bins=30, label='empirical',alpha=0.7)
pd.Series([xx[0][0] for xx in res2]).hist(bins=30, label='cvd',alpha=0.7)
plt.axvline(a, ls='--', color='black')
plt.title('variance')
plt.legend()
plt.show()

pd.Series([xx[0][1] for xx in res]).hist(bins=30, label='empirical',alpha=0.7)
pd.Series([xx[0][1] for xx in res2]).hist(bins=30, label='cvd',alpha=0.7)
plt.axvline(covar, ls='--', color='black')
plt.title('covariance')
plt.legend()
plt.show()
```

COMMENT:
I believe @hongfei0224 is correct in that the issue is in the corrections being applied. To test his proposed solution, I slightly adjusted the code, see below. Basically, based on my understanding of the source code, I applied the `correct_covariance` method again and created the same plots but comparing the corrected (co)variances to the ones calculated by `np.cov`.

```python
#!/usr/bin/env python3
import pandas as pd
import numpy as np
from sklearn.covariance import EmpiricalCovariance, MinCovDet
import matplotlib.pyplot as plt

res = []
res2 = []
res2_cor = []
a = 2
b = 5
covar = 3

for i in range(1000):
    x=np.random.multivariate_normal(mean=[0,0],cov=[[a,covar],[covar,b]], size=250).T

    res.append(np.cov(x))

    mcd = MinCovDet().fit(x.T)
    res2.append(mcd.covariance_)

    res2_cor.append(mcd.correct_covariance(x.T))

plt.figure()
pd.Series([xx[0][0] for xx in res]).hist(bins=30, label='empirical',alpha=0.7)
pd.Series([xx[0][0] for xx in res2]).hist(bins=30, label='cvd',alpha=0.7)
plt.axvline(a, ls='--', color='black')
plt.title('variance')
plt.legend()
plt.savefig('variance')

plt.figure()
pd.Series([xx[0][0] for xx in res]).hist(bins=30, label='empirical',alpha=0.7)
pd.Series([xx[0][0] for xx in res2_cor]).hist(bins=30, label='cvd',alpha=0.7)
plt.axvline(a, ls='--', color='black')
plt.title('variance (after correction)')
plt.legend()
plt.savefig('variance (after correction)')

plt.figure()
pd.Series([xx[0][1] for xx in res]).hist(bins=30, label='empirical',alpha=0.7)
pd.Series([xx[0][1] for xx in res2]).hist(bins=30, label='cvd',alpha=0.7)
plt.axvline(covar, ls='--', color='black')
plt.title('covariance')
plt.legend()
plt.savefig('covariance')

plt.figure()
pd.Series([xx[0][1] for xx in res]).hist(bins=30, label='empirical',alpha=0.7)
pd.Series([xx[0][1] for xx in res2_cor]).hist(bins=30, label='cvd',alpha=0.7)
plt.axvline(covar, ls='--', color='black')
plt.title('covariance (after correction)')
plt.legend()
plt.savefig('covariance (after correction)')
```

The result aren't as expected:
![image](https://github.com/scikit-learn/scikit-learn/assets/29434548/de3192ff-8f67-42cf-af54-24b217cd41de)

![image](https://github.com/scikit-learn/scikit-learn/assets/29434548/b15e221d-2e68-407d-9aee-e4c59e8f4370)

I believe the problem is in that there are two (well in fact, three, including a small sample correction) different consistency corrections to consider in the fast MCD algorithm by Rousseuw and Van Driessen.

Firstly, the MCD algorithm (even before making it fast) needs a correction to be Fisher consistent. This is the one that (I believe) is necessary to correct for removing outliers even though they could occur with a normal distribution, and it is the one that @hongfei0224 noted can be found in the linked R code:
```R
[.MCDcons](https://rdrr.io/cran/robustbase/man/covMcd.html) <- [function](https://rdrr.io/r/base/function.html)(p, alpha)
{
    qalpha <- [qchisq](https://rdrr.io/r/stats/Chisquare.html)(alpha, p)
    caI <- [pgamma](https://rdrr.io/r/stats/GammaDist.html)(qalpha/2, p/2 + 1) / alpha
    1/caI
}
```
Source: https://rdrr.io/cran/robustbase/src/R/covMcd.R

The correction is

$$
c_\alpha = \frac{\alpha}{F_{\Gamma(\frac{p}{2} + 1, 1)}(\frac{\chi_\alpha^2}{2})}
$$

Here $\alpha$ denotes the `support_fraction` in SKLearn's implementation, and $\chi_\alpha^2$ is the $\alpha$ th percentile of the Chi-squared distribution with $p$ (`n_features`) degrees of freedom, and $F_{\Gamma(\frac{p}{2} + 1, 1)}$ denotes the CDF of the Gamma distribution with parameters $p / 2$ and $1$.

I sadly haven't been able to find a source for this, and I'm not allowed to share my lecturer's slides; I'll contact him to see if I can find a source.

Additionally, there is a small sample correction $c_{np}$ which is determined empirically by Pison et al, 2002.
Lastly, there is the correction factor for the fast MCD, described in the paper by Rousseuw and Van Driessen in 1999.

If I'm understanding the source code right, only the last correction is implemented in SKLearn. I believe the way it should be, is that when calculating the raw MCD, both $c_\alpha$ (and $c_{np}$, although it is less important) should be applied. Then the correction from the paper on the fast MCD algorithm should be applied; I believe this is currently correct. Then, when calculating the reweighted MCD, $c_\alpha$ and $c_{np}$ should be applied again.

Labels: module:covariance, Needs Investigation