- PVSM.RU - https://www.pvsm.ru -

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

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

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

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

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])

Прежде чем приступить к отрисовке эллипса, нужно определиться с тем, какого размера будет фигура. Вот несколько примеров:
Визуализация двумерного гауссиана на плоскости
Для определения размера эллипса вспомним расстояние Махаланобиса [4] между двумя случайными векторами из одного вероятностного распределения с ковариационной матрицей Σ:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

И вот мы пришли к ответу на вопрос о размере эллипса — его размер мы будем детерминировать квантилями [7]распределения хи-квадрат, это легко делается в 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), 

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

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

Почитать

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

Автор: mephistopheies

Источник [11]


Сайт-источник PVSM.RU: https://www.pvsm.ru

Путь до страницы источника: https://www.pvsm.ru/matematika/46701

Ссылки в тексте:

[1] гауссиан: http://en.wikipedia.org/wiki/Multivariate_normal_distribution

[2] ковариационной матрице: http://ru.wikipedia.org/wiki/%D0%9A%D0%BE%D0%B2%D0%B0%D1%80%D0%B8%D0%B0%D1%86%D0%B8%D0%BE%D0%BD%D0%BD%D0%B0%D1%8F_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D0%B0

[3] собственные числа и собственные векторы: http://ru.wikipedia.org/wiki/%D0%A1%D0%BE%D0%B1%D1%81%D1%82%D0%B2%D0%B5%D0%BD%D0%BD%D1%8B%D0%B9_%D0%B2%D0%B5%D0%BA%D1%82%D0%BE%D1%80

[4] расстояние Махаланобиса: http://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D1%81%D1%82%D0%BE%D1%8F%D0%BD%D0%B8%D0%B5_%D0%9C%D0%B0%D1%85%D0%B0%D0%BB%D0%B0%D0%BD%D0%BE%D0%B1%D0%B8%D1%81%D0%B0

[5] квантили: http://ru.wikipedia.org/wiki/%D0%9A%D0%B2%D0%B0%D0%BD%D1%82%D0%B8%D0%BB%D1%8C

[6] распределение хи-квадрат: http://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D1%85%D0%B8-%D0%BA%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82

[7] квантилями: http://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%BE%D1%86%D0%B5%D0%BD%D1%82%D0%B8%D0%BB%D0%B8_%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D1%8F_%D1%85%D0%B8-%D0%BA%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82

[8] Gaussians: http://www.cs.ubc.ca/~murphyk/Teaching/CS340-Fall06/reading/gaussians.pdf

[9] Some notes on ellipses: http://people.mmci.uni-saarland.de/~ihrke/software/ellipse.pdf

[10] у меня на гитхабе: https://github.com/mephistopheies/ml-r

[11] Источник: http://habrahabr.ru/post/199060/