Кластеризация последовательностей
By Firepush
Как-то раз меня спросили на работе: “А есть ли у нас некие шаблоны поведения, которым следуют игроки при прохождении определенного этапа игры? И если есть – сколько таких шаблонов и какие они?”
Вопрос интересный и на первый взгляд простой. Но, как это обычно бывает, пришлось повозиться.
Первое решение, которое приходит в голову – кластеризация на основе алгоритмов расстояния между словами (расстояние Левенштейна, Хэмминга и т.д.). Но самое неприятное – это то, что готовые библиотеки применяют алгоритмы только к “словам”. Но не к последовательности элементов в массиве. (Очень странно, что очевидное решение этой проблемы не пришло мне в голову сразу, но об этом – в конце статьи). Расстояние Хэмминга – очень прост в реализации, но слишком чувствителен. Левенштайн мне переписать было не под силу.
Так что я решил попробовать написать свой велосипед)
Итак, определимся, чем для нас является шаблон поведения в игре:
Шаблон поведения – это последовательность из уникальных событий.
В контексте игры последовательность может выглядеть так:
установка – завершение таска №1 квеста №1 – завершение №2 квеста №1 – завершение №1 – получение левела №1 – трата рубинов №1 и т.д.
Почему уникальных? Для простоты алгоритма, признаюсь. Но для большинства событий это условие и так выполнятся: “квест №1”, например, можно завершить только один раз. А для таких, как “трата рубинов”, к названию будем приписывать порядковый номер события данного типа: трата рубинов №1, трата рубинов №2 и т.д.
Подготовка данных для анализа потребовала некоторых усилий. В конечном итоге в R мы получаем Named List, где имя элемента листа – это ID соответствующего игрока, а значение – вектор строк с уникальными элементами последовательности.
На время забудем об игре. Представим, что у нас есть две последовательности из всех букв английского алфавита:
a <- LETTERS[1:6]
b <- LETTERS[c(1,3,2,4,5,6)]
> a
[1] “A” “B” “C” “D” “E” “F”
> b
[1] “A” “C” “B” “D” “E” “F”
1
2
3
4
5
6
a <- LETTERS[1:6]
b <- LETTERS[c(1,3,2,4,5,6)]
> a
[1] “A” “B” “C” “D” “E” “F”
> b
[1] “A” “C” “B” “D” “E” “F”
Как можно охарактеризовать “расстояние” между двумя этими последовательностями? Первое, что приходит на ум – сравнить индексы позиций соответствующих элементов и посчитать среднюю “разницу”. Чем она больше, тем “дальше” они друг от друга.
Функция для подсчета данного коэффициента в R:
rangdiff <- function(x,y){
xdf <- as.data.frame(cbind(x, 1:length(x)), stringsAsFactors = F)
colnames(xdf) <- c(“elem”, “rangx”)
ydf <- as.data.frame(cbind(y, 1:length(y)), stringsAsFactors = F)
colnames(ydf) <- c(“elem”, “rangy”)
xydf <- inner_join(xdf, ydf, “elem”)
xydf <- mutate(xydf, rangx = as.numeric(rangx), rangy = as.numeric(rangy))
xydf <- mutate(xydf, diff = abs(rangy-rangx)/(max(length(x),length(y))/2))
diffres <- round(mean(xydf$diff),4)
diffres
}
1
2
3
4
5
6
7
8
9
10
11
12
13
rangdiff <- function(x,y){
xdf <- as.data.frame(cbind(x, 1:length(x)), stringsAsFactors = F)
colnames(xdf) <- c(“elem”, “rangx”)
ydf <- as.data.frame(cbind(y, 1:length(y)), stringsAsFactors = F)
colnames(ydf) <- c(“elem”, “rangy”)
xydf <- inner_join(xdf, ydf, “elem”)
xydf <- mutate(xydf, rangx = as.numeric(rangx), rangy = as.numeric(rangy))
xydf <- mutate(xydf, diff = abs(rangy-rangx)/(max(length(x),length(y))/2))
diffres <- round(mean(xydf$diff),4)
diffres
}
Мы получаем таблицу такого вида:
elem rangx rangy diff
1 A 1 1 0.0000000
2 C 3 2 0.3333333
3 B 2 3 0.3333333
4 D 4 4 0.0000000
5 E 5 5 0.0000000
6 F 6 6 0.0000000
inner_join – позволяет нам игнорировать элементы, которые не встретились в обеих последовательностях. rangx и rangy – индексы соответствующих элементов в первом и втором векторе. diff считается для каждой пары этих индексов, как абсолютная разность индексов, деленная на половину длины наибольшей последовательности. В конечном итоге получаем diffres, как среднее значение всех diff-ов, в нашем случае – 0.1111.
Если сделать две перестановки, то расстояние уже будет равно 0.2222:
> a
[1] “A” “B” “C” “D” “E” “F”
> c
[1] “A” “C” “B” “D” “F” “E”
1
2
3
4
> a
[1] “A” “B” “C” “D” “E” “F”
> c
[1] “A” “C” “B” “D” “F” “E”
Ну а если сравнить a с “совсем” противоположным d, то получим 1:
> d
[1] “F” “E” “D” “C” “B” “A”
1
2
> d
[1] “F” “E” “D” “C” “B” “A”
Таким образом, самые близкие последовательности будут иметь расстояние близкое к нулю, а самые далекие – к единице.
Как ни странно, даже с этим простым алгоритмом уже можно делать кластеризацию!
Но давайте поразмыслим еще. Мы ведь совсем не учитываем “липкость” элементов. Давайте сравним три такие последовательности:
> a
[1] “A” “B” “C” “D” “E” “F”
> b
[1] “B” “A” “D” “C” “F” “E”
> c
[1] “D” “E” “F” “A” “B” “C”
1
2
3
4
5
6
> a
[1] “A” “B” “C” “D” “E” “F”
> b
[1] “B” “A” “D” “C” “F” “E”
> c
[1] “D” “E” “F” “A” “B” “C”
В b – элементы попарно переставлены. А в с – первую тройку элементов поменяли местами с последней тройкой. Теперь сравним расстояния:
> rangdiff(a, b)
[1] 0.3333
> rangdiff(a, c)
[1] 1
1
2
3
4
> rangdiff(a, b)
[1] 0.3333
> rangdiff(a, c)
[1] 1
Но резонно ли говорить, что c гораздо дальше от a, чем b? Ведь в с есть такие же подпоследовательности, как и в a, тогда как в b – ни одной?
Поэтому нам нужно придумать коэффициент, учитывающий наличие совпадающих “лоскутков”, как я их для себя назвал.
Упрощенно он рассчитывается так: из двух рассматриваемых последовательностей берется та, что с меньшей длиной. Смотрим первый элемент этой последовательности, смотрим какой элемент за ним следует и проверяем повторяется ли эта подпоследовательность у соседа. Если да, идем дальше – попутно запоминая длину рассматриваемого “лоскутка”. Если лоскуток “обрывается” – запоминаем длину этого лоскутка и начинаем считать длину для нового. При этом в алгоритме есть понятие нулевой длины, когда лоскуток – из одного элемента.
Поделив среднюю длину “лоскутков”, на длину более длинной последовательности и применив некоторые магические арифметические преобразования, получим другой коэффициент расстояния.
> meanlen(a,b)
[1] “B” “A” “D” “C” “F” “E”
[1] “A” “B” “C” “D” “E” “F”
[1] 1
> meanlen(a,c)
[1] “D” “E” “F” “A” “B” “C”
[1] “A” “B” “C” “D” “E” “F”
[1] 0.1308987
1
2
3
4
5
6
7
8
> meanlen(a,b)
[1] “B” “A” “D” “C” “F” “E”
[1] “A” “B” “C” “D” “E” “F”
[1] 1
> meanlen(a,c)
[1] “D” “E” “F” “A” “B” “C”
[1] “A” “B” “C” “D” “E” “F”
[1] 0.1308987
Как и ожидалось, этот коэффициент считает по-другому. И по-своему правильно 🙂
Сама функция выглядит так (простите за отстойный код):
meanlen <- function(x,y) {
if(length(y) > length(x)){
mn <- x
mx <- y
}else {
mn <- y
mx <- x
}
seqLenVec <- c();
seqLength = 1
for(q in 1:(length(mn)-1)) {
xindex <- q
xletter <- mn[xindex];
yindex <- which(mx == xletter)
nexXIndex <- xindex + 1
nextXLetter <- mn[nexXIndex]
nextYindex <- which(mx == nextXLetter)
if(length(yindex) != 0 & length(nextYindex) != 0) {
if(nextYindex – yindex == 1){
seqLength <- seqLength + 1
}else{
seqLenVec <- c(seqLenVec, seqLength)
seqLength <- 1
}
}else{
seqLenVec <- c(seqLenVec, seqLength)
seqLength <- 1
}
if(q == length(mn)-1){
seqLenVec <- c(seqLenVec, seqLength)
}
}
seqLenVec <- seqLenVec – 1
# -7 – это магия логарифма
uper <- ifelse(mean(seqLenVec) == 0, -7, log(mean(seqLenVec)/(length(mn)-1)))
#нормализация “на глазок”. тоже магия
sim <- uper/7 + 1
1 – sim
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
meanlen <- function(x,y) {
if(length(y) > length(x)){
mn <- x
mx <- y
}else {
mn <- y
mx <- x
}
seqLenVec <- c();
seqLength = 1
for(q in 1:(length(mn)-1)) {
xindex <- q
xletter <- mn[xindex];
yindex <- which(mx == xletter)
nexXIndex <- xindex + 1
nextXLetter <- mn[nexXIndex]
nextYindex <- which(mx == nextXLetter)
if(length(yindex) != 0 & length(nextYindex) != 0) {
if(nextYindex – yindex == 1){
seqLength <- seqLength + 1
}else{
seqLenVec <- c(seqLenVec, seqLength)
seqLength <- 1
}
}else{
seqLenVec <- c(seqLenVec, seqLength)
seqLength <- 1
}
if(q == length(mn)-1){
seqLenVec <- c(seqLenVec, seqLength)
}
}
seqLenVec <- seqLenVec – 1
# -7 – это магия логарифма
uper <- ifelse(mean(seqLenVec) == 0, -7, log(mean(seqLenVec)/(length(mn)-1)))
#нормализация “на глазок”. тоже магия
sim <- uper/7 + 1
1 – sim
}
Добавляем шумы и строим графики
Давайте проверим на графиках, как работают обе функции определения расстояния в зависимости от силы шума, добавляемой к одной и той же последовательности. Пусть дана последовательность a <- LETTERS (26 букв алфавита).
Создадим функцию shuffle, которая принимает на вход два аргумента: кол-во перестановок и саму последовательность:
shuffle <- function(x, seq){
temp <- seq
indexesToShuffle <- sample(26, x)
temp[indexesToShuffle] <- temp[sample(indexesToShuffle)]
temp
}
1
2
3
4
5
6
shuffle <- function(x, seq){
temp <- seq
indexesToShuffle <- sample(26, x)
temp[indexesToShuffle] <- temp[sample(indexesToShuffle)]
temp
}
Чем больше x, тем сильнее последовательность на выходе будет не похожа на исходную.
Создадим вектор “иксов” с нарастанием, которые будем скармливать этой функции:
numb <- unlist(lapply(2:26, rep, 10))
1
numb <- unlist(lapply(2:26, rep, 10))
Теперь нагенеририруем последовательности с нарастающей шумностью и сохраним в листе gro:
gro <- lapply(numb, shuffle, a)
1
gro <- lapply(numb, shuffle, a)
Ну и наконец посчитаем двумя способами расстояния между a и каждой из последовательности листа gro
inde1 <- sapply(gro, rangdiff, x)
inde2 <- sapply(gro, meanlen, x)
1
2
3
inde1 <- sapply(gro, rangdiff, x)
inde2 <- sapply(gro, meanlen, x)
plot(inde1, main=”Расстояние ~ шумность (rangdiff)”, xlab=”шум”, ylab=”расстояние”)
1
plot(inde1, main=”Расстояние ~ шумность (rangdiff)”, xlab=”шум”, ylab=”расстояние”)
iumG2NA
plot(inde2, main=”Расстояние ~ шумность (meanlen)”, xlab=”шум”, ylab=”расстояние”)
1
plot(inde2, main=”Расстояние ~ шумность (meanlen)”, xlab=”шум”, ylab=”расстояние”)
8pEqUWy
Что же, вполне неплохо. Чтобы учитывать и тот и другой фактор, будем считать среднее геометрическое двух коэффициентов:
geomMean <- function(x,y){ sqrt(rangdiff(x,y) * meanlen(x,y)) }
inde3 <- sapply(gro, geomMean , x)
1
2
geomMean <- function(x,y){ sqrt(rangdiff(x,y) * meanlen(x,y)) }
inde3 <- sapply(gro, geomMean , x)
plot(inde3, main=”Расстояние ~ шумность (geomMean )”, xlab=”шум”, ylab=”расстояние”)
1
plot(inde3, main=”Расстояние ~ шумность (geomMean )”, xlab=”шум”, ylab=”расстояние”)
geomMean
Приступаем к кластеризации
В R есть функция, позволяющая построить матрицу соответствий для векторов. Но, чтобы применить ее для листов, нужно немного переработать:
f1 <- function(a,b, fun) {
outer(a, b, function(x,y) vapply(seq_along(x), function(i) fun(x[[i]], y[[i]]), numeric(1)))
}
1
2
3
f1 <- function(a,b, fun) {
outer(a, b, function(x,y) vapply(seq_along(x), function(i) fun(x[[i]], y[[i]]), numeric(1)))
}
Теперь создадим четыре исходные тестовые последовательности, вокруг которых будем плясать.
a <- LETTERS
b <- LETTERS[26:1]
c <- LETTERS[sample(26)]
d <- LETTERS[sample(26)]
1
2
3
4
a <- LETTERS
b <- LETTERS[26:1]
c <- LETTERS[sample(26)]
d <- LETTERS[sample(26)]
Генерация шума:
numb <- unlist(lapply(2:7, rep, 10))
gro <- lapply(numb, shuffle, a)
gro2 <- lapply(numb, shuffle, b)
gro3 <- lapply(numb, shuffle, c)
gro4 <- lapply(numb, shuffle, d)
1
2
3
4
5
numb <- unlist(lapply(2:7, rep, 10))
gro <- lapply(numb, shuffle, a)
gro2 <- lapply(numb, shuffle, b)
gro3 <- lapply(numb, shuffle, c)
gro4 <- lapply(numb, shuffle, d)
Соберем в кучу и перемешаем:
bigGro <- c(gro, gro2, gro3, gro4)
bigGro <- BigGro[sample(length(bigGro ))]
1
2
bigGro <- c(gro, gro2, gro3, gro4)
bigGro <- BigGro[sample(length(bigGro ))]
А теперь посчитаем матрицу расстояний. Внимание! Эта операция ресурсоемкая. Поэтому для больших объемов данных лучше делать выборки.
ptm <- proc.time()
res <- f1(bigGro,bigGro, meanGeom)
proc.time() – ptm
1
2
3
ptm <- proc.time()
res <- f1(bigGro,bigGro, meanGeom)
proc.time() – ptm
Примерно за 80 секунд получаем такую матрицу (первые 5 строк и столбцов):
KaJkikZ
Каждый элемент означает дистанцию между последовательностями с индексом соответствующей строки и столбца. На основе этой матрицы расстояний и строится кластерный анализ. Я использовал R-библиотеку “cluster”.
Итак:
res <- f1(bigGro,bigGro, meanGeom)
pamx <- pam(res, 4)
clusplot(pamx)
1
2
3
res <- f1(bigGro,bigGro, meanGeom)
pamx <- pam(res, 4)
clusplot(pamx)
LWR0teq
res2 <- f1(bigGro,bigGro, rangdiff)
pamx <- pam(res2, 4)
clusplot(pamx)
1
2
3
res2 <- f1(bigGro,bigGro, rangdiff)
pamx <- pam(res2, 4)
clusplot(pamx)
rangdiff
res3 <- f1(bigGro,bigGro, meanlen)
pamx <- pam(res3, 4)
clusplot(pamx)
1
2
3
res3 <- f1(bigGro,bigGro, meanlen)
pamx <- pam(res3, 4)
clusplot(pamx)
meanlen
Вывод: применение двух коэффициентов дает лучшее разбиение на кластеры, чем по отдельности, так как близкие группы не накладываются друг на друга, как это видно выше.
Мысль о том, что каждый элемент можно сопоставить отдельному символу (что по сути и было сделано) и “склеить” в слово до меня дошла уже после проделанной работы. Применив код ниже (который отработал в разы быстрее):
library(stringdist)
bigGro2 <- lapply(bigGro, paste, collapse = “”)
res4 <- f1(bigGro2,bigGro2, stringdist)
pamx <- pam(res4, 4)
clusplot(pamx, main=”Кластеризация (Левенштайн)”)
1
2
3
4
5
library(stringdist)
bigGro2 <- lapply(bigGro, paste, collapse = “”)
res4 <- f1(bigGro2,bigGro2, stringdist)
pamx <- pam(res4, 4)
clusplot(pamx, main=”Кластеризация (Левенштайн)”)
Получим вот такую замечательную картинку:
lev
Очевидно, что применение этого алгоритма намного эффективнее. Как в плане точности, так и в скорости. При этом элементы последовательности могут повторяться. Недостаток этого подхода – кол-во символов ограничено, тогда как уникальных элементов в наборе может быть очень много. Я пытался использовать список разных символов из юникода, но R у меня отказывается их обрабатывать.
Так а что там с игрой? Кластеризация не позволила выделить какие-то обособленные группы:
seq
Что же, плохой опыт тоже опыт)