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:
Thanks @Encephala. I am now working on a PR to fix this issue. I think I figured out the issues, which are a combination of all what's been mentioned above. I would appreciate feedback.

**Background:**

First, I found a source for what @Encephala and @hongfei0224 are saying, in the [documentation of fastMcd](https://rdrr.io/cran/robustbase/f/inst/doc/fastMcd-kmini.pdf):

> ... the raw MCD estimate of scatter is their covariance matrix, multiplied by a consistency
factor `.MCDcons(p, h/n))` and (by default) a finite sample correction factor `.MCDcnp2(p, n, alpha)`,
to make it consistent at the normal model and unbiased at small samples.

The paper by [Pison 2002](https://link.springer.com/article/10.1007/s001840200191) explains it nicely, and cites [C. Croux & G. Haesbroeck 1999](https://www.sciencedirect.com/science/article/pii/S0047259X99918390?ref=cra_js_challenge&fr=RR-1) for the correction.

To understand what is happening, we need to consider that the MinCovDet algorithm has 2 steps that require correction, and that there are 2 different covariance matrices:

1) The **minimum determinant covariance** (MCD) is first estimated in `fast_mcd`, using a subset of the data. MCD is then used to compute the Mahalanobis distances of all the points in the datset (the variable `dist`), which is used to define outliers.
2) A new **robust covariance** is estimated without outliers in the method `reweight_covariance`. This step is called "reweighting".

Both steps 1) and 2) above first compute the empirical covariance from a subset of points, and in both cases this should be followed by two types of weighting:
**A)** consistency correction
**B)** finite sample correction.

**MinCovDet problems:**

There are different issues with the `MinCovDet` implementation:
1) Consistency correction A) is only applied after the **MCD** step, in the `correct_covariance` call, but not after the **robust covariance** step
2) The consistency correction that is applied is different from the one in the R version. The `scikit-learn` version uses the correction from the Fast-MCD 1999 paper, based on the median distance (which seems somewhat adhoc). R uses a correction from Croux & Haesbroeck.
3) Finite sample correction B) is not applied in either step.

Also, there is something misleading in the implementation. The method `correct_covariance` says that it returns the corrected "robust covariance". However, that's not correct, it returns the corrected **MCD** covariance, which is the `raw_covariance` computed earlier than the call to `reweight_covariance`.

**Solution:**

I was able to reduce the bias by solving point 1) above. I added The consistency correction for the **robust covariance**, as used by the R version, [in these lines](https://github.com/scikit-learn/scikit-learn/pull/32117/commits/cee4252f80821a2bd746a5944187a84bb9e62f66#diff-a5e40ac2b0d6abbfa016939827fc8b1b3747486549cc008129a56a2ee6231ef9R908). The new outputs look like this:

<img width="400" height="479" alt="Image" src="https://github.com/user-attachments/assets/659bb90a-e799-4802-bdce-856c91def6bd" />

<img width="400" height="479" alt="Image" src="https://github.com/user-attachments/assets/ab35ce67-7b01-428b-ac64-bcdf4853c391" />

There is still some bias, tending to underestimate the covariance elements.

**Action plan and new suggestions:**

I think it is best to submit a PR with the partial fix that reduces the bias. A future PR should add the finite sample correction.

I also suggest changing the first consistency correction to match that from the R implementation. The correction from the Fast-MCD 1999 paper seems adhoc, while according [to this paper](https://link.springer.com/article/10.1007/s001840200191) , the correction used by the R implementation "makes the MCD scatter estimator consistent at the normal model".

I also suggest changing the `correct_covariance` method to be less misleading. I suggest changing it to not return anything, since there does not seem to be any use for the corrected **MCD** matrix. Alternatively, make the docstrings clearer, saying that it returns the "Corrected Minimum Covariance Determinant estimate", rather than "Corrected robust covariance estimate" as is now.


**Future improvements:**

I think that solving point 3) and adding the Finite sample correction will remove the bias. Implementing this is not as straight forward though, I tried to follow the R implementation but it's pretty messy and it uses a lot of hard-coded values for the empirical corrections.

Any comments @hongfei0224 @Encephala @glemaitre?

Labels: module:covariance, Needs Investigation