Несоответствие между расчетом устойчивой ковариационной матрицы Python и R
В настоящее время я разрабатываю статистический пакет на Python, используя код R в качестве образца, и заметил разные результаты между двумя программами при расчете устойчивой ковариационной матрицы в Python и R.
Используя следующий код с эквивалентными входными данными (x) в Python:
# 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_
и в Р:
# 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]])
и в Р:
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
Возможно, я что-то упускаю, но не могу понять, в чем несоответствие. Сталкивались ли вы с чем-то подобным или, возможно, у вас есть какие-нибудь предложения для меня? Спасибо!
@ Грегор Томас, я тоже об этом думал, но единственные аргументы, которые, кажется, имеют значение, - это альфа в R и support_fraction в Python. К сожалению, установка этих равных ничего не меняет. Только Python позволяет указать, центрированы ли данные, и это также не может объяснить несоответствие. Я здесь в тупике...
Установка значения 1 для альфа и support_fraction делает выходные данные намного ближе (но все же не идентичными, хотя изменение начального числа, похоже, не меняет выходные данные ни для одного из них).
Я предполагаю, что вы в курсе, но напоминаю, что «надежная ковариационная матрица», используемая в методе определения минимальной ковариации, — это не то же самое, что «надежная ковариационная матрица», используемая, например, в методе определения минимальной ковариации. «надежные стандартные ошибки» и не обязательно гарантируют какую-либо последовательность, несмещенность или что-то еще, что вы получаете от некоторых других оценщиков.
Ответы 1
Программное обеспечение не идентично, поэтому простая передача данных с аргументами по умолчанию вряд ли приведет к идентичным результатам. Функция R имеет обширную документацию по своим аргументам и значениям по умолчанию, но функция Python предоставляет гораздо меньше подробностей, поэтому вам, возможно, придется обратиться к ее исходному коду, чтобы быть абсолютно уверенным в причине. Однако параметры, влияющие на оценку, включают как минимум следующее, в зависимости от используемого алгоритма R и ваших данных (т. е. если при решении достигается максимальное количество шагов):
Питон
- Предположим_центрированный
- support_fraction
- случайное_состояние
Р
- альфа
- нсамп
- Нмини
- масштабирование
- Максимальное количество шагов
- initHsets
- семя
- tolSolve
- использование.коррекция
- wgtFUN
Как вы можете видеть, они не совсем четко сопоставляются друг с другом, и это особенно проблематично, поскольку в документации Python отсутствуют уравнения, необходимые для прямого сравнения методов, например. определение размера подмножества (по крайней мере, без чтения нескольких статей).
Я бы порекомендовал точно определить, когда возникает расхождение. Вы сказали, что дело не в центрировании (хотя я бы проверил это вручную с предварительно центрированными данными). В случае с MCD, по-видимому, существует несколько возможных точек расхождения:
Отбор выборки: метод определения того, какие наблюдения используются.
Взвешивание: как взвешиваются эти наблюдения
Поправки: конечная выборка и метод коррекции согласованности.
Алгоритмы: инструкции, передаваемые решателю, метод решения проблемы (например, «Fast MCD» или «DetMcd» в R).
Проще всего начать с проверки использованных наблюдений. Это best в R и похоже, что это support в Python. Если использовались разные наблюдения, то нужно посмотреть количество и размер выборок, проверенных на наличие различий. Если бы использовались те же наблюдения, я бы посмотрел на необработанные оценки перед повторным взвешиванием и коррекцией. Python содержит это (с помощью raw_location_ и raw_covariance), и вы можете либо напрямую извлечь это в R (с помощью raw.center и raw.cov). Значения R для них могут быть пост-коррекцией (я в этом сомневаюсь), поэтому вы также можете использовать raw.cnp2, чтобы «исправить» оценки, если это необходимо. Если оценки совпадают, значит, дело в взвешивании или поправках.
Спасибо! Глядя на best в R и support_ в Python, я получаю совершенно другой список использованных наблюдений, так что это, должно быть, проблема. Что странно, так это то, что установка доли наблюдений с помощью support_fraction в Python на самом деле не контролирует, какая пропорция используется, похоже, что все сводится к случайному состоянию.
Проанализировали ли вы аргументы для каждой реализации, чтобы увидеть, совпадают ли они? Часто подобные расхождения возникают из-за разных настроек по умолчанию, например, когда предполагается, что данные уже центрированы или нормализованы, а другой нет.