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:
@glemaitre @thomasjpfan thanks for picking this up.

I have personally spent quite some time on this, and I think I have found the issue.

First the quick review of what MCD is doing:
(a) it tries to give robust estimation of covariance which is not sensitive to outlier
(b) in order to achieve (a), it labels some data point as outlier, and then calculate covariance on non-outliers
(c) to detect the outilers in (b), it just say data points with too large Mahalanobis distance are outliers (similar idea of in 1-dimension you can say data points with too extreme zscore are outliers)
(d) however Mahalanobis distance need covariance as input, which is sensitive to outliers
(e) it tries to select a subset of original data points, and make sure the subset has minimum determinant, and then use the subset to calculate covariance, and then calculate Mahalanobis distance for all data points, which is consider as a robust Mahalanobis data points
(f) then data points with too large robust Mahalanobis in (e) are "deleted" since they are considered to be outliers, this step is called as a "reweight" step in the original paper (step 7 on page 218)
(g) also need to adjust covariance estimation by "consistency" factor, the idea is you have removed some data points far away from center since you think they are outliers, but even if the data are drawn from normal distribution you also have chance to remove some data points near the tail, thus you will underestimate the covariance, and then you need chi-sqaured distribution to rescale the covariance (step 6 on page 218)

see the snapshot of original paper (A Fast Algorithm for the Minimum Covariance Determinant Estimator)

![image](https://user-images.githubusercontent.com/51205365/181595492-e68e698b-35da-4af3-946f-f89a6f4ef559.png)

I have checked the implementation of sklearn: 
the step 6, consistency factor implemented as method `correct_covariance`
https://github.com/scikit-learn/scikit-learn/blob/e5c7fb1e1f8d947b380f30da0c251b2e955ef725/sklearn/covariance/_robust_covariance.py#L779
the step7, reweight, implemented as method `reweight_covariance`
https://github.com/scikit-learn/scikit-learn/blob/e5c7fb1e1f8d947b380f30da0c251b2e955ef725/sklearn/covariance/_robust_covariance.py#L819

Now about the issue:
We use consistency factor to correct the under estimation when calculate robust Mahalanobis Distance, since we only take a subsample, this is correct.
https://github.com/scikit-learn/scikit-learn/blob/e5c7fb1e1f8d947b380f30da0c251b2e955ef725/sklearn/covariance/_robust_covariance.py#L772-L775
However, after the reweighting step, we have deleted some data points with too large Mahalanbis Distance, we would need another round of consistency factor correction.

For reference, I found another implementation of MCD, in the "Details" part they mentioned they did have two rounds of consistency correction.
https://rdrr.io/cran/robustbase/man/covMcd.html

Labels: module:covariance, Needs Investigation