MDS.sagews

Многомерное шкалирование

Решение примера 9.2 из учебно-методического пособия (стр. 116).

1. Подготовка данных к работе

# Создаём матрицу исходных данных (без заголовков)
data = [['Firm 1', 10, 320, 50], ['Firm 2', 1.5, 80, 10], ['Firm 3', 2.0, 120, 10], ['Firm 4', 2.5, -20, 30]]
# Простой вывод на экран
for row in data :
    print row

['Firm 1', 10, 320, 50] ['Firm 2', 1.50000000000000, 80, 10] ['Firm 3', 2.00000000000000, 120, 10] ['Firm 4', 2.50000000000000, -20, 30]

2. Стандартизация данных

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

# Стандартизируем наши данные
sc_data = scale_data(data, 0, 1)
show(sc_data)

[[Firm 1, $\displaystyle 1.49224801012$, $\displaystyle 1.36638986991$, $\displaystyle 1.30558241967$], [Firm 2, $\displaystyle -0.621770004217$, $\displaystyle -0.31532073921$, $\displaystyle -0.783349451801$], [Firm 3, $\displaystyle -0.497416003374$, $\displaystyle -0.03503563769$, $\displaystyle -0.783349451801$], [Firm 4, $\displaystyle -0.37306200253$, $\displaystyle -1.01603349301$, $\displaystyle 0.261116483934$]]

# Выделим только числовой массив данных
num_data = [row[1:] for row in sc_data]
num_data = matrix(RDF, num_data)
show(r'$Z='+latex(num_data)+'$')

$Z= \left(\begin{array}{rrr} 1.49224801012 & 1.36638986991 & 1.30558241967 \\ -0.621770004217 & -0.31532073921 & -0.783349451801 \\ -0.497416003374 & -0.03503563769 & -0.783349451801 \\ -0.37306200253 & -1.01603349301 & 0.261116483934 \end{array}\right)$

3. Построение матрицы $\gamma$ расстояний между объектами

n_obj = num_data.nrows() # количество объектов
gamma = matrix(RDF, n_obj) # пустая матрица
for i in range(n_obj) :
    for j in range(n_obj) :
        x = num_data.row(i)
        y = num_data.row(j)
        gamma[i,j] = norm(x-y) # воспольуемся тем, что евклидово расстояние - это норма разности двух векторов
show(r'$\gamma='+latex(gamma)+'$')

$\gamma= \left(\begin{array}{rrrr} 0.0 & 3.41480000899 & 3.20724066817 & 3.20097354166 \\ 3.41480000899 & 0.0 & 0.306632770036 & 1.28209325884 \\ 3.20724066817 & 0.306632770036 & 0.0 & 1.43830796444 \\ 3.20097354166 & 1.28209325884 & 1.43830796444 & 0.0 \end{array}\right)$

4. Построение матрицы $B$

$$b_{ij}=\frac{1}{2}\left(-\gamma_{ij}^2+\frac{1}{n}\sum\limits_{i=1}^{n}\gamma_{ij}^2+\frac{1}{n}\sum\limits_{j=1}^{n}\gamma_{ij}^2-\frac{1}{n^2}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\gamma_{ij}^2\right)$$

B = matrix(RDF, n_obj) # пустая матрица
ss = sum([sum([x^2 for x in gamma.row(i)]) for i in range(n_obj)]) # общая сумма квадратов
for i in range(n_obj) :
    si = sum([x^2 for x in gamma.row(i)]) # сумма квадратов по i-й строке
    for j in range(n_obj) :
        sj = sum([x^2 for x in gamma.column(j)]) # сумма квадратов по j-му столбцу
        B[i,j] = 0.5*(-gamma[i,j]^2+(si+sj)/n_obj-ss/(n_obj^2))
show(r'$B='+latex(B)+'$')

$B= \left(\begin{array}{rrrr} 5.79837085485 & -2.3814133881 & -1.81286765439 & -1.60408981236 \\ -2.3814133881 & 1.09966147036 & 0.933962177327 & 0.34778974042 \\ -1.81286765439 & 0.933962177327 & 0.862286539957 & 0.0166189371059 \\ -1.60408981236 & 0.34778974042 & 0.0166189371059 & 1.23968113483 \end{array}\right)$

5. Нахождение собственных значений и векторов матрицы $B$

L = B.eigenvectors_right() # получаем СЗ и нормированные СВ
show(L)
# values = [x[0] for x in L]  # отдельно СЗ
# vectors = [x[1] for x in L]  # отдельно СВ

[($\displaystyle 7.7908609798$, [$\displaystyle \left(0.862446891087,\,-0.357241703659,\,-0.274369363339,\,-0.230835824089\right)$], $\displaystyle 1$), ($\displaystyle 1.20899831927$, [$\displaystyle \left(0.0531391112958,\,0.295127073868,\,0.477569131312,\,-0.825835316476\right)$], $\displaystyle 1$), ($\displaystyle -8.39453994591 \times 10^{-16}$, [$\displaystyle \left(-0.5,\,-0.500000000002,\,-0.499999999998,\,-0.5\right)$], $\displaystyle 1$), ($\displaystyle 0.00014070092767$, [$\displaystyle \left(-0.0579792627207,\,-0.7316272107,\,0.668318170693,\,0.121288302727\right)$], $\displaystyle 1$)]

L = sorted(L, key=lambda x: x[0], reverse=True) # сортируем в порядке убывания СЗ
show(L)

[($\displaystyle 7.7908609798$, [$\displaystyle \left(0.862446891087,\,-0.357241703659,\,-0.274369363339,\,-0.230835824089\right)$], $\displaystyle 1$), ($\displaystyle 1.20899831927$, [$\displaystyle \left(0.0531391112958,\,0.295127073868,\,0.477569131312,\,-0.825835316476\right)$], $\displaystyle 1$), ($\displaystyle 0.00014070092767$, [$\displaystyle \left(-0.0579792627207,\,-0.7316272107,\,0.668318170693,\,0.121288302727\right)$], $\displaystyle 1$), ($\displaystyle -8.39453994591 \times 10^{-16}$, [$\displaystyle \left(-0.5,\,-0.500000000002,\,-0.499999999998,\,-0.5\right)$], $\displaystyle 1$)]

6. Вычисление матрицы $X$

$$X=V{\Lambda}^{\frac{1}{2}}$$

# Матрица V
V = []
for i in range(2) : # берём только две шкалы
    V.append(L[i][1][0])
V = matrix(V).transpose()
show(V)

$\displaystyle \left(\begin{array}{rr} 0.862446891087 & 0.0531391112958 \\ -0.357241703659 & 0.295127073868 \\ -0.274369363339 & 0.477569131312 \\ -0.230835824089 & -0.825835316476 \end{array}\right)$

# Матрица sqrt(Lambda)
Lam = matrix(RDF, 2)
for i in range(2) :
    Lam[i,i] = sqrt(L[i][0])
show(Lam)

$\displaystyle \left(\begin{array}{rr} 2.79121138214 & 0.0 \\ 0.0 & 1.09954459631 \end{array}\right)$

# Матрица X
X = V*Lam
show(r'$X='+latex(X)+'$')

$X= \left(\begin{array}{rr} 2.40727157889 & 0.058428822678 \\ -0.997137109429 & 0.324505379296 \\ -0.765822889862 & 0.525108557698 \\ -0.644311579604 & -0.908042759671 \end{array}\right)$

7. Рисунок

pts = X.rows() # преобразуем Х в список точек (строк)
g = Graphics()
for i in range(len(pts)) :
    g += point(pts[i], size=200, rgbcolor=(random(), random(), random())) + text('Firm '+str(i+1), (pts[i][0], pts[i][1]+0.1))
g.show(axes_labels=['Scale 1', 'Scale 2'], xmin=-1.5, xmax=2.7, ymin=-1, ymax=0.7)