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

Один из подходов к решению данной задачи (чтобы метод кластеризации автоматически отсеивал выбросы) получил название "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 points
y <- 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)

Получаем вот такую картинку

Пакет 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$code
clus_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) # принадлежность к кластерам

Тут мы можем посмотреть центра кластеров и вероятности отнесения наблюдений к шуму или иным кластерам. Теперь еще посмотрим на то, как алгоритм классифицировал наши наблюдения.

Можно заметить, что даже если угадать с числом кластеров (3), то результаты кластеризации, вообще-то, очень сильно будут отличаться от истинных - и более адекватную картину будут давать разбиения с большим числом кластеров, чем на самом деле. Это произошло из-за отличной от эллиптической формы кластеров, но в целом алгоритм работает неплохо даже в условиях шума.