Python 和 R 计算鲁棒协方差矩阵之间的差异

Question

我目前正在使用 R 代码作为参考，在 Python 中开发一个统计包，并且我注意到在用 Python 和 R 计算鲁棒协方差矩阵时，两个程序之间存在不同的结果。

在 Python 中使用以下代码和等效输入数据 (

x

)：

# dummy data
x = np.array([
    [0.5488135, 0.71518937, 0.60276338, 0.54488318],
    [0.4236548, 0.64589411, 0.43758721, 0.891773],
    [0.96366276, 0.38344152, 0.79172504, 0.52889492],
    [0.56804456, 0.92559664, 0.07103606, 0.0871293],
    [0.0202184, 0.83261985, 0.77815675, 0.87001215],
    [0.97861834, 0.79915856, 0.46147936, 0.78052918],
    [0.11827443, 0.63992102, 0.14335329, 0.94466892],
    [0.52184832, 0.41466194, 0.26455561, 0.77423369],
    [0.45615033, 0.56843395, 0.0187898, 0.6176355],
    [0.61209572, 0.616934, 0.94374808, 0.6818203],
    [0.3595079, 0.43703195, 0.6976312, 0.06022547],
    [0.66676672, 0.67063787, 0.21038256, 0.1289263],
    [0.31542835, 0.36371077, 0.57019677, 0.43860151],
    [0.98837384, 0.10204481, 0.20887676, 0.16130952],
    [0.65310833, 0.2532916, 0.46631077, 0.24442559],
    [0.15896958, 0.11037514, 0.65632959, 0.13818295],
    [0.19658236, 0.36872517, 0.82099323, 0.09710128],
    [0.83794491, 0.09609841, 0.97645947, 0.4686512],
    [0.97676109, 0.60484552, 0.73926358, 0.03918779],
    [0.28280696, 0.12019656, 0.2961402, 0.11872772]
])

# fit mcd
mcd = MinCovDet().fit(x)

# define robust covariance variable
x_cov = mcd.covariance_

在 R 中：

# dummy data
x <- matrix(c(
  0.5488135, 0.71518937, 0.60276338, 0.54488318,
  0.4236548, 0.64589411, 0.43758721, 0.891773,
  0.96366276, 0.38344152, 0.79172504, 0.52889492,
  0.56804456, 0.92559664, 0.07103606, 0.0871293,
  0.0202184, 0.83261985, 0.77815675, 0.87001215,
  0.97861834, 0.79915856, 0.46147936, 0.78052918,
  0.11827443, 0.63992102, 0.14335329, 0.94466892,
  0.52184832, 0.41466194, 0.26455561, 0.77423369,
  0.45615033, 0.56843395, 0.0187898, 0.6176355,
  0.61209572, 0.616934, 0.94374808, 0.6818203,
  0.3595079, 0.43703195, 0.6976312, 0.06022547,
  0.66676672, 0.67063787, 0.21038256, 0.1289263,
  0.31542835, 0.36371077, 0.57019677, 0.43860151,
  0.98837384, 0.10204481, 0.20887676, 0.16130952,
  0.65310833, 0.2532916, 0.46631077, 0.24442559,
  0.15896958, 0.11037514, 0.65632959, 0.13818295,
  0.19658236, 0.36872517, 0.82099323, 0.09710128,
  0.83794491, 0.09609841, 0.97645947, 0.4686512,
  0.97676109, 0.60484552, 0.73926358, 0.03918779,
  0.28280696, 0.12019656, 0.2961402, 0.11872772
), nrow = 20, ncol = 4, byrow = TRUE)

# fit mcd
x.mcd <- covMcd(x)

# define robust covariance variable
x_cov <- x.mcd$cov

我得到的

x_cov

值非常不同。我知道这些算法都涉及随机性，但差异太大，无法归因于此。

例如，在Python中：

x_cov = array([[ 0.06669275,  0.01987514,  0.01294049,  0.0235569 ],
              [ 0.01987514,  0.0421388 , -0.00541365,  0.0462657 ],
              [ 0.01294049, -0.00541365,  0.06601437, -0.02931285],
              [ 0.0235569 ,  0.0462657 , -0.02931285,  0.08961389]])

在 R 中：

x_cov =             [,1]       [,2]        [,3]        [,4]
          [1,]  0.15762177 0.01044705 -0.04043184 -0.01187968
          [2,]  0.01044705 0.09957141  0.03036312  0.08703770
          [3,] -0.04043184 0.03036312  0.06045952 -0.01781794
          [4,] -0.01187968 0.08703770 -0.01781794  0.16634435

也许我在这里遗漏了一些东西，但我似乎无法找出其中的差异。您是否遇到过类似的情况或者您对我有什么建议？谢谢！

Answer 1

软件并不相同，因此简单地使用默认参数传递数据不太可能产生相同的结果。 R 函数对其参数和默认值有大量文档，但 Python 函数提供的详细信息要少得多，因此您可能必须查阅其源代码才能绝对确定原因。然而，影响估计的参数至少包括以下内容，具体取决于所使用的 R 算法和您的数据（即求解时是否达到最大步数）：

Python

• 假设_中心
• 支持分数
• 随机状态

R

• 阿尔法
• nsamp
• n迷你
• scalefn
• 最大c步数
• 初始化Hsets
• 种子
• tol求解
• 使用修正
• wgtFUN

正如您所看到的，这些并没有清楚地相互映射，这尤其成问题，因为 Python 文档省略了直接比较方法所需的方程，例如确定子集大小（至少在不阅读多篇文章的情况下）。

我的建议是精确定位何时出现分歧。 MCD 似乎存在一些可能的分歧点：

1. 样本选择：确定使用哪些观测值的方法

2. 加权：如何对这些观察值进行加权

3. 修正：有限样本和一致性修正法

4. 算法：传递给求解器的指令、解决问题的方法（例如 R 中的“Fast MCD”与“DetMcd”）

最简单的起点是检查所使用的观察结果。这是 R 中的

best

，看起来像是 Python 中的

support

。如果使用不同的观察结果，那么您需要查看检查差异的样本数量和大小。如果使用相同的观察结果，那么我会在重新加权和校正之前查看原始估计。 Python 包含这个（使用

raw_location_

和

raw_covariance

），并且您似乎可以直接在 R 中提取它（使用

raw.center

和

raw.cov

）。这些的 R 值可能是校正后的（我对此表示怀疑），因此如果需要，您还可以使用

raw.cnp2

来“取消校正”估计值。如果估计值匹配，则必须进行加权或修正。