Примечание

Перейти к концу, чтобы скачать полный код примера или запустить этот пример в браузере через JupyterLite или Binder.

Робастная оценка ковариации и значимость расстояний Махаланобиса

В этом примере показана оценка ковариации с использованием расстояний Махаланобиса для данных с гауссовым распределением.

Для данных с гауссовым распределением расстояние от наблюдения \(x_i\) до моды распределения можно вычислить с помощью расстояния Махаланобиса:

\[d_{(\mu,\Sigma)}(x_i)^2 = (x_i - \mu)^T\Sigma^{-1}(x_i - \mu)\]

где \(\mu\) и \(\Sigma\) — это положение и ковариация исходных гауссовых распределений.

На практике \(\mu\) и \(\Sigma\) заменяют некоторыми оценками. Стандартная оценка ковариации методом максимального правдоподобия (MLE) очень чувствительна к наличию выбросов в наборе данных, а значит, чувствительны к ним и вычисляемые далее расстояния Махаланобиса. Лучше использовать робастную оценку ковариации, чтобы гарантировать устойчивость оценки к «ошибочным» наблюдениям в наборе данных и обеспечить точное отражение вычисленными расстояниями Махаланобиса истинной структуры наблюдений.

Оценка минимального определителя ковариации (MCD) — это робастная оценка ковариации с высокой точкой разрушения (то есть её можно использовать для оценки матрицы ковариации сильно загрязнённых наборов данных, содержащих до \(\frac{n_\text{samples}-n_\text{features}-1}{2}\) выбросов). Идея MCD состоит в том, чтобы найти \(\frac{n_\text{samples}+n_\text{features}+1}{2}\) наблюдений, эмпирическая ковариация которых имеет наименьший определитель. Таким образом получается «чистое» подмножество наблюдений, на основе которого вычисляются стандартные оценки положения и ковариации. MCD был предложен P. J. Rousseeuw в работе [1].

В этом примере показано, как выбросы влияют на расстояния Махаланобиса. При использовании расстояний Махаланобиса, основанных на стандартной оценке ковариации MLE, наблюдения из загрязняющего распределения невозможно отличить от наблюдений из исходного гауссова распределения. При использовании расстояний Махаланобиса на основе MCD эти две группы становятся различимы. Такие методы применяются, в частности, для обнаружения выбросов, ранжирования наблюдений и кластеризации.

Примечание

См. также Робастная и эмпирическая оценки ковариации

Литература

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

Генерация данных

Сначала создадим набор данных из 125 образцов с двумя признаками. Оба признака имеют гауссово распределение со средним значением 0, но стандартное отклонение первого признака равно 2, а второго — 1. Затем 25 образцов заменяются гауссовыми выбросами, у которых стандартное отклонение первого признака равно 1, а второго — 7.

import numpy as np

# for consistent results
np.random.seed(7)

n_samples = 125
n_outliers = 25
n_features = 2

# generate Gaussian data of shape (125, 2)
gen_cov = np.eye(n_features)
gen_cov[0, 0] = 2.0
X = np.dot(np.random.randn(n_samples, n_features), gen_cov)
# add some outliers
outliers_cov = np.eye(n_features)
outliers_cov[np.arange(1, n_features), np.arange(1, n_features)] = 7.0
X[-n_outliers:] = np.dot(np.random.randn(n_outliers, n_features), outliers_cov)

Сравнение результатов

Ниже мы применяем к нашим данным оценки ковариации на основе MCD и MLE и выводим оценённые матрицы ковариации. Обратите внимание, что оценённая дисперсия второго признака значительно выше при использовании оценки на основе MLE (7.5), чем при использовании робастной оценки MCD (1.2). Это показывает, что робастная оценка на основе MCD гораздо устойчивее к выбросам, которые были созданы так, чтобы иметь значительно большую дисперсию второго признака.

import matplotlib.pyplot as plt

from sklearn.covariance import EmpiricalCovariance, MinCovDet

# fit an MCD robust estimator to data
robust_cov = MinCovDet().fit(X)
# fit an MLE estimator to data
emp_cov = EmpiricalCovariance().fit(X)
print(
    "Estimated covariance matrix:\nMCD (Robust):\n{}\nMLE:\n{}".format(
        robust_cov.covariance_, emp_cov.covariance_
    )
)

Estimated covariance matrix:
MCD (Robust):
[[ 3.60075119 -0.07640781]
 [-0.07640781  1.51855963]]
MLE:
[[ 3.23773583 -0.24640578]
 [-0.24640578  7.51963999]]

Чтобы нагляднее показать различия, построим контуры расстояний Махаланобиса, вычисленных обоими методами. Обратите внимание, что робастные расстояния Махаланобиса на основе MCD гораздо лучше соответствуют чёрным точкам без выбросов, тогда как расстояния на основе MLE сильнее зависят от красных точек-выбросов.

import matplotlib.lines as mlines

fig, ax = plt.subplots(figsize=(10, 5))
# Plot data set
inlier_plot = ax.scatter(X[:, 0], X[:, 1], color="black", label="inliers")
outlier_plot = ax.scatter(
    X[:, 0][-n_outliers:], X[:, 1][-n_outliers:], color="red", label="outliers"
)
ax.set_xlim(ax.get_xlim()[0], 10.0)
ax.set_title("Mahalanobis distances of a contaminated data set")

# Create meshgrid of feature 1 and feature 2 values
xx, yy = np.meshgrid(
    np.linspace(plt.xlim()[0], plt.xlim()[1], 100),
    np.linspace(plt.ylim()[0], plt.ylim()[1], 100),
)
zz = np.c_[xx.ravel(), yy.ravel()]
# Calculate the MLE based Mahalanobis distances of the meshgrid
mahal_emp_cov = emp_cov.mahalanobis(zz)
mahal_emp_cov = mahal_emp_cov.reshape(xx.shape)
emp_cov_contour = plt.contour(
    xx, yy, np.sqrt(mahal_emp_cov), cmap=plt.cm.PuBu_r, linestyles="dashed"
)
# Calculate the MCD based Mahalanobis distances
mahal_robust_cov = robust_cov.mahalanobis(zz)
mahal_robust_cov = mahal_robust_cov.reshape(xx.shape)
robust_contour = ax.contour(
    xx, yy, np.sqrt(mahal_robust_cov), cmap=plt.cm.YlOrBr_r, linestyles="dotted"
)

# Add legend
ax.legend(
    [
        mlines.Line2D([], [], color="tab:blue", linestyle="dashed"),
        mlines.Line2D([], [], color="tab:orange", linestyle="dotted"),
        inlier_plot,
        outlier_plot,
    ],
    ["MLE dist", "MCD dist", "inliers", "outliers"],
    loc="upper right",
    borderaxespad=0,
)

plt.show()

Mahalanobis distances of a contaminated data set

Наконец, покажем, как расстояния Махаланобиса на основе MCD позволяют различать выбросы. Извлечём кубический корень из расстояний Махаланобиса, получив приблизительно нормальные распределения (как предложили Wilson и Hilferty [2]), а затем построим диаграммы размаха для значений образцов без выбросов и образцов-выбросов. При использовании робастных расстояний Махаланобиса на основе MCD распределения образцов-выбросов и образцов без выбросов разделены сильнее.

fig, (ax1, ax2) = plt.subplots(1, 2)
plt.subplots_adjust(wspace=0.6)

# Calculate cubic root of MLE Mahalanobis distances for samples
emp_mahal = emp_cov.mahalanobis(X - np.mean(X, 0)) ** (0.33)
# Plot boxplots
ax1.boxplot([emp_mahal[:-n_outliers], emp_mahal[-n_outliers:]], widths=0.25)
# Plot individual samples
ax1.plot(
    np.full(n_samples - n_outliers, 1.26),
    emp_mahal[:-n_outliers],
    "+k",
    markeredgewidth=1,
)
ax1.plot(np.full(n_outliers, 2.26), emp_mahal[-n_outliers:], "+k", markeredgewidth=1)
ax1.axes.set_xticklabels(("inliers", "outliers"), size=15)
ax1.set_ylabel(r"$\sqrt[3]{\rm{(Mahal. dist.)}}$", size=16)
ax1.set_title("Using non-robust estimates\n(Maximum Likelihood)")

# Calculate cubic root of MCD Mahalanobis distances for samples
robust_mahal = robust_cov.mahalanobis(X - robust_cov.location_) ** (0.33)
# Plot boxplots
ax2.boxplot([robust_mahal[:-n_outliers], robust_mahal[-n_outliers:]], widths=0.25)
# Plot individual samples
ax2.plot(
    np.full(n_samples - n_outliers, 1.26),
    robust_mahal[:-n_outliers],
    "+k",
    markeredgewidth=1,
)
ax2.plot(np.full(n_outliers, 2.26), robust_mahal[-n_outliers:], "+k", markeredgewidth=1)
ax2.axes.set_xticklabels(("inliers", "outliers"), size=15)
ax2.set_ylabel(r"$\sqrt[3]{\rm{(Mahal. dist.)}}$", size=16)
ax2.set_title("Using robust estimates\n(Minimum Covariance Determinant)")

plt.show()

Using non-robust estimates (Maximum Likelihood), Using robust estimates (Minimum Covariance Determinant)

Общее время выполнения скрипта: (0 минут 0.252 секунды)

Робастная оценка ковариации и значимость расстояний Махаланобиса

Литература

Генерация данных

Сравнение результатов

Похожие примеры