planets-hierarchical-clustering.sagews

Иерархическая кластеризация

def scale_data(data, start_row=1, start_col=1) :
    '''
    Стандартизация (масштабирование) данных.
    
    Аргументы:
    data - матрица (список списков) исходных данных
    start_row - строка, с которой начинаются ЧИСЛОВЫЕ данные (по умолчанию: 2-я)
    start_col - столбец, с которого начинаются ЧИСЛОВЫЕ данные (по умолчанию: 2-й)
    
    Результат:
    матрица (список списков) стандартизированных значений
    '''
    num_data = [] # выделяем только числовые значения из data
    for i in range(start_row, len(data)) :
        num_data.append([float(x) for x in data[i][start_col:]])
    num_data = matrix(num_data) # для удобства преобразуем в Sage-матрицу
    # Процедура стандартизации
    means = [mean(x) for x in num_data.columns()]
    stds = [std(list(x)) for x in num_data.columns()]
    for i in range(num_data.nrows()) :
        for j in range(num_data.ncols()) :
            num_data[i,j] = (num_data[i,j]-means[j])/stds[j]
    # Строим матрицу-результат, сохраняя все заголовки
    data_scaled = data[:start_row] # эта часть матрицы не изменилась
    for i in range(start_row, len(data)) :
        new_row = data[i][:start_col] + list(num_data.row(i-start_row))
        data_scaled.append(new_row)
    return data_scaled

# Чтение из файла
f = open('planets.csv', 'r')
planets = [f.readline()[:-1].split(',')]
for line in f :
    planets.append(line[:-1].split(','))
f.close()
sc_planets = scale_data(planets)
for row in sc_planets :
    show(row)

[Name, Distance, Diameter, Period, Mass]

[Mercury, $\displaystyle -0.794766516527$, $\displaystyle -0.754826546633$, $\displaystyle -0.71803204524$, $\displaystyle -0.475055000305$]

[Venus, $\displaystyle -0.771570371374$, $\displaystyle -0.613954355503$, $\displaystyle -0.714130146724$, $\displaystyle -0.467837334642$]

[Earth, $\displaystyle -0.752466789775$, $\displaystyle -0.607002619115$, $\displaystyle -0.710118838069$, $\displaystyle -0.466023371049$]

[Mars, $\displaystyle -0.716369042275$, $\displaystyle -0.721445962213$, $\displaystyle -0.700935521903$, $\displaystyle -0.474520358404$]

[Jupiter, $\displaystyle -0.462570681765$, $\displaystyle 1.8647837206$, $\displaystyle -0.596968538067$, $\displaystyle 2.55031166807$]

[Saturn, $\displaystyle -0.163390173678$, $\displaystyle 1.41950509907$, $\displaystyle -0.413220157144$, $\displaystyle 0.433320683788$]

[Uranus, $\displaystyle 0.502115623854$, $\displaystyle 0.130218306337$, $\displaystyle 0.998384689104$, $\displaystyle -0.311359106842$]

[Neptune, $\displaystyle 1.25458654003$, $\displaystyle 0.0873135368213$, $\displaystyle 0.991962789687$, $\displaystyle -0.313268542202$]

[Pluto, $\displaystyle 1.90443141151$, $\displaystyle -0.804591179369$, $\displaystyle 1.86305776836$, $\displaystyle -0.475568638417$]

Вычисляем матрицу расстояний

Не вычисляйте матрицу расстояний для нестандартизированных данных!

def dist(data, metric=2, start_row=1, start_col=1) :
    '''
    Вычисление матрицы расстояний между объектами.
    
    Аргументы:
    data - матрица СТАНДАРТИЗИРОВАННЫХ исходных данных; объекты - строки; признаки - столбцы
    metric - вид метрики, целое число (1 - манхэттенская, 2 - евклидова, любое целое число p - p-метрика Минковского)
    start_row - строка, с которой начинаются ЧИСЛОВЫЕ данные
    start_col - столбец, с которого начинаются ЧИСЛОВЫЕ данные
    
    Результат:
    матрица расстояний между объектами в указанной метрике
    '''
    p = metric # для удобства
    num_data = matrix([[float(x) for x in row[start_col:]] for row in data[start_row:]])
    n_obj = num_data.nrows() # число объектов
    d = lambda x, y : numerical_approx((sum([abs(x[i]-y[i])^p for i in range(len(x))]))^(1/p)) # мини-функция для расстояния
    dist_matrix = zero_matrix(RR, n_obj)
    for i in range(n_obj) :
        for j in range(n_obj) :
            if i != j :
                dist_matrix[i, j] = d(num_data[i], num_data[j])
    return dist_matrix

# Проверяем результат
show(dist(sc_planets))

$\displaystyle \left(\begin{array}{rrrrrrrrr} 0.000000000000000 & 0.143004737311175 & 0.154225061533540 & 0.0869080222846028 & 4.01748825458545 & 2.45853835894471 & 2.33197136116222 & 2.80344521331895 & 3.73499275274658 \\ 0.143004737311175 & 0.000000000000000 & 0.0208003492699203 & 0.121739037199194 & 3.91951189604482 & 2.32539740738687 & 2.26567020189974 & 2.74440031804796 & 3.72011898510131 \\ 0.154225061533540 & 0.0208003492699203 & 0.000000000000000 & 0.120651805805302 & 3.91214104056470 & 2.31316089955092 & 2.24952818664257 & 2.72594091982160 & 3.70398140534323 \\ 0.0869080222846028 & 0.121739037199194 & 0.120651805805302 & 0.000000000000000 & 3.98916225568661 & 2.40757399988922 & 2.26370220688759 & 2.72592455440577 & 3.66736565358439 \\ 4.01748825458545 & 3.91951189604482 & 3.91214104056470 & 3.98916225568661 & 0.000000000000000 & 2.19161956603874 & 3.83061984384455 & 4.10278241728561 & 5.28544627478770 \\ 2.45853835894471 & 2.32539740738687 & 2.31316089955092 & 2.40757399988922 & 2.19161956603874 & 0.000000000000000 & 2.15692712741840 & 2.51342931183290 & 3.90256471542637 \\ 2.33197136116222 & 2.26567020189974 & 2.24952818664257 & 2.26370220688759 & 3.83061984384455 & 2.15692712741840 & 0.000000000000000 & 0.753722883872243 & 1.90131077079040 \\ 2.80344521331895 & 2.74440031804796 & 2.72594091982160 & 2.72592455440577 & 4.10278241728561 & 2.51342931183290 & 0.753722883872243 & 0.000000000000000 & 1.41525268515808 \\ 3.73499275274658 & 3.72011898510131 & 3.70398140534323 & 3.66736565358439 & 5.28544627478770 & 3.90256471542637 & 1.90131077079040 & 1.41525268515808 & 0.000000000000000 \end{array}\right)$

Запускаем агломеративный метод

def mmin(M) :
    '''
    Поиск максимального элемента и его позиции в симметричной квадратной матрице М.
    '''
    imin = 0
    jmin = 1
    for i in range(M.nrows()) :
        for j in range(M.ncols()) :
            if i < j and M[i,j] < M[imin,jmin] :
                imin = i
                jmin = j
    return M[imin,jmin], imin, jmin

def mreduce(M, i0) :
    '''
    Удаление строки и столбца с номером i0 из матрицы М.
    
    Аргументы:
    M - квадратная Sage-матрица
    i0 - целое число в пределах допустимых индексов
    
    Результат: новая матрица
    '''
    n = M.nrows()
    return M.matrix_from_rows_and_columns([i for i in range(n) if i != i0], [j for j in range(n) if j != i0])

c = [[row[0]] for row in planets[1:]] # кластеры
planets_dist = dist(sc_planets)
n = len(c) # число кластеров
while n > 1 : # пока число кластеров больше одного
    dmin, imin, jmin = mmin(planets_dist) # находим два кластера для объединения
    print 'Объединяем', c[imin], 'и', c[jmin], 'на расстоянии', dmin
    new = min(imin, jmin) # номер нового кластера
    old = max(imin, jmin) # этот номер удаляется
    c[new].extend(c[old])
    del c[old]
    for j in range(n) : # пересчёт матрицы расстояний
        planets_dist[new, j] = min(planets_dist[imin, j], planets_dist[jmin, j]) # метод ближайшего соседа
    planets_dist = mreduce(planets_dist, old)
    n = len(c)

Объединяем ['Venus'] и ['Earth'] на расстоянии 0.0208003492699203 Объединяем ['Mercury'] и ['Mars'] на расстоянии 0.0869080222846028 Объединяем ['Mercury', 'Mars'] и ['Venus', 'Earth'] на расстоянии 0.121739037199194 Объединяем ['Uranus'] и ['Neptune'] на расстоянии 0.753722883872243 Объединяем ['Uranus', 'Neptune'] и ['Pluto'] на расстоянии 1.41525268515808 Объединяем ['Saturn'] и ['Uranus', 'Neptune', 'Pluto'] на расстоянии 2.15692712741840 Объединяем ['Jupiter'] и ['Saturn', 'Uranus', 'Neptune', 'Pluto'] на расстоянии 2.19161956603874 Объединяем ['Mercury', 'Mars', 'Venus', 'Earth'] и ['Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto'] на расстоянии 3.91214104056470

Используем R

hclust()

%r
# Без масштабирования
library(cluster)
planets = read.table('planets.csv', header=TRUE, sep=",", row.names=1) # чтение из файла
planets.dist <- daisy(planets) # вычисление матрицы расстояний
planets.h <- hclust(planets.dist, method="ward.D") # кластеризация
plot(planets.h, main="Solar System Planets + Pluto") # рисуем картинку

%r
# С масштабированием
planets = read.table('planets.csv', header=TRUE, sep=",", row.names=1)
planets.dist <- daisy(scale(planets))
planets.h <- hclust(planets.dist, method="ward.D")
plot(planets.h, main="Solar System Planets + Pluto")

agnes()

%r
planets = read.table('planets.csv', header=TRUE, sep=",", row.names=1)
planets.a <- agnes(planets, diss=FALSE, metric="ward", stand=TRUE, trace.lev=2)
plot(planets.a, main="Solar System Planets + Pluto")

C agnes(n=9, method = 1, ..): 8 merging steps nmerge=0, j=2, d_min = D(2,3) = 0.0257017; last=3; size(A_new)= 2 nmerge=1, j=2, d_min = D(1,4) = 0.107353; last=4; upd(n,b); size(A_new)= 2 nmerge=2, j=2, d_min = D(1,2) = 0.172696; last=3; size(A_new)= 4 nmerge=3, j=5, d_min = D(7,8) = 0.926559; last=8; size(A_new)= 2 nmerge=4, j=5, d_min = D(7,9) = 2.04475; last=9; size(A_new)= 3 nmerge=5, j=5, d_min = D(1,6) = 3.12381; last=6; upd(n,b); size(A_new)= 5 nmerge=6, j=5, d_min = D(1,7) = 3.52082; last=7; upd(n,b); size(A_new)= 8 nmerge=7, j=5, d_min = D(1,5) = 5.46740; last=9; size(A_new)= 9

diana()

%r
planets = read.table('planets.csv', header=TRUE, sep=",", row.names=1)
planets.d <- diana(planets, diss=FALSE, metric="ward", stand=TRUE, trace.lev=2)
plot(planets.d, main="Solar System Planets + Pluto")

C diana(): ndist= 37, diameter = 7.02223

Проверка качества кластеризации

pvclust()

%r
library(pvclust)
planets.pv <- pvclust(scale(t(planets)), method.dist="euclid", method.hclust="ward.D", nboot=1000) # транспонирование!!!
plot(planets.pv, main="Planets")

Bootstrap (r = 0.5)... Done. Bootstrap (r = 0.5)... Done. Bootstrap (r = 0.5)... Done. Bootstrap (r = 0.75)... Done. Bootstrap (r = 0.75)... Done. Bootstrap (r = 1.0)... Done. Bootstrap (r = 1.0)... Done. Bootstrap (r = 1.0)... Done. Bootstrap (r = 1.25)... Done. Bootstrap (r = 1.25)... Done.

Warning message in a$p[] <- c(1, bp[r == 1]): “number of items to replace is not a multiple of replacement length”