自己組織化マップ（SOM: Self-Organizing Map）試作記
自己組織化マップ（SOM）の簡明な解説（神野ら、2009）があったので、Pythonで実装してみた。以下は忘備録である。
なお、プログラムファイルなどは圧縮ファイルfiles_som.zipにまとめた。ダウンロード解凍すれば自由に利用できる。
SOMは、高次元空間のデータ
を低次元（２次元など）上のニューロンと呼ばれる点に対応させる方法である。ニューロンは、２次元の場合は平面上に配置される。各ニューロンは参照ベクトルと呼ばれるベクトル
をもち、参照ベクトルは入力データと同次元である。入力データがいずれかの参照ベクトルとよく対応するようにの値が調整される。具体的なアルゴリズムは以下のようである（神野ら、2009）。
SOMのアルゴリズム
入力データの次元をn、その個数をJMとおく。入力データの成分の値は、区間[0,1]の範囲になるようにスケーリングしておく。２次元の格子
上にニューロンを配置し、各ニューロンにn次元参照ベクトル
を付ける。
ステップ０：
参照ベクトル
を一様乱数で初期化する。
ステップ数tの値を０とおく。
更新の繰り返し最大数をt_maxとおく。
ステップ１：
各入力データ
に対して以下の参照ベクトルの更新手続きを行う。
に最も近い参照ベクトルを
とおく。
すべての参照ベクトルを次式で更新する。
ここで、
である。
すべての入力データによる参照ベクトルの更新が終わったら、
とおく。
ステップ２：
ならば、ステップ１に戻る。
ステップ３：
各入力データを、それぞれの最も近い参照ベクトルの位置
にプロットする。
上のアルゴリズムをリスト１のようにPythonで実装した。
リスト１のプログラムを実行すると、入力データファイル名の入力が求められる（図）。図１では、ExcelデータファイルAnimals.xlsxが入力されている。
図１
データファイルは、図２の形式で用意する。図２は古川（2017）の表1.1をExcelで用意したものである。
図２ ファイル名Animals.xlsx
第１行目に変数名、第1列にケース名が設定されている。データ値は２行目以降に入力されている。データ値の１は「該当する」、０は「該当しない」を表す。
図１において、データファイル名を設定するとファイルからデータが読み込まれて計算が始まるが、モジュールxlrdがないなどのエラーが出たときは、無いというエラーの出たモジュールをpipコマンドなどでインストールしてから実行し直せばよい。
計算が終了すると、図３のようなマップが表示される。
図３
各入力データのケース名が、一番近い参照ベクトルのニューロンの格子上の位置に表示されている。文字列の左下の角がその位置である。色は、参照ベクトルの近傍の傾きを表し、最も緩やかなところが緑、最も急峻なところが赤であるが、アルファ値で透明度の調節を行っている。
図３では、上に鳥類、下に哺乳類がまとまっているが、SOMにおける計算において初期値を乱数でランダムに与えているので、鳥類と哺乳類のグループ分けが左右に分かれることもある。
図２のデータの場合、図３に示されるように鳥類と哺乳類の２つのグループに分かれたマップが得られた。次に、ニューロンを２個とした場合のSOMを求めてみる。ニューロンが２個の場合のSOMを求めるプログラムをリスト３に示す。ニューロンの個数が２個であることを除いてリスト１とアルゴリズムは同じである。
リスト３のプログラムを実行すると、入力データファイル名の設定が求められる（図４）。
図４
図２のファイル名を設定してEnterキーを押すと計算が始まる。
計算が終了すると図５のグラフが表示される。
図５
横軸が各ケースとニューロン０の参照ベクトルとの距離の２乗であり、縦軸がニューロン１の参照ベクトルとの距離の２乗を表す。哺乳類がニューロン０の参照ベクトルに近く、鳥類がニューロン１の参照ベクトルに近いことが分かる。
古川徹生（2017）．自己組織化マップ入門．Ver. 1.00.00.
神野健二、横田いずみ、井芹慶彦、井料隆太（2009）．自己組織化マップ（Self-Organizing Maps）法と応用例紹介．地下水学会誌、51、151−162．
リスト１ SOMのプログラム例（MainSom.py）：本プログラムはリスト２のモジュールdrawmap.pyを用いている。
import pandas as pd
import numpy as np
import scipy.stats as ss
from drawmap import *
#
# Excelファイルの読み込み
#
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
X = data.values
var_names = data.columns.values
print(var_names)
n = len(X[0]) - 1
JM = len(X)
print('n = ', n)
print('JM = ', JM)
case_labels = X[:,0]
print('case_labels = ', case_labels)
xi = np.empty((JM, n))
# 入力データ
for i in range(JM):
for j in range(n):
xi[i][j] = X[i][j + 1]
print('xi =\n', xi)
g_size = 10
m = np.empty((g_size, g_size, n)) # 参照ベクトル
for i in range(g_size):
for j in range(g_size):
for
k in range(n):
m[i][j][k] = ss.uniform.rvs()
t_max = 100
diffs = np.empty((g_size, g_size))
a_0 = 0.5
sgm_t = g_size / 2
for t in range(t_max):
print('t = ', t)
a_t = a_0 * (t_max - t) /
t_max
for IN in range(JM):
for
i in range(g_size):
for j in range(g_size):
diffs[i][j] = ((m[i][j] - xi[IN]) ** 2).sum()
tc =
np.argmin(diffs) # 勝者ベクトル
tc_i
= tc // g_size
# 勝者ベクトルの行位置
tc_j
= tc % g_size
# 勝者ベクトルの列位置
#
#
参照ベクトルの更新
#
for
i in range(g_size):
for j in range(g_size):
h_cr_IN = a_t * np.exp(-((i - tc_i)**2 + (j - tc_j)**2) /
(2 * (sgm_t**2)))
m[i][j] = m[i][j] + h_cr_IN * (xi[IN] - m[i][j])
sgm_t = 1 + (sgm_t - 1) *
(t_max - t) / t_max
#
# マップの描画
#
draw_som(case_labels, xi, m)
import matplotlib.pyplot as plt
import numpy as np
def draw_som(case_labels, xi, m):
JM = len(xi)
n = len(xi[0])
g_size = len(m)
#
# 参照ベクトルの傾き
#
slopes = np.empty((g_size,
g_size))
for i in range(0, g_size):
for
j in range(0, g_size):
count = 0
v = 0.0
for ix in range(-1, 2):
for jy in range(-1, 2):
try:
if ix != 0 or jy != 0:
if (i+ix) >= 0 and (j+jy) >= 0:
v += np.absolute(m[i+ix][j+jy] - m[i][j]).sum()
count += 1
except:
pass
slopes[i][j] = v / count
slope_max =
np.ma.max(slopes[0:g_size, 0:g_size])
slope_min =
np.ma.min(slopes[0:g_size, 0:g_size])
# 日本語表記の設定
plt.rcParams['font.family']
= ['Yu Gothic', 'Noto Sans CJK JP',
'Hiragino Maru Gothic Pro']
#
# 傾きに応じた色付け
#
plt.figure(figsize = (7, 7))
plt.axis('off')
plt.title('自己組織化マップ\nSOM', fontsize = 20)
for i in range(0, g_size):
for j in range(0, g_size):
vc = (slopes[i][j] - slope_min) / (slope_max - slope_min)
plt.fill_between([i-0.5, i+0.5], [j-0.5, j-0.5], [j+0.5, j+0.5],
color = (vc, 1.0-vc, 0.0), alpha = 0.5)
m_c = []
diffs = np.empty((g_size,
g_size))
plt.xticks([-1, g_size+1],
[' ', ' '])
plt.yticks([-1, g_size+1],
[' ', ' '])
for IN in range(JM):
#
#
勝者ベクトル
#
for
i in range(g_size):
for j in range(g_size):
diffs[i][j] = ((m[i][j] - xi[IN]) ** 2).sum()
tc =
np.argmin(diffs)
tc_i
= tc // g_size
tc_j
= tc % g_size
#
#
勝者ベクトルの位置にプロット
#
plt.text(tc_i, tc_j, case_labels[IN])
plt.show()
リスト３ 2個のニューロンに対するSOMを求めるプログラム（MainTwoClusters.py）。このプログラムはリスト３に示す描画用モジュールmap2clusters.pyを用いている。
import pandas as pd
import numpy as np
import scipy.stats as ss
from map2clusters import *
#
# Excelファイルの読み込み
#
flnm = input('Data File Name (*.xlsx)
= ')
xlsx = pd.ExcelFile(flnm)
data = pd.read_excel(xlsx)
X = data.values
var_names = data.columns.values
print(var_names)
n = len(X[0]) - 1
JM = len(X)
print('n = ', n)
print('JM = ', JM)
case_labels = X[:,0]
print('case_labels = ', case_labels)
xi = np.empty((JM, n)) # 入力データ
for i in range(JM):
for j in range(n):
xi[i][j] = X[i][j + 1]
print('xi =\n', xi)
g_size = 2
m = np.empty((g_size, n)) # 参照ベクトル
for i in range(g_size):
for k in range(n):
m[i][k] = ss.uniform.rvs()
#print('m[0] = n', m[0])
t_max = 100
diffs = np.empty(g_size)
a_0 = 0.5
sgm_t = g_size / 2
for t in range(t_max):
print('t = ', t)
a_t = a_0 * (t_max - t) /
t_max
for IN in range(JM):
for
i in range(g_size):
diffs[i] = ((m[i] - xi[IN]) ** 2).sum()
tc =
np.argmin(diffs)
for
i in range(g_size):
h_cr_IN = a_t * np.exp(-((i - tc)**2) /
(2 * (sgm_t**2)))
m[i] = m[i] + h_cr_IN * (xi[IN] - m[i])
sgm_t = 1 + (sgm_t - 1) *
(t_max - t) / t_max
draw_som_2clusters(case_labels, xi, m)
import matplotlib.pyplot as plt
import numpy as np
def draw_som_2clusters(case_labels,
xi, m):
JM = len(xi)
n = len(xi[0])
g_size = len(m)
diffs0 = []
diffs1 = []
for i in range(JM):
v =
((m[0] - xi[i]) ** 2).sum()
diffs0.append(v)
v =
((m[1] - xi[i]) ** 2).sum()
diffs1.append(v)
plt.plot(diffs0, diffs1,
'bo')
for j in range(JM):
plt.text(diffs0[j], diffs1[j], case_labels[j])
plt.xlabel('$Dist.^2$ from
Unit-0', fontsize = 12)
plt.ylabel('$Dist.^2$ from
Unit-1', fontsize = 12)
plt.show()