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

    Доброго времени суток. В процессе разработки одного из методов кластеризации, возникла у меня потребность визуализировать гауссиан (нарисовать эллипс по сути) на плоскости по заданной ковариационной матрице. Но я как-то сразу и не задумался, что за простой отрисовкой обычного эллипса по 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), 
    
    


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



    Почитать




    Код можно найти у меня на гитхабе.
    Поделиться публикацией

    Комментарии 18

      +10
      Немного офтоп, но было бы неплохо, если бы авторы постов в разделе Математика менее грубо рисовали бы формулы. Вот например у вас



      а можно было бы нарисовать как



      Кстати, инлайн формулы тоже можно вставлять картинками, т.е. вместо Λ^(1/2) можно было написать что, согласитесь, выглядит намного лучше.

      Я уже давно просил администрацию Хабра добавить поддержку MathJax для этого блога — надеюсь что рано или поздно она появится.
        +2
        ага я тоже отсылал запрос на поддержку латеха

        так я что не понял то, вас размер не устраивает? или качество картинки? я для генерации картинок юзаю www.forkosh.com/mimetextutorial.html с размером шрифта \Large

        на счет инлайна я в курсе, иногда это нарушает межстрочный интервал, но согласен, что для мелких формул типа image это лучше
          0
          Качество картинки. Посмотрите как убого рисуются прямые линии на U и Λ, например. Мне кажется стоит использовать механизм, который сохраняет формулы в хорошем качестве. Я посмотрел на ссылку что вы дали выше и там действительно плохо все рендерится. Я бы не пользовался таким сервисом.

          А нарушение интервала — вещь некритичная, особенно если не пытаться туда воткнуть что-нибудь с \displaystyle.
            0
            спс за отзыв, нада будет потестить другие сервисы
          +1
          абрабрабрабрабрабрабрабрабрабр
          абрабрабрабрабрабрабрабрабрабр
          абрабрабрабраimageбрабрабрабраб
          абрабрабрабрабрабрабрабрабрабр
          абрабрабрабрабрабрабрабрабрабр

          вот например, можно варировать с размером картинки конечно, но думаю не просто будет не нарушить интервал, оставить картинку качественной и читаемой, тут палка о двух концах, либо так либо так
            +1
            можно поиграться с юникодом ещё,
            абрабрабрабрабрабрабрабрабрабр
            абрабрабрабрабрабрабрабрабрабр
            абрабрабрабра Λ½ или Λ¹´² брабр
            абрабрабрабрабрабрабрабрабрабр
            абрабрабрабрабрабрабрабрабрабр
        • НЛО прилетело и опубликовало эту надпись здесь
            0
            Элипс показывает в какой области находится двумерная случайная величина, в соответствии с заданной ошибкой первого рода.
            • НЛО прилетело и опубликовало эту надпись здесь
                0
                Нет, с заданной ошибкой первого рода вполне верно. Или с доверительным интервалом. Когда говорят, что велина равна 7±3 с вероятностью 95% подразумевают доверительный интервал. Ошибка первого рода равна 100 — 95 = 5%.
                • НЛО прилетело и опубликовало эту надпись здесь
                    0
                    «5% — это же вероятность того, что мы совершим ошибку первого рода.»
                    А, вы об этом. Ну да, конечно, вероятность ошибки.
              +6
              В данный период времени, одна из проблем которой я занимаюсь — это кластеризация с использованием расстояния Махаланобиса. Планирую кстати пост запилить как будут какие то интересные результаты. Если использовать EM алгоритм, как в k-means, то слишком быстро алгоритм впадает в локальный минимум, и результат на много хуже чем в k-means.



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



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

              Возможно понимание природы этих эллипсоидов поможет придумать новые эвристики.
              • НЛО прилетело и опубликовало эту надпись здесь
                  0
                  Вы же используете ЕМ для смеси гауссовских распределений? Картинка очень похожа на такую смесь, я немного удивлен, что ЕМ так сошелся. Если Вам не трудно, скажите, какие Вы данные используете, очень хочу посмотреть, откуда такой эффект.
                    0
                    полностью искусственные данные, сгенерил несколько 2д гауссиан, черным отмечена инициализация для каминса
                      +1
                      А все-таки при чем здесь каминс, если Вы кластеризуете смесь гауссиан? Вот пример с ЕМ алгоритмом.

                      Было:
                      Стало:

                      Плюс, есть оценки матожиданий и матриц ковариации для всех гауссиан.

                      Код (используются пакеты mclust и mvtnorm)
                      set.seed(1)
                      
                      x1<-rmvnorm(1000,c(0,0),matrix(c((40/3)^2,0,0,5^2),nrow=2))
                      x2<-rmvnorm(1000,c(0,-30),matrix(c(6^2,0,0,6^2),nrow=2))
                      x3<-rmvnorm(1000,c(-33,-25),matrix(c(6^2,-34,-34,6^2),nrow=2))
                      x4<-rmvnorm(1000,c(43,0),matrix(c(6^2,-34,-34,6^2),nrow=2))
                      
                      x <- rbind(x1,x2,x3,x4)
                      
                      plot(x,xlim=c(-60,60),ylim=c(-45,20),pch=16)
                      
                      emcl <- Mclust(x,G=4)
                      cls <- unlist(apply(emcl$z,1,which.max))
                      
                      plot(x,xlim=c(-60,60),ylim=c(-45,20),col=cls,pch=16)
                      
                      

            Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

            Самое читаемое