Основной из проблем кластерного анализа практических данных является наличие выбросов. Большинство существующих методов кластеризации не учитывают их существование, из-за этого явно аномальные наблюдения включаются в состав каких-то кластеров, что может серьезно смещать их центры и влиять на качество классификации. Разумеется, можно сначала исходные данные проверить на выбросы, их отсеять и т.д., но тогда задача превратиться в двухэтапную, а хотелось бы, чтобы было "все и сразу".
Один из подходов к решению данной задачи (чтобы метод кластеризации автоматически отсеивал выбросы) получил название "optimally tuned robust improper maximum likelihood estimator" и был описан вот в этой статье 2017 года (http://dx.doi.org/10.1080/01621459.2015.1100996), а недавно и получил реализацию на R. Поговорим о нем.
Минутка теории
Если совсем вкратце - кластеризация с помощью EM - алгоритмов основана на предположении, что конкретные кластеры имеют форму, близкую к форме многомерного нормального распределения. Тогда задачу кластеризации можно рассматривать как задачу выделения компонентов смеси и кластеры определяются согласно вероятности попадания в тот или иной компонент.
Математический апгрейд классического алгоритма состоит в том, что берется чуть иное предположение. Считается, что функция вероятности распределения точек является суммой смеси гауссиан и равномерного распределения шумовых точек (считается, что выбросы распределены равномерно), то есть
- несобственная постоянная плотность (плотность
шума)
Это приводит к некоторому усложнению решения задачи оптимизации, которая теперь выглядит как решение задачи на максимум
но вот с такими ограничениями
Эта задача оптимизации (особенно с условиями, накладываемыми на отношение собственных значений) далеко не всегда имеет решение. Авторы честно об этом говорят, но также они утверждают, что их алгоритм достаточно спокойно справляется к кластерами довольно сложной формы. Проверим
Эксперимент
Создадим какой-нибудь двумерный датасет, в котором форма кластеров будет не особо отличаться от нормальных.
library(ggplot2)# Создаем базу данныхset.seed(739)n <- 500 # numer of pointsy <- abs(c(rnorm(n,5,2)))x <- c(rnorm(n,5,1))/sqrt(y)class<-rep(1,n)dat1 <- as.data.frame(cbind(x,y,class)) # Первый кластер - что-то типа кривоугольного прямоугольного треугольникаy <- c(rnorm(n,8,0.4))x <- rlnorm(n, meanlog = log(7), sdlog = 0.2) class<-rep(2,n)dat2 <- as.data.frame(cbind(x,y,class)) # Второй кластер - больше похож на форме на горизонтальный прямоугольникy <- c(runif(n, min = 2, max = 7))x <- c(rnorm(n,8,0.4))class<-rep(3,n)dat3 <- as.data.frame(cbind(x,y,class)) # Третий кластер - вертикальный прямоугольникy <- c(runif(n/10, min = 0, max = 10))x <- c(runif(n/10, min = 0, max = 10))class<-rep(0,n/10)dat4 <- as.data.frame(cbind(x,y,class)) # Шумdat <- rbind(dat1,dat2,dat3,dat4)colnames(dat) <- c("x","y", "class")dat$class <- as.factor(dat$class)dat_test <- as.data.frame(dat)p_base <- ggplot(dat,aes(x=x,y=y,color=class)) + geom_point()ggExtra::ggMarginal(p_base, groupColour = TRUE, groupFill = TRUE)
Получаем вот такую картинку
Здесь и далее меткой "0" обозначаются
наблюдения, определенные как шум
Пакет otrimle достаточно современен в том плане, что в нем уже есть функция, позволяющая сравнить результаты кластеризации в случае задания нескольких кластеров. Зададим число кластеров от 2 до 5, и попробуем посмотреть на ход выполнения.
library(otrimle)clus <- otrimleg(dat[,c(1,2)], G=2:5, monitor=1) # параметр monitor позволяет видеть ход выполнения
В результате выведется вот что-то такое
В теории, этого уже достаточно для определения оптимального числа кластеров (мы же решаем задачу на максимум, и, в теории, нужно взять число кластеров, для которого iloglik максимален. С другой же стороны, видно, что значимого преимущества в разбиении на 4 или 5 кластеров по сравнению с 3 нет - поэтому выбор количества кластеров в значительной степени эмпирика). Но желательно посмотреть еще и на количество наблюдений внутри кластеров. Поэтому, не отходя далеко, жмем
clus$solution
и получаем следующую информацию
Для каждого из построенных разбиений Noise - количество элементов, классифицированных как выбросы. Совсем аномальных кластеров (в которых 1,2,10 наблюдений) вроде бы нет - поэтому более подробно посмотрим на разбиение с 5 кластерами.
clus_2 <-otrimle(dat[,c(1,2)], G = 5, npr.max = 0.01, erc = 20, monitor = TRUE)# npr.max - максимальная доля выбросов в выборке# erc - соотношение максимального/минимального собственных значенияclus_2$codeclus_2$flag
Если clus_2$code возвращает 0, это значит, что задача оптимизации вообще не решилась, если clus_2$code = 1, это означает, что EM-алгоритм не сошелся (и надо поиграться с параметрами), если clus_2$code = 2, то все хорошо, алгоритм все посчитал.
Параметр clus_2$flag дает информацию о внутренней структуре работы алгоритма. Если
clus_2$flag = 1, то присутствует вырождение апостерирорных вероятностей
Фактически, это говорит о наличии лишних кластеров, так как вероятность одного из компонент смеси близка к 0.
clus_2$flag = 2, то для сходимости алгоритма было нарушено предположение о доли шума
clus_2$flag = 3, то для при расчетах ограничение на долю шума было использовано
clus_2$flag = 4, то для при расчетах было использовано ограничение на отношение собственных значений
В нашем случае все хорошо (clus_2$code = 2, clus_2$flag = 4), продолжаем разбираться дальше.
clus_2$mean # центра кластеровhead(clus_2$tau) # вероятности принадлежности к кластерамhead(clus_2$cluster) # принадлежность к кластерам
Тут мы можем посмотреть центра кластеров и вероятности отнесения наблюдений к шуму или иным кластерам. Теперь еще посмотрим на то, как алгоритм классифицировал наши наблюдения.
Слева - истинное разбиение, справа -
предложенное алгоритмом для случая 5 кластеров
Слева - истинное разбиение, справа -
предложенное алгоритмом для случая 4 кластеров
Слева - истинное разбиение, справа -
предложенное алгоритмом для случая 3 кластеров
Слева - истинное разбиение, справа -
предложенное алгоритмом для случая 2 кластеров
Можно заметить, что даже если угадать с числом кластеров (3), то результаты кластеризации, вообще-то, очень сильно будут отличаться от истинных - и более адекватную картину будут давать разбиения с большим числом кластеров, чем на самом деле. Это произошло из-за отличной от эллиптической формы кластеров, но в целом алгоритм работает неплохо даже в условиях шума.
