K-Means: Ллойд, Форджи, Маккуин, Хартиган-Вонг

Ответ 1

R обеспечивает алгоритм Ллойда как вариант для kmeans(); алгоритм по умолчанию, Хартиган и Вонг (1979) гораздо умнее. Как и алгоритм MacQueen (MacQueen, 1967), он обновляет центроиды в любое время, когда перемещается точка; он также делает умные (экономящие время) выбор при проверке ближайшего кластера. С другой стороны, алгоритм Lloyd k-mean является первым и самым простым из всех этих алгоритмов кластеризации.

Алгоритм Ллойда (Ллойд, 1957) принимает ряд наблюдений или случаев (думаю: строки nxp-матрица или точки в Reals) и группирует их в группы k. Он пытается свести к минимуму внутриклассовая сумма квадратов

где u_i - среднее значение всех точек в кластере S_i. Алгоритм выполняется следующим образом (я освободите вас от формальной исчерпывающей нотации):

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

(Здесь я вставлю некоторые изображения, так как мы не можем писать математические формулы с латексом)

Обратите внимание, что решение очень похоже на решение, достигнутое ранее, хотя порядок кластеры произвольны. Что еще более важно, работа заняла только 0.199 секунд! конечно это слишком хорошо, чтобы быть правдой: использование 3 процессорных ядер должно в лучшем случае заняло треть времени нашего первого (последовательного) запуска. Это проблема? Это звучит как бесплатный обед. Здесь нет проблема с бесплатным обедом раз в то время, есть?

Это не всегда работает с R-функциями, но иногда у нас есть возможность смотреть напрямую в коде. Это как-раз тот случай. Я собираюсь поместить этот код в файл, mykmeans.R, и отредактировать его вручную, вставив инструкции cat() в разные места. Здесь умный способ сделать это, используя sink() (хотя это, похоже, не работает в Sweave, оно будет работать в интерактивном режиме):

> sink("mykmeans.R")
> kmeans
> sink()

Теперь отредактируйте файл, изменив имя функции и добавив инструкции cat(). Заметка что вам также нужно удалить завершающую строку::

Затем мы можем повторить наши исследования, но используя mykmeans():

> source("mykmeans.R")
> start.kmeans <- proc.time()[3]
> ans.kmeans <- mykmeans(x, 4, nstart = 3, iter.max = 10, algorithm = "Lloyd")
JJJ statement 1: 0 elapsed time.
JJJ statement 5: 2.424 elapsed time.
JJJ statement 6: 2.425 elapsed time.
JJJ statement 7: 2.52 elapsed time.
JJJ statement 6: 2.52 elapsed time.
JJJ statement 7: 2.563 elapsed time.

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

Вот код:

#######################################################################
# kmeans()

N <- 100000
x <- matrix(0, N, 2)
x[seq(1,N,by=4),] <- rnorm(N/2)
x[seq(2,N,by=4),] <- rnorm(N/2, 3, 1)
x[seq(3,N,by=4),] <- rnorm(N/2, -3, 1)
x[seq(4,N,by=4),1] <- rnorm(N/4, 2, 1)
x[seq(4,N,by=4),2] <- rnorm(N/4, -2.5, 1)
start.kmeans <- proc.time()[3]
ans.kmeans <- kmeans(x, 4, nstart=3, iter.max=10, algorithm="Lloyd")
ans.kmeans$centers
end.kmeans <- proc.time()[3]
end.kmeans - start.kmeans

these <- sample(1:nrow(x), 10000)
plot(x[these,1], x[these,2], pch=".")
points(ans.kmeans$centers, pch=19, cex=2, col=1:4)

library(foreach)
library(doMC)
registerDoMC(3)
start.kmeans <- proc.time()[3]
ans.kmeans.par <- foreach(i=1:3) %dopar% {
  return(kmeans(x, 4, nstart=1, iter.max=10, algorithm="Lloyd"))
}
TSS <- sapply(ans.kmeans.par, function(a) return(sum(a$withinss)))
ans.kmeans.par <- ans.kmeans.par[[which.min(TSS)]]
ans.kmeans.par$centers
end.kmeans <- proc.time()[3]
end.kmeans - start.kmeans

sink("mykmeans.Rfake")
kmeans
sink()

source("mykmeans.R")
start.kmeans <- proc.time()[3]
ans.kmeans <- mykmeans(x, 4, nstart=3, iter.max=10, algorithm="Lloyd")
ans.kmeans$centers
end.kmeans <- proc.time()[3]
end.kmeans - start.kmeans

#######################################################################
# Diving

x <- read.csv("Diving2000.csv", header=TRUE, as.is=TRUE)
library(YaleToolkit)
whatis(x)

x[1:14,c(3,6:9)]

meancol <- function(scores) {
  temp <- matrix(scores, length(scores)/7, ncol=7)
  means <- apply(temp, 1, mean)
  ans <- rep(means,7)
  return(ans)
}
x$panelmean <- meancol(x$JScore)

x[1:14,c(3,6:9,11)]

meancol <- function(scores) {
  browser()
  temp <- matrix(scores, length(scores)/7, ncol=7)
  means <- apply(temp, 1, mean)
  ans <- rep(means,7)
  return(ans)
}

x$panelmean <- meancol(x$JScore)

Вот описание:

Number of cases: 10,787 scores from 1,541 dives (7 judges score each
dive) performed in four events at the 2000 Olympic Games in Sydney,
Australia.

Number of variables: 10.

Description: A full description and analysis is available in an
article in The American Statistician (publication details to be
announced).

Variables:

Event       Four events, men and women 3M and 10m.
Round       Preliminary, semifinal, and final rounds.
Diver       The name of the diver.
Country     The country of the diver.
Rank        The final rank of the diver in the event.
DiveNo      The number of the dive in sequence within round.
Difficulty  The degree of difficulty of the dive.
JScore      The score provided for the judge on this dive.
Judge       The name of the judge.
JCountry    The country of the judge.

И набор данных, чтобы поэкспериментировать с ним https://www.dropbox.com/s/urgzagv0a22114n/Diving2000.csv