Визуализация двумерного гауссиана на плоскости

Доброго времени суток. В процессе разработки одного из методов кластеризации, возникла у меня потребность визуализировать гауссиан (нарисовать эллипс по сути) на плоскости по заданной ковариационной матрице. Но я как-то сразу и не задумался, что за простой отрисовкой обычного эллипса по 4 числам скрываются какие то трудности. Оказалось, что при расчете точек эллипса используются собственные числа и собственные векторы ковариационной матрицы, расстояние Махаланобиса, а так же квантили распределение хи-квадрат, которое я, честно говоря, не использовал со времен университета ни разу.

Давайте начнем с начальных условий. Итак, мы имеем некоторый массив двумерных данных

для которого мы можем легко узнать ковариационную матрицу и средние значения (центр будущего эллипса):

Прежде чем приступить к отрисовке эллипса, нужно определиться с тем, какого размера будет фигура. Вот несколько примеров:

Для определения размера эллипса вспомним расстояние Махаланобиса между двумя случайными векторами из одного вероятностного распределения с ковариационной матрицей Σ:

Так же можно определить расстояние от случайного вектора x до множества со средним значением μ и ковариационной матрицей Σ:

Стоит заметить, что в случае, когда Σ равна единичной матрице, расстояние Махаланобиса вырождается в Евклидово расстояние. Смысл расстояния Махаланобиса в том, что оно учитывает корреляцию между переменными; или другими словами, учитывается разброс данных относительно центра масс (предполагается, что разброс имеет форму эллипсоида). В случае же использования Евклидова расстояния, используется предположение, что данные распределены сферически (равномерно по всем измерениям) вокруг центра масс. Проиллюстрируем это следующим графиком:

Желтым цветом отмечен центр масс, а две красные точки из набора данных, расположенные на главных осях эллипса, находятся на одинаковом, в смысле Махаланобиса, расстоянии от центра масс.

Эллипс является центральной невырожденной кривой второго порядка, его уравнение можно записать в общем виде (ограничения выписывать не будем):

С другой стороны, можно записать уравнение эллипса в матричной форме (в однородных двумерных координатах):

и получить следующее выражение, для того, чтобы показать, что в матричной форме действительно задан эллипс:

Теперь вспомним расстояние Махаланобиса, и рассмотрим его квадрат:

Легко заметить, что это представление идентично записи уравнения эллипса в матричной форме. Таким образом, мы убедились, что расстояние Махаланобиса описывает эллипс в Евклидовом пространстве. Расстояние Махаланобиса — это просто расстояние между заданной точкой и центром масс, делённое на ширину эллипсоида в направлении заданной точки.

Наступает тонкий момент для понимания, у меня это заняло некоторое время, что бы осознать: квадрат расстояния Махаланобиса это сумма квадратов k-ого количества нормально распределенных случайных величин, где n — это размерность пространства.

Вспомним, что такое распределение хи-квадрат — это распределение суммы квадратов k независимых стандартных нормальных случайных величин (это распределение параметризируется количеством степеней свободы k). А это как раз и есть расстояние Махаланобиса. Таким образом вероятность того, что x находится внутри эллипса выражается следующей формулой:

И вот мы пришли к ответу на вопрос о размере эллипса — его размер мы будем детерминировать квантилями распределения хи-квадрат, это легко делается в R (где q из (0, 1) и k — количество степеней свободы):

Идея генерации контура нужного нам эллипса очень проста, мы просто возьмем ряд точек на единичной окружности, сместим эту окружность в центр масс массива данных, затем масштабируем и растянем эту окружность в нужных направлениях. Рассмотрим геометрическую интерпретацию многомерного гауссового распределения: как мы знаем, это эллипсоид, у которого направление главных осей задано собственными значениями ковариационной матрицы, а относительная длина главных осей задана корнем из соответствующих собственных значений.

На следующем графике изображены собственные векторы ковариационной матрицы масштабированные на корни из соответствующих собственных значений, направления соответствуют главным осям эллипса:

Рассмотрим разложение ковариационной матрицы в следующем виде:

где U — матрица образованная единичными собственными векторами матрицы Σ, а Λ — диагональная матрица, составленная из соответствующих собственных значений.

А также рассмотрим следующее выражение для случайного вектора из многомерного нормального распределения:

Таким образом, распределение N(μ, Σ) — это, по сути, стандартное многомерное нормальное распределение N(0, I) масштабированное на Λ^(1/2), повернутое на U и смещенное на μ.

• m.x, m.y — координаты центра масс
• sigma — ковариационная матрица
• q — квантиль распределения хи-квадрат
• n — плотность дискретизации эллипса (количество точек по которым будет строится эллипс)

Следующий код рисует множество доверительных эллипсов вокруг центра масс датасета: