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

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

Данные и их взаимное расположение

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



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

plot(data[, 1], data[, 2], pch=19, asp=1, 
     col=rgb(0, 0.5, 1, 0.2),
     xlab="x", ylab="y")

sigma <- cov(data)
m.x <- mean(data[, 1])
m.y <- mean(data[, 2])

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



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

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

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



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

Эллипс и распределение хи-квадрат

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

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

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

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

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

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

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

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

v <- qchisq(q, k)

Получение контура эллипса

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

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

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

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

e <- eigen(sigma)

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

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

Давайте теперь напишем функцию, которая рисует эллипс, со следующими входными данными:

• m.x, m.y — координаты центра масс


• sigma — ковариационная матрица


• q — квантиль распределения хи-квадрат


• n — плотность дискретизации эллипса (количество точек по которым будет строится эллипс)

GetEllipsePoints <- function(m.x, m.y, sigma, q = 0.75, n = 100)
{
  k <- qchisq(q, 2) # вычисляем значение квантиля
  sigma <- k * sigma # масштабирование ковариационной матрицы на значение квантиля
  e <- eigen(sigma) # вычисление собственных значений масштабированной ковариационной матрицы
  angles <- seq(0, 2*pi, length.out=n) # разбиваем круг на n-ое количество углов
  cir1.points <- rbind(cos(angles), sin(angles)) # генерируем точки на единичной окружности
  ellipse.centered <- (e$vectors %*% diag(sqrt(abs(e$values)))) %*% cir1.points # масштабируем и поворачиваем полученный датасет
  ellipse.biased <- ellipse.centered + c(m.x, m.y) # смещаем его до центра масс
  return(ellipse.biased) # готово
}

Результат

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

points(m.x, m.y, pch=20, col="yellow")

q <- seq(0.1, 0.95, length.out=10)
palette <- cm.colors(length(q))
for(i in 1:length(q))
{
  p <- GetEllipsePoints(m.x, m.y, sigma, q = q[i])
  points(p[1, ], p[2, ], type="l", col=palette[i])
}

e <- eigen(sigma)
v <- (e$vectors %*% diag(sqrt(abs(e$values))))
arrows(c(m.x, m.x), c(m.y, m.y), 
       c(v[1, 1] + m.x, v[1, 2] + m.x), c(v[2, 1] + m.y, v[2, 2] + m.y), 

В итоге получаем такую картину:

Почитать

Код можно найти у меня на гитхабе.

This entry passed through the Full-Text RSS service — if this is your content and you're reading it on someone else's site, please read the FAQ at fivefilters.org/content-only/faq.php#publishers. Five Filters recommends:

• Massacres That Matter - Part 1 - 'Responsibility To Protect' In Egypt, Libya And Syria


• Massacres That Matter - Part 2 - The Media Response On Egypt, Libya And Syria


• National demonstration: No attack on Syria - Saturday 31 August, 12 noon, Temple Place, London, UK