2017-06-24

pytorchでgmmの最尤推定

はじめに

今まではKerasを使っていたけど、最近になってpytorchを覚えようとしている。 “Define by Run"と"Define and Run"の違いとかはよくわかっていないのでそのへんは適当。

普通にtutorialだけやっていると、 “なんとかネットワークは作れるけど、自分が考えた新しい層を追加できない” ということになりそうだったので、ネットにあまり情報のなかったgmmを勾配法(最尤推定)で解くプログラムを作って、pytorchを理解することにした。

gaussian mixture model

適当にデータを作る

%matplotlib inline
import pylab as plt
import seaborn as sns
sns.set_style("white")
from scipy.stats import norm

import numpy as np
import torch
from torch.autograd import Variable


K = 2
nb_data = 2000
nb_steps = 2000

true_μ_K = [-3, 3]
true_σ_K = [1, 1]
true_π_K = [0.5, 0.5]

np_data = []
for μₖ, σₖ, πₖ in zip(true_μ_K, true_σ_K, true_π_K):
    for i in range(int(nb_data * πₖ)):
        np_data.append(np.random.normal(μₖ, σₖ))
np_data = np.array(np_data)

def gmm_plot(list_mean, list_std, list_pi, **kwargs):
    x = np.linspace(-10,10,500)
    y = np.sum([norm.pdf(x, mean, np.abs(std))*pi 
                for mean,std,pi in zip(list_mean, list_std, list_pi)], axis=0)
    return plt.plot(x,y, **kwargs)
gmm_plot(true_μ_K, true_σ_K, true_π_K)

f:id:ksknw:20170624233221p:plain

最尤推定

勾配法でパラメータを推定するための誤差関数として、今回は単純な尤度を使った。

gmmの尤度は ${ \displaystyle p(x|\theta) = \prod_n \prod_k \pi_k \mathcal{N}(x | \mu_k, \sigma_k) }$ 。これと、πの合計が1になるように、適当に制約を加えて、以下のように誤差関数を定義した。

def get_normal_lpdf(x_N, μ, σ):
    μ_N = μ.expand(x_N.size())
    σ_N = σ.expand(x_N.size())
    return -0.5 * torch.log(2 * np.pi * σ_N ** 2) - 0.5 * (x_N - μ_N)**2 / σ_N ** 2

def get_loss(normal_lpdf_K_N, π_K):
    gmm_lpdf_N = 0
    for normal_lpdfₖ_N, πₖ in zip(normal_lpdf_K_N, π_K):
        πₖ_N = πₖ.expand(normal_lpdfₖ_N.size())
        gmm_lpdf_N += (torch.exp(normal_lpdfₖ_N) * πₖ_N) # TODO logsumexpを実装したほうがいいかも
    gmm_lpdf = torch.mean(torch.log(gmm_lpdf_N))
    
    Σπ = torch.sum(π_K)
    gmm_lpdf -= torch.abs(1 - Σπ) # 制約条件
    return -gmm_lpdf

pytorchではautograd.Variableで変数を定義しておくと、勝手に微分を計算してくれるらしい。ので、以下のようにデータと求めたいパラメータを定義する。

x_N = Variable(torch.from_numpy(np_data), requires_grad=False).float()

lr = 0.05
μ_K = Variable(torch.randn(K), requires_grad=True)
σ_K = Variable(torch.randn(K)**2, requires_grad=True)
π_K = Variable(torch.abs(torch.randn(K)), requires_grad=True)

あとは以下のように誤差を伝搬させて、パラメータを更新する。 grad.zero_をしないといけないと知らなくて苦労した。

history_loss = []
history_μ_K = []
history_σ_K = []
history_π_K  = []

for i in range(nb_steps):
    normal_lpdf_K_N = []
    for k in range(K):
        normal_lpdf_K_N.append(get_normal_lpdf(x_N, μ_K[k], σ_K[k]))
    loss = get_loss(normal_lpdf_K_N, π_K)
    loss.backward()
    
    μ_K.data -= μ_K.grad.data * lr
    μ_K.grad.data.zero_()
    
    σ_K.data -= σ_K.grad.data * lr
    σ_K.grad.data.zero_()
    
    π_K.data -= π_K.grad.data * lr
    π_K.grad.data.zero_()
    π_K.data = torch.abs(π_K.data)
    
    
    history_loss.extend(loss.data.tolist())
    history_μ_K.append(μ_K.data.tolist())
    history_σ_K.append(σ_K.data.tolist())
    history_π_K.append(π_K.data.tolist())

うまく収束してくれた。

gmm_plot(μ_K.data, σ_K.data, π_K.data)
gmm_plot(true_μ_K, true_σ_K, true_π_K)
plt.show()

f:id:ksknw:20170624233252p:plain

plt.plot(history_loss)
plt.show()

f:id:ksknw:20170624233305p:plain

収束のアニメーションをかいてみる

import matplotlib.animation as animation

plts = []
fig = plt.figure()

for μ_K, σ_K, π_K in zip(history_μ_K[::10],
                         history_σ_K[::10],
                         history_π_K[::10]):   
    plts.append(gmm_plot(μ_K, σ_K, π_K, c="b"))

ani = animation.ArtistAnimation(fig, plts, interval=100)
ani.save('anim.gif', writer="imagemagick")
ani.save('anim.mp4', writer="ffmpeg")
plt.show()

f:id:ksknw:20170624233331g:plain

これはうまく収束した結果だけど、何度か実行すると以下のような微妙な結果に収束することも多かった。 (収束したあとに"脈打ってる"のはなんだ…)

f:id:ksknw:20170624233734p:plain

まとめ

pytorchの使い方を覚えるために、gmmをやった。あまり安定していないけど、動くものができた。なんとなく使い方はわかってきたけど、"勝手に中身が書き換えられている"という印象をもってしまう部分がある。

やる気があればMAP推定、SGDもしくはHMCに続くかも。

2017-04-23

関係データ学習の実装 PythonでStochastic Block Modelの実装

概要

関係データ学習を見ながらpythonでSBMの実装をした。
Twitterのフォローフォロワー関係データに適用して、それっぽい結果を得た。はてなの数式がいまいちわからないので、外部でレンダリングをしていて表示が遅い。

はじめに

ksknw.hatenablog.com

以前StanでStochastic Block Modelをやろうとして失敗した．すっかり忘れていたけど，ふと思い出したので，Stanではなく，普通にpythonで実装することにした．更新式の導出などは前回と同様に「関係データ学習」を参考にした．

データ

前回と同じツイッターのフォローフォロワー関係のデータを使って、アルゴリズムをテストする。使うのは以下のような、非対称な関係データ。

import pandas as pd

data = pd.read_csv("./combinationTable.csv")
uname = data[:1].get_values()[0]
data.drop(0).head()

	11213962	1603589724	68746721	267765193	10985942	14009672	167346791	2896013873	17364190	31442147	...	118320586	218493756	232520574	385409365	419425806	102227818	301210136	175163526	252996913	152543735
1	False	False	True	False	False	False	False	True	False	True	...	False	False	False	False	False	False	False	False	False	False
2	False	False	False	False	False	False	False	False	False	True	...	False	False	False	False	False	False	False	False	False	False
3	False	False	False	False	False	True	False	False	False	False	...	False	False	False	False	False	False	False	False	False	False
4	False	False	False	False	False	False	False	False	False	True	...	False	False	False	False	False	False	False	False	False	False
5	False	False	True	False	False	True	True	False	False	False	...	False	False	False	False	False	False	False	False	False	False

5 rows × 120 columns

%matplotlib inline
import pylab as plt
import seaborn
import numpy as np

X = (data.get_values()[1:]=="True")

def plot_matrix(matrix, z1=None, z2=None):
    f, ax = plt.subplots(figsize=(10, 10))
    plt.pcolor(matrix, cmap=plt.cm.Blues )
    
    if not z1 is None:
        z1_diff = np.r_[[0], np.diff(z1)]
        z2_diff = np.r_[[0], np.diff(z2)]

        for i,c in enumerate(z1_diff):
            if c!=0:
                ax.axhline(i, c="grey")#, linewidth=1)
        for i,c in enumerate(z2_diff):
            if c!=0:
                ax.axvline(i, c="grey") #, linewidth=1)
    plt.show()
    
plot_matrix(X)

f:id:ksknw:20170423192850p:plain

SBMの事後確率

本の内容を参考に、周辺化ギブスサンプラーによって、クラスタの割り当て $z_{1,i}$ , $z_{2,j}$ をサンプリングする。

他の変数がgivenだとした時の、 $z_{1,i}$ の事後確率は、

$p(z_{1,i}=k|{\bf X},{\bf Z}_1^{\backslash(i)}}, {\bf Z}_2)$ $\varpropto \prod_{l=1}^L \hat{\alpha}_{1,k} \frac{\Gamma(\hat{a}_{k,l}+\hat{b}_{k,l})}{\Gamma(\hat{a}_{k,l})\Gamma(\hat{b}_{k,l})} \frac{\Gamma\left(\hat{a}_{k,l}+\sum_{j=1}^{N_2}x_{i,j}\mathbb{I}(z_{2,j}=l)\right)\Gamma\left(\hat{b}_{k,l}+\sum_{j=1}^{N_2}(1-x_{i,j})\mathbb{I}(z_{2,j}=l)\right)}{\Gamma(\hat{a}_{k,l}+\hat{b}_{k,l}+\sum_{j=1}^{N_2}\mathbb{I}(z_{2,j}=l))}$ ここで、

$\hat{\alpha}_{1,k}=\alpha_{1,k}+\hat{m}_{1,k}$

$\hat{a}_{k,l}=a_0 + \hat{n}_{k,l}^{(+)}$

$\hat{b}_{k,l}=b_0 + \hat{n}_{k,l}^{(-)}$

$\hat{m}_{1,k} = \sum_{i'\neqi, i'=1}^{N_1}\mathbb{I}(z_{1,i'}=k)$

$\hat{n}_{k,l}^{(+)} = \sum_{i'\neqi, i'=1}^{N_1}\sum_{j=1}^{N_2}x_{i',j}\mathbb{I}(z_{1,i'}=k)\mathbb{I}(z_{2,j}=l)$

$\hat{n}_{k,l}^{(-)} = \sum_{i'\neqi, i'=1}^{N_1}\sum_{j=1}^{N_2}(1-x_{i',j})\mathbb{I}(z_{1,i'}=k)\mathbb{I}(z_{2,j}=l)$

$\Gamma$ はガンマ関数で、 $K,L,a_0,b_0,\alpha_{1,k}$ はパラメータ

$z_{2,j}$ については対称なので省略。

以上をそのままpythonのプログラムにした。 z1とz2はまとめられそうだけど(というかXを転置して同じ関数に突っ込むだけでいいはずだけど)、ややこしくなるのが嫌だったので2つバラバラの関数として実装した。プログラムは以下。

サンプリングの回数や独立なサンプルを得るためにサンプルを何回おきに保存するかなど、よくわからなかったので、適当に決めた。

実行結果

core-i5-7200U(2.5GHz)でだいたい5時間半かかった。一切並列に計算していないので、むちゃくちゃ重い。ラベルがスイッチしないようにするのが難しいかもしれないが、burn-inが終わった後から、chainを生成して並列にサンプリングしたりしてもいいのかもしれない。

import pickle
with open("./sample_z.pkl", "rb") as f:
    samples_z1, samples_z2 = pickle.load(f)
nb_k = 8

def onehot(i, nb_k):
    ret = np.zeros(nb_k)
    ret[i] = 1
    return ret

z1 = np.array([[onehot(i, nb_k) for i in sample] for sample in samples_z1])
z2 = np.array([[onehot(i, nb_k) for i in sample] for sample in samples_z2])

ちゃんと収束しているかを確認するためにいくつかのヒストグラムを書く。

%matplotlib inline
import pylab as plt

N=6
i=0
plt.hist(np.array(samples_z1)[:,N*i:N*(i+1)], linewidth=0)

([array([   0.,    0.,    0.,    0.,    0.,   38.,    0.,    0.,  953.,    0.]),
 array([ 555.,  430.,    0.,    0.,    0.,    6.,    0.,    0.,    0.,    0.]),
 array([   0.,  988.,    0.,    0.,    0.,    3.,    0.,    0.,    0.,    0.]),
 array([ 971.,   17.,    0.,    0.,    0.,    2.,    0.,    0.,    0.,    1.]),
 array([   0.,  988.,    3.,    0.,    0.,    0.,    0.,    0.,    0.,    0.]),
 array([   0.,  964.,   27.,    0.,    0.,    0.,    0.,    0.,    0.,    0.])],
 array([ 0. ,  0.7,  1.4,  2.1,  2.8,  3.5,  4.2,  4.9,  5.6,  6.3,  7. ]),
<a list of 6 Lists of Patches objects>)

f:id:ksknw:20170423192908p:plain

怪しい部分(緑)もある。収束していないのかもしれない。どのようにクラスタができたかをプロットする。

z1 = z1.mean(axis=0)
z2 = z2.mean(axis=0)

def sort_by_cluster(matrix, z1, z2):
    sorted_mat = list(zip(z1, matrix))
    sorted_mat.sort(key=lambda x:x[0])
    sorted_z1,sorted_mat = zip(*sorted_mat)
    
    sorted_mat = list(zip(z2, np.array(sorted_mat).T))
    sorted_mat.sort(key=lambda x:x[0])
    sorted_z2,sorted_mat = zip(*sorted_mat)
    return np.array(sorted_mat).T, sorted_z1, sorted_z2

%matplotlib inline
import pylab as plt
import seaborn
import numpy as np

X = (data.get_values()[1:]=="True")

plot_matrix(*sort_by_cluster(X, np.argmax(z1, axis=1), np.argmax(z2, axis=1)))

f:id:ksknw:20170423192917p:plain

ここには載せないけど、以下のようにユーザ名とクラスタを確認した。
概ねどっちも大丈夫に思うけれど、z1(フォローする人でクラスタリング)のほうは直感とやや異なる部分もあった。

temp = list(zip(uname, np.argmax(z1, axis=1)))
temp.sort(key=lambda x:x[1])
#temp

temp = list(zip(uname, np.argmax(z2, axis=1)))
temp.sort(key=lambda x:x[1])
#temp

おわりに

PythonでSBMを実装した。
概ねできていると思うが、収束しているのか、そもそもサンプルがちゃんと独立になっているかなど不安な面もいくつかある。

自分で実装を書けばStan(もしくはEdward)でも書けるようになるかなと思ったけど、今のところよくわからない。
事後確率まで書いてサンプリングだけ頼っても意味ない気がするし。

今回はギリシャ文字とか使いまくってプログラムを書いてみた。
数式とできるだけ同じように書くとわかりやすくていいと思う。
EmacsでTeX書式で入力できるようにしたり、半角にしたりすると、特に違和感もなかった。
ただ、デバッグしたり、別の環境でソースコードを見るときは大変かもしれない。
またpythonだと₁とか∇とかを使えないので、z1とか統一できない部分もいくつかあって微妙だった。

参考

機械学習プロフェッショナルシリーズ | 書籍情報 | 株式会社講談社サイエンティフィク

2017-03-26

DTW (動的時間伸縮法)の実装

概要

自分の勉強のために、Dynamic Time Warpingを実装した。正弦波データでいろいろプロットした後、気温のデータに適用した。たぶんバグってないと思う。

はじめに

時系列データの類似度(距離)を計算するとき、単純には例えば各時刻での二乗誤差の平均などを求めることを思いつくが、これは以下のようなデータに対して、直感に反する結果になる。

%matplotlib inline
import numpy as np
import pylab as plt
import seaborn as sns

T = 150
t = .4

A = np.sin(np.array(range(T))/10)
B = np.sin((np.array(range(T))/10 + t*np.pi))
C = np.zeros((T))

plt.plot(A)
plt.plot(B)
plt.plot(C)
plt.show()

f:id:ksknw:20181001220748p:plain

これら3つの時系列データ間の平均二乗誤差を求めると以下のようになる。

mse = lambda a,b: ((a-b)**2).mean()
print("A-B: ", mse(A,B))
print("A-C: ", mse(A,C))

A-B:  0.663947047674
A-C:  0.515002685544

この結果から、青とオレンジよりも青と緑のほうが近いということになる。これは直感に反するので、もう少し上手く距離(的なもの)を定義したい。

Dynamic Time Warping

Dynamic Time Warping (DTW)は2つの時系列データ間の距離(的なもの)を求めるアルゴリズム。
具体的には2つ時系列データの時刻をt1, t2とすると、t1,t2全ての組の誤差を計算し、それらの合計が最小になるような経路を求めるアルゴリズムになっている。
こちらのブログにアニメーションがあって良かった。

以下では、これを参考にして実装する。書いてあったアルゴリズムをできるだけそのまま書こうと思ったのでギリシャ文字とかを使っている。

δ = lambda a,b: (a - b)**2
first = lambda x: x[0]
second = lambda x: x[1]

def minVal(v1, v2, v3):
    if first(v1) <= min(first(v2), first(v3)):
        return v1, 0
    elif first(v2) <= first(v3):
        return v2, 1
    else:
        return v3, 2 

def calc_dtw(A, B):
    S = len(A)
    T = len(B)

    m = [[0 for j in range(T)] for i in range(S)]
    m[0][0] = (δ(A[0],B[0]), (-1,-1))
    for i in range(1,S):
        m[i][0] = (m[i-1][0][0] + δ(A[i], B[0]), (i-1,0))
    for j in range(1,T):
        m[0][j] = (m[0][j-1][0] + δ(A[0], B[j]), (0,j-1))

    for i in range(1,S):
        for j in range(1,T):
            minimum, index = minVal(m[i-1][j], m[i][j-1], m[i-1][j-1])
            indexes = [(i-1,j), (i,j-1), (i-1,j-1)]
            m[i][j] = (first(minimum)+δ(A[i], B[j]), indexes[index])
    return m

print("A-B: ", calc_dtw(A, B)[-1][-1][0])
print("A-C: ", calc_dtw(A, C)[-1][-1][0])

A-B:  9.67371629972
A-C:  77.2504028315

無事にA-Bのほうがだいぶ近そうという結果が得られた。

なにをやっているのかをプロットを書いて確かめる。

def backward(m):
    path = []
    path.append([len(m)-1, len(m[0])-1])
    while True:
        path.append(m[path[-1][0]][path[-1][1]][1])
        if path[-1]==(0,0):
            break
    path = np.array(path)
    return path

import matplotlib.gridspec as gridspec

def plot_path(path, A, B):
    gs = gridspec.GridSpec(2, 2,
                       width_ratios=[1,5],
                       height_ratios=[5,1]
                       )
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])
    ax4 = plt.subplot(gs[3])
    
    list_δ = [[t[0] for t in row] for row in m]
    list_δ = np.array(list_δ)
    ax2.pcolor(list_δ, cmap=plt.cm.Blues)
    ax2.plot(path[:,1], path[:,0], c="C3")
    
    ax1.plot(A, range(len(A)))
    ax1.invert_xaxis()
    ax4.plot(B, c="C1")
    plt.show()
    
    for line in path:
        plt.plot(line, [A[line[0]], B[line[1]]], linewidth=0.2, c="gray")
    plt.plot(A)
    plt.plot(B)
    plt.show()

m = calc_dtw(A, B)
path = backward(m)
plot_path(path, A, B)

f:id:ksknw:20181001220942p:plain

f:id:ksknw:20181001221013p:plain

上のグラフは、t1,t2の位置での蓄積した誤差(青)と選択した経路(赤)を表している。横と下のは元の時系列データ。
また、下のグラフは選択した経路によって、2つの時系列データがどのように対応付けられたかを示している。DTWは結局、この対応する2点の(縦方向の)誤差を求めている。最初と最後は対応づけができないので、無理に合わせるしかない。

m = calc_dtw(A, C)
path = backward(m)
plot_path(path, A, C)

f:id:ksknw:20181001221044p:plain

f:id:ksknw:20181001221055p:plain

どこに対応づいていても同じだが、実装の関係からこんな感じになっていた。

長さが異なる時系列の距離(的なもの)

DTWは2つの時系列間の対応付けを重複を許す形で行うので、長さの異なる時系列データに対しても問題なく適用できる。

T = 150
S = 200
t = .4

A = np.sin(np.array(range(T))/10)
D = np.sin((np.array(range(S))/10 + t*np.pi))

m = calc_dtw(A, D)
path = backward(m)
plot_path(path, A, D)

f:id:ksknw:20181001221210p:plain

f:id:ksknw:20181001221217p:plain

位相をだんだんずらしていったときの挙動

2つの正弦波の位相を徐々にずらしていったときに、MSEおよびDTWが距離(的なもの)をどのように評価するのかをみてみる。

list_dtw_AB = []
list_mse_AB = []
list_dtw_AC = []
list_mse_AC = []
for t in np.arange(0, 2, 0.1):
    T = 150
    A = np.sin(np.array(range(T))/10)
    B = np.sin((np.array(range(T))/10 + t*np.pi))
    list_dtw_AB.append(calc_dtw(A, B)[-1][-1][0])
    list_mse_AB.append(mse(A,B))
    list_dtw_AC.append(calc_dtw(A, C)[-1][-1][0])
    list_mse_AC.append(mse(A,C))

plt.plot(list_mse_AB, label="mse_AB")
plt.plot(list_mse_AC, label="mse_AC")
plt.legend()
plt.show()

plt.plot(list_dtw_AB, label="dtw_AB")
plt.plot(list_dtw_AC, label="dtw_AC")
plt.legend()
plt.show()

f:id:ksknw:20181001221302p:plain

f:id:ksknw:20181001221315p:plain

はじめの例で見たようにMSEは位相が少しずれるとy=0のほうが近いと評価してしまう。
DTWも、おそらく端の対応付けから、値は変化するが、y=0のほうが近いと評価されることはない。
(なんで尖ってるんだ...?)

実データ（国内１週間の気温データ）

正弦波ばっかりいじっててもあれなので、気象庁のサイトから適当な地点の1週間分の気温データを落として、各2点間のDTWを求める。
はじめは北海道と沖縄も入れてやったけど、その2つだけ他と離れすぎていて、グラフが見にくかったのでやめた。

また、本当は海外の「同じ時刻」のデータを持ってきて、時差があるけど云々をやりたかったが、めんどくさくなってやめた。

import pandas  as pd

data = pd.read_csv("./data2.csv", encoding='SHIFT-JIS',skiprows=[0,1, 3,4])
data.index = data["Unnamed: 0"]
data = data[data.columns[1::3]]
print(data.columns.get_values())
data.head()

['名古屋' '豊中' '八戸' '鹿児島' '福岡' '秋田' '盛岡' '山形' '仙台' '福島' 'つくば（館野）' '宇都宮' 'さいたま'
 '千葉' '東京' '新潟' '前橋' '河口湖' '甲府' '静岡' '横浜' '軽井沢' '岐阜' '富山' '金沢' '津' '大津'
 '敦賀' '京都' '奈良' '和歌山' '神戸' '岡山' '鳥取' '松江' '呉' '山口' '松山' '高知' '徳島' '高松' '大分'
 '宮崎' '熊本' '佐賀' '佐世保']

	名古屋	豊中	八戸	鹿児島	福岡	秋田	盛岡	山形	仙台	福島	...	山口	松山	高知	徳島	高松	大分	宮崎	熊本	佐賀	佐世保
Unnamed: 0
2017/3/18 1:00:00	6.4	3.8	-0.4	9.8	7.8	1.4	-0.1	2.0	3.6	4.3	...	3.4	6.1	6.6	7.2	5.4	6.2	6.8	7.8	6.4	7.5
2017/3/18 2:00:00	5.3	3.4	-0.4	9.5	7.0	0.9	0.0	2.2	3.3	4.7	...	3.0	5.6	6.0	6.4	4.7	5.9	6.8	6.1	5.9	6.8
2017/3/18 3:00:00	5.7	2.7	0.0	9.0	7.0	0.2	-0.6	2.0	3.1	5.0	...	2.3	5.4	5.4	6.5	4.0	4.9	6.0	5.8	4.9	6.2
2017/3/18 4:00:00	5.2	2.1	0.5	8.4	6.3	0.2	-1.3	1.5	3.0	5.8	...	1.9	5.3	5.2	6.0	4.1	4.2	5.6	5.3	4.9	6.1
2017/3/18 5:00:00	4.5	2.0	0.8	9.1	6.1	-0.2	-0.6	0.3	3.8	5.7	...	1.7	4.9	4.8	4.5	3.0	3.8	5.1	4.5	4.2	5.7

5 rows × 46 columns

import matplotlib.font_manager as fm
font = {'family' : 'TakaoGothic'}
plt.rc('font', **font)

data.plot(figsize=(15,6), legend=False)
plt.show()

f:id:ksknw:20170326233352p:plain

dtw = []
from itertools import product
for i,(key1, key2) in enumerate(product(data.columns, data.columns)):
    dtw.append(calc_dtw(data[key1].get_values(),
                        data[key2].get_values())[-1][-1][0])

dtw = np.array(dtw).reshape(len(data.columns), -1)

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
ax.set_aspect("equal")
plt.pcolor(dtw, cmap=plt.cm.Blues)
plt.xlim(xmax=46)
plt.ylim(ymax=46)
plt.xticks(np.arange(0.5, len(data.columns)+0.5), list(data.columns), rotation=90)
plt.yticks(np.arange(0.5, len(data.columns)+0.5), list(data.columns))
plt.colorbar()
plt.show()

f:id:ksknw:20170326233404p:plain

軽井沢、河口湖は他の地点と結構違う。リゾート地ということなんだろうか。

なんとなく東北と九州は気候が結構違うのかなと思った。

おまけ

こういう行列を見ると、ここから色々できそうな気がしてしまう。

tSNE

距離行列あるならtSNEできる気がするので、scikit-learnに突っ込んでみる。

from sklearn.manifold import TSNE
tsne = TSNE(metric="precomputed")

tsned = tsne.fit_transform(dtw)

plt.scatter(tsned[:,0], tsned[:,1])
for d,u in zip(tsned, data.columns):
    plt.text(d[0], d[1], u)

f:id:ksknw:20170326233419p:plain

思ってたのと違う図が得られたけど、なんか間違ってる??

スペクトルクラスタリング

類似度という関係データ、しかも対称行列だということでスペクトルクラスタリングをやる。以前書いたものからコードをコピペしてくる。(みんなノートブックのコードの再利用ってどうしてんの)

import scipy.linalg
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA

def plot_matrix(matrix, clusters=None):
    f, ax = plt.subplots(figsize=(10, 10))
    plt.pcolor(matrix, cmap=plt.cm.Blues )

    if not clusters is None:
        clusters_diff = np.r_[[0], np.diff(clusters)]
        for i,c in enumerate(clusters_diff):
            if c==1:
                ax.axhline(i, c="grey") #, linewidth=1)
                ax.axvline(i, c="grey")#, linewidth=1)
    plt.show()
    
def kmeans(eigen_vector, unames, n_clusters):
    kmean = KMeans(n_clusters=n_clusters)
    clusters = kmean.fit_predict(eigen_vector)
    cluster_uname = list(zip(clusters, unames))
    cluster_uname.sort(key=lambda x:x[0])
    return clusters

D = np.eye(dtw.shape[0]) * dtw.sum(axis=0)
X = np.mat(dtw)
K = 5
hi = K-1
lo = 0

L = D - X
L_rw = np.dot(np.linalg.inv(D) , L)  # D**(-1) * L 
eigen_value, eigen_vector = scipy.linalg.eigh(L_rw,eigvals=(lo,hi))

clusters = kmeans(eigen_vector, data.columns, K)

decomp = PCA()
decomped = decomp.fit_transform(eigen_vector)

plt.figure(figsize=(7,7))
plt.scatter(decomped[:,0], decomped[:,1], c=clusters, linewidths=0, cmap=plt.cm.jet)
for d,u in zip(decomped, data.columns):
    plt.text(d[0], d[1], u)
plt.show()

f:id:ksknw:20170326233430p:plain

近い地域が同じクラスタになっているかと思ったら、そうでもない部分もあったり。

この１週間は似てたのかもしれないし、本当はデータを見ないとなぁと思いつつ、面倒なので終わり。

おわりに

自分の勉強のためにDTWの実装をやって、色々なデータを突っ込んでみた。今回はスペクトルクラスタリングをやったけど、参考にした論文やブログでは、重心を求めてkmeansでクラスタリングとかをやっているので、そのあたりもやってみたいなと思う。

参考

2017-03-11

numpyのarrayとmatrixの書き方の違い

はじめに
- 環境
ベクトルと行列ベクトル
トレース，行列式と逆行列
ベクトルの内積
行列とベクトルの積
行列の積
要素ごとの積 (アダマール積)
ベクトルから行列 (テンソル積 ?)
その他，関係ありそうだけど使ってない関数
おわりに

はじめに

numpyの行列関連の演算がわからなくなってググることがよくある．また，動くことは動くがもう少し綺麗にならないかと思うこともよくある．

行列を表すために，numpyではarrayとmatrixを使うことができる．
しかし，掛け算の挙動などが，これら2つで異なるためにさらにややこしい印象がある．

自分用備忘録のためにarray，matrixそれぞれで特定の演算をするためにはどうすればいいかをまとめる．

色々な関数(特にnp.tensordotとか)でこれらの演算はできるが，自分が使うだろうなと思うものだけ書く．

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

環境

import sys
version = sys.version_info
print("Python ",version.major,".",version.minor)
import numpy as np
print("numpy",np.__version__)

Python  3 . 6
numpy 1.11.3

ベクトルと行列ベクトル ${\bf x},A$

以下では，ベクトル $\bf{x}$ と行列 $A$ を使う．
${\bf x} \in \mathbb{R}^2$
$A \in \mathbb{R}^{2\times 2}$

x_a = np.array([i for i in range(1,3)])
x_m = np.mat([i for i in range(1,3)] ).T

A_a = np.array([[i*2+j for i in range(2)] for j in range(2)])
A_m = np.mat([[i*2+j for i in range(2)] for j in range(2)])

"----np.array:----"
A_a
x_a

"----np.matrix:----"
A_m
x_m

'----np.array:----'
array([[0, 2],
       [1, 3]])
array([1, 2])
'----np.matrix:----'
matrix([[0, 2],
        [1, 3]])
matrix([[1],
        [2]])

np.matrix型のベクトルは転置することで列ベクトルにできる．
np.array型もnp.transpose([x_a])と書くと列ベクトルっぽくなるが，後の処理が面倒になるのでしない．

トレース，行列式と逆行列 $\mathrm{trace(A)}, \mathrm{det}(A), A^{-1}$

ありそうだけど，matrix型にdet関数がない．(ないよね?)

"----np.array:----"
np.trace(A_a) 
np.linalg.det(A_a)
np.linalg.inv(A_a) 

"----np.matrix:----"
A_m.trace()
np.linalg.det(A_m)
A_m.I
A_m**-1

'----np.array:----'
3
-2.0
array([[-1.5,  1. ],
       [ 0.5,  0. ]])

'----np.matrix:----'
matrix([[3]])
-2.0
matrix([[-1.5,  1. ],
        [ 0.5,  0. ]])
matrix([[-1.5,  1. ],
        [ 0.5,  0. ]])

ベクトルの内積 $x・x$

array型はnp.dotもしくはx_a.dotで明示的に書かないといけない．

"----np.array:----"
np.dot(x_a, x_a) # arrayはnp.dotを使う．

"----np.matrix:----"
x_m.T * x_m # matrixは掛けるだけでいい

'----np.array:----'
5

'----np.matrix:----'
matrix([[5]])

行列とベクトルの積 $Ax$

上と同じ． matrix型は列ベクトルと行ベクトルがちゃんと結果に反映される．

"----np.array:----"
np.dot(A_a, x_a)
np.dot(x_a ,A_a)

"----np.mat:----"
A_m  * x_m
x_m.T * A_m

'----np.array:----'
array([4, 7])
array([2, 8])

'----np.mat:----'
matrix([[4],
        [7]])
matrix([[2, 8]])

行列の積 $AA, A^\top A$

転置が入ったりすると，arrayではやや面倒になる．
matrixではべき乗はこの積になる．

"----np.array:----"
np.dot(A_a, A_a)
np.dot(np.transpose(A_a), A_a) # 転置が入ると面倒

"----np.mat:----"
A_m  * A_m
A_m**2  # matrixのべき乗はこの積になる
A_m.T * A_m

'----np.array:----'
array([[ 2,  6],
       [ 3, 11]])
array([[ 1,  3],
       [ 3, 13]])

'----np.mat:----'
matrix([[ 2,  6],
        [ 3, 11]])
matrix([[ 2,  6],
        [ 3, 11]])
matrix([[ 1,  3],
        [ 3, 13]])

要素ごとの積 (アダマール積) $C_{ij} = A_{ij}B_{ij}$

array型は*が要素ごとの積になっている．
一方でmatrix型はnp.multiply関数を使う必要がある．
特にべき乗などがある場合には一旦array型に変換するほうが楽かもしれない．

"----np.array:----"
A_a * A_a
A_a ** 3 # arrayのべき乗はこの積

"----np.mat:----"
np.multiply(A_m, A_m) # matrixでは面倒
np.multiply(np.multiply(A_m, A_m), A_m)

np.mat(A_m.A**3) # .Aで一旦arrayに変換したほうが短い

'----np.array:----'
array([[0, 4],
       [1, 9]])
array([[ 0,  8],
       [ 1, 27]])

'----np.mat:----'
matrix([[0, 4],
        [1, 9]])
matrix([[ 0,  8],
        [ 1, 27]])
matrix([[ 0,  8],
        [ 1, 27]])

ベクトルから行列 (テンソル積 ?) $B_{ij} = x_ix_j, B_{ijk} = A_{ij}x_{k}$

(3次元以上の多次元配列のことをテンソルと呼ぶ文化圏にいます．)

ベクトルから行列を作ったり，行列からテンソルを作ったりする．
これはテンソル積なのか?

array型はtensordot(axes=0)を使う． matrix型はベクトル同士だとただ掛けるだけでできる．一方で，出力がテンソルになる計算はできない．

"----np.array:----"
np.tensordot(x_a, x_a, axes=0) # axesによって別の積になる
np.tensordot(x_a, A_a,  axes=0)

"----np.mat:----"
x_m * x_m.T # ベクトルはこれでいける
x_m * A_m # 答えがテンソルになる演算はできない．

'----np.array:----'
array([[1, 2],
       [2, 4]])
array([[[0, 2],
        [1, 3]],
       [[0, 4],
        [2, 6]]])

'----np.mat:----'
matrix([[1, 2],
        [2, 4]])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-ede3e147096b> in <module>()
      5 "----np.mat:----"
      6 x_m * x_m.T
----> 7 x_m * A_m
/home/kei/anaconda3/envs/tensorflow/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py in __mul__(self, other)
    341         if isinstance(other, (N.ndarray, list, tuple)) :
    342             # This promotes 1-D vectors to row vectors
--> 343             return N.dot(self, asmatrix(other))
    344         if isscalar(other) or not hasattr(other, '__rmul__') :
    345             return N.dot(self, other)
ValueError: shapes (2,1) and (2,2) not aligned: 1 (dim 1) != 2 (dim 0)

その他，関係ありそうだけど使ってない関数

なんか色々あるけど，よく知らないので使ってない．

おわりに

記述が楽なので，matrixを使って書いてしまうことが多いけど，3次元以上はarrayで書くことになるので，慣れてないとバグを生む可能性が高い気がする．皆さんどうしてるんですか．

2017-03-05

関係データ学習の実装 MovieLensの行列分解

概要

関係データ学習の後半のうち，行列分解を実装した． lossが下がっているので，たぶん合っていると思うが，分解したあとのベクトルを見るとよくわからない．

はじめに

以前こんなものをやった． ksknw.hatenablog.com

前半部分は終わったので(SBMのことは置いといて)，後半部分の実装をやる．今回は行列分解を実装する．

www.kspub.co.jp

MovieLensデータ

MovieLens 100K Dataset | GroupLens

これといって分解したい行列データもなかったので，本に書いてあったデータを使うことにする． MovieLensにはユーザが映画を評価したデータが含まれている．
データを作成した時期によってデータのサイズが違うらしいが，今回は一番小さいサイズのデータを使う．
このデータには943人のユーザの映画の評価(1から5)が含まれている．
映画は1682本で各ユーザはこれらのうち少なくとも5本を評価している．
とりあえずデータを見てプロットする．

%matplotlib inline
import pylab as plt
import seaborn
import numpy as np
import pandas as pd
data = pd.read_csv("./ml-100k/u.data", names=["uid", "mid", "rating", "timestamp"], sep="\t")
data

	uid	mid	rating	timestamp
0	196	242	3	881250949
1	186	302	3	891717742
2	22	377	1	878887116
3	244	51	2	880606923
4	166	346	1	886397596
29	160	234	5	876861185
...	...	...	...	...
99995	880	476	3	880175444
99996	716	204	5	879795543
99997	276	1090	1	874795795
99998	13	225	2	882399156
99999	12	203	3	879959583

100000 rows × 4 columns

X = np.zeros((data["uid"].max(), data["mid"].max()))
for i,item in data.iterrows():
    X[item["uid"]-1, item["mid"]-1] = item["rating"]

f, ax = plt.subplots(figsize=(7, 7))

import numpy.ma as ma
masked = ma.masked_where((X==0),X)


cmap = plt.cm.jet
cmap.set_bad('white',1.)
plt.pcolormesh(masked, cmap=cmap)
plt.xlabel("movie id")
plt.ylabel("user id")
plt.show()

f, ax = plt.subplots(figsize=(7, 7))

plt.pcolormesh(masked[:100, :100], cmap=cmap)
plt.xlabel("movie id")
plt.ylabel("user id")
plt.show()

f:id:ksknw:20170305215433p:plain

f:id:ksknw:20170305215506p:plain

人気の映画やよく映画を見ている人などがなんとなくわかる．多くの人が好む映画，やたら低評価ばかりつけるユーザなど色々いる．なんらかソートされていることもわかる．
また大半は欠損データである．

まずは，(練習なので)欠損データは0とレーティングされているとして行列分解する．
その後，欠損値を除外したり，補完をするアルゴリズムを実装する．
また，本に載っているアルゴリズムでは収束判定をしていたが省略する．

import copy 
original_X = copy.deepcopy(X) 

# X = X+(X==0)*3
#X = X[:100, :100]
f, ax = plt.subplots(figsize=(7, 7))
plt.pcolormesh(X[:100, :100], cmap=cmap)

f:id:ksknw:20170305215522p:plain

1次交互勾配降下法

正則化つきの行列分解は特異値分解では求まらないので，勾配法を用いる必要がある．

それぞれ5つのベクトルでユーザ，映画が表現できると考えて，行列分解する．
動けばいいやの精神で書いているので，行列の計算はあれな感じだけど許してほしい．というかnumpy.matとarray関連の正しい使い方がわからない．

1次交互勾配法は不安定で，learning_rateを少しあげたら発散した．

nb_epoch = 100
reg = 0.01

R = 5
learning_rate = 0.001

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.random.random((I, R)))
V = np.mat(np.random.random((J, R)))

dU = lambda U,V: -X*V + U * (V.T*V + reg*np.mat(np.identity(R)))
dV = lambda U,V: -X.T*U + V * (U.T*U + reg*np.mat(np.identity(R)))
error = lambda U,V: 0.5*np.sum(np.array(X-U*V.T)**2) + 0.5*reg*np.sum( np.sum(np.array(U)**2,axis=0) + np.sum(np.array(V)**2,axis=0))
err_hist = []

for i in range(nb_epoch):
    e = error(U,V)
    U -= learning_rate * dU(U,V)
    V -= learning_rate * dV(U,V)
    err_hist.append(e)

def hist_plot(err_hist):
    plt.plot(err_hist)
    plt.text(len(err_hist) -5, err_hist[-1]+500 ,"%.1f"%err_hist[-1])
hist_plot(err_hist)

f:id:ksknw:20170305215543p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)

(10.693971705841749, -2.4352853936451466)

def compare_plot(X, predict, cmap=plt.cm.Blues, vmin=None, vmax=None):
    f, ax = plt.subplots(figsize=(14, 7))
    plt.subplot(121)
    plt.pcolormesh(X[:100, :100], cmap=cmap, vmin=vmin, vmax=vmax)
    plt.xlabel("movie id")
    plt.ylabel("user id")
    
    plt.subplot(122)
    plt.pcolormesh(predict[:100, :100], cmap=cmap, vmin=vmin, vmax=vmax)
    plt.xlabel("movie id")
    plt.ylabel("user id")

    plt.show()
compare_plot(X, predict)

f:id:ksknw:20170305215555p:plain

最大値，最小値がずれているのでJetだと見にくかったので，青だけでプロットしている．
なんとなく出来ている気がする．

擬似2次交互勾配法

1次勾配だけだと収束が不安定で学習率をあまり大きくすることができず，収束が遅い．
2次勾配を使うことで収束を早くすることができる．

learning_rate = 1.

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.random.random((I, R)))
V = np.mat(np.random.random((J, R)))
err_hist2 = []

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    V = (1-learning_rate)*V + learning_rate * X.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    err_hist2.append(e)

hist_plot(err_hist)
hist_plot(err_hist2)

f:id:ksknw:20170305215610p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)

(10.414641991236321, -2.2410733139172789)

compare_plot(X, predict)

f:id:ksknw:20170305215631p:plain

学習率を1にしても発散せず，1次勾配(青)よりも2次勾配(緑)を利用したものは早く収束するようになった．

制約付きの最適化

ここまでの実装では，推定した評価値が負の値になっている． U，Vに制約をつけることで0以上の値にすることができる． (評価値は1から5までなので本当はもう少し制約しないといけないがどうするのかはわからない．)

具体的には，負の値をとった要素を0に射影しながら最適化する．

learning_rate = 1.

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.random.random((I, R)))
V = np.mat(np.random.random((J, R)))
err_hist_con = []

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    U = np.maximum(U, 0, U)
    V = (1-learning_rate)*V + learning_rate * X.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    V = np.maximum(V, 0, V)

    err_hist_con.append(e)

hist_plot(err_hist_con)

f:id:ksknw:20170305215653p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)

(8.9972633401464464, 0.0)

ちゃんと0以上に制約されている．

compare_plot(X, predict)

f:id:ksknw:20170305215708p:plain

欠損値をちゃんとやる

ここまでは欠損値を全て0だと仮定して実装した．

以下では欠損値を除外して行列分解する方法と補完する方法を行う．

除外する方法

lossの計算およびU，Vの勾配の計算のときに，欠損した値は利用しない．
実装が楽なので1次勾配を利用する方法で実装した．

M = np.array((X!=0)*1)

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.abs(np.random.random((I, R))))
V = np.mat(np.abs(np.random.random((J, R))))

dU = lambda U,V: -(np.array(X - U*V.T)*M)*V   + reg*U
dV = lambda U,V: -(np.array(X - U*V.T)*M).T*U + reg*V
error = lambda U,V: 0.5*np.sum((np.array(X-U*V.T)*M)**2) + 0.5*reg*np.sum( np.sum(np.array(U)**2,axis=0) + np.sum(np.array(V)**2,axis=0))
err_hist = []

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    V = (1-learning_rate)*V + learning_rate * X.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    err_hist.append(e)
hist_plot(err_hist)

f:id:ksknw:20170305215722p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)

(10.414841938424392, -2.2413625198086238)

import numpy.ma as ma
masked_predict = ma.masked_where((X==0),predict)

compare_plot(masked[:100, :100], masked_predict, cmap=cmap, vmin=1, vmax=5)
compare_plot(masked[:100, :100], predict, cmap=cmap, vmin=1, vmax=5)

f:id:ksknw:20170305215738p:plain

f:id:ksknw:20170305215755p:plain

全体的に評価が低めな感じになってしまっている．また，欠損していた部分も低めの値で埋まっている感じ．

補完する方法

こちらの方法では欠損値を毎回，予測値で補完する．欠損値を除外する方法と等価らしいので，こちらは2次勾配を使って，制約もつける方法で実装する．

M_hat = np.array((X==0)*1)

I = X.shape[0]
J = X.shape[1]
U = np.mat(np.abs(np.random.random((I, R))))
V = np.mat(np.abs(np.random.random((J, R))))
error = lambda U,V: 0.5*np.sum((np.array(X-U*V.T)*M)**2) + 0.5*reg*np.sum( np.sum(np.array(U)**2,axis=0) + np.sum(np.array(V)**2,axis=0))

err_hist_comp = []

X_hat = X + np.mat(np.array(U*V.T)*M_hat)

for i in range(nb_epoch):
    e = error(U,V)
    U = (1-learning_rate)*U + learning_rate * X_hat*V*((V.T*V+reg*np.mat(np.identity(R))).I)
    U = np.maximum(U, 0, U)
    X_hat = X + np.mat(np.array(U*V.T)*M_hat)
    V = (1-learning_rate)*V + learning_rate * X_hat.T*U*((U.T*U+reg*np.mat(np.identity(R))).I)
    V = np.maximum(V, 0, V)
    X_hat = X + np.mat(np.array(U*V.T)*M_hat)
    err_hist_comp.append(e)
hist_plot(err_hist_comp)

f:id:ksknw:20170305215813p:plain

predict = np.array(U*V.T)
np.max(predict), np.min(predict)

(7.7494405681324903, 0.0082254500224862269)

import numpy.ma as ma
masked_predict = ma.masked_where((X==0),predict)

compare_plot(masked[:100, :100], masked_predict, cmap=cmap, vmin=1, vmax=5)
compare_plot(masked[:100, :100], predict, cmap=cmap, vmin=1, vmax=5)

f:id:ksknw:20170305215827p:plain

f:id:ksknw:20170305215841p:plain

だいたいできているような気がする．

ベクトルをみる

行列分解で得られたU，Vでは近いユーザもしくは映画が近いベクトルで表現される．
MovieLensには評価だけでなく，どんな映画かやどんなユーザかといった情報も含まれている．
これらとU,Vを比較することで，正しい行列分解ができたかを評価する．

普段あまり映画を見ないのと，古い映画のデータしかないので，Toy StoryとGolden Eyeぐらいしかわからなかった．
ユーザの属性ならなんとなくわかるので，そちらを使う．ユーザの情報は以下．

user_info = pd.read_csv('ml-100k/u.user', sep='|',
                         encoding='latin-1',header=None)
user_info

	0	1	2	3	4
0	1	24	M	technician	85711
1	2	53	F	other	94043
2	3	23	M	writer	32067
3	4	24	M	technician	43537
...	...	...	...	...	...
937	938	38	F	technician	55038
938	939	26	F	student	33319
939	940	32	M	administrator	02215
940	941	20	M	student	97229
941	942	48	F	librarian	78209
942	943	22	M	student	77841

943 rows × 5 columns

職業がわかりやすそうなので，プロットしてみる．

from sklearn.manifold import TSNE
tsne = TSNE()
tsned = tsne.fit_transform(U) 
plt.scatter(tsned[:,0], tsned[:,1], c=(user_info[3]=="scientist"), cmap=plt.cm.Accent, marker=".")
plt.show()
plt.scatter(tsned[:,0], tsned[:,1], c=(user_info[3]=="engineer"), cmap=plt.cm.Accent, marker=".")
plt.show()
plt.scatter(tsned[:,0], tsned[:,1], c=(user_info[3]=="student"), cmap=plt.cm.Accent, marker=".")
plt.show()

f:id:ksknw:20170305215854p:plain

f:id:ksknw:20170305215903p:plain

f:id:ksknw:20170305215911p:plain

そうでもない感じだったので，うまくいっているのかよくわからない．

まとめ

勉強のために行列分解の実装をした．収束はしたけど，うまくいっているのかはよくわからない．

参考

2017-01-08

EmacsでPythonとTeXとMarkdownを書いてる動画と使ってるパッケージ

はじめに

今年は頑張ってアウトプットを増やすことを目標にした。とはいえ、特に書くこともなかったので、Emacsについて書く。

正月やることがなかったので、またEmacsの設定を見なおした。 Emacs Rocksという動画に触発されて、動画で説明することにした。普段pythonとTeXとMarkdownしか書かないので、それらについての設定しかできていない。

環境

Emacs 24.5
Ubuntu 14.04
.emacs

github.com

python

youtu.be

python modeはこのころとあまり変わっていない。

yasnippet
予め登録されているnpとかmainとかの後にtabを打つと、スニペットを挿入してくれる。自分でスニペットを登録できるので、独自のライブラリに関するスニペットとかも作れて便利。
jedi
jediはオムニ補完、つまり文法的な部分をある程度汲んだ上で補完をしてくれる。動画ではやらなかったが関数定義にジャンプしたりもできる。 importした直後とかはたまに補完候補をだすのに数秒かかることもある。
flymake-python-pyflakes
文法エラーやwarningの表示。 flymake-cursorを入れないとカーソル位置のエラーが表示されなかったようなそうでもなかったような気がする。
py-autopep8
pep8準拠のコードに自動的に修正。演算子の両側にスペースを入れたり、改行を入れたり。 after-save-hookに登録して、保存時に修正するといい感じ。
iedit
複数箇所の同時編集。変数名を変えたいときに便利。範囲を関数内に限定したり、手動で範囲を決めたりできる。 multiple-cursorsとどっちがいいんでしょうか。
smartparens
括弧などの補完。補完もそうだけど、範囲選択して括弧を入力すると、その範囲を括弧で囲ってくれるのがいい。
highlight-parentheses
括弧の色を深さで変える。趣味。

TeX

youtu.be

AUCTeX
以前までYaTeXを使っていたけど、どちらがいいかわからなかったので、とりあえず使ってみている。数式に関するショートカットがなくなったけど、yasnippet使うしどっちでもいい気がする。 C-c C-c でlatexmkを呼ぶようにできるらしいので、その点はこっちのほうがいいかも。
yasnippet
機能はpythonと同じだが、複雑な構文が多いのでかなり役に立つ。 figureの構文とか覚えてられないし。 TeX固有の機能として、cite + Tabとかref+Tabでreftexを呼ぶようになっている。デフォルトのままなので、自分の環境に合わせてもうちょっといじったほうがいいかなと思っている。
reftex
bibファイルは使い回しなので、どういう名前で登録したか覚えていないし、長い文章を書いていると、\labelで設定した名前なんか忘れる。 reftexは参考文献入れるときや\refを使うときに、名前を補完してくれる。参考文献は特に著者名前でもキーワードでもヒットするので便利。
同じbibファイルを使いまわすときは、default-bibliographyに登録しておくと便利。
magic latex buffer
\beginとか\endを▽△に変えたり、数式を表示したり、見やすいように色々やってくれる。コマンドのタイポに気づきやすい。数式以外の場所で"_"を打っても添字扱いされるところが少し困る。

Markdown

youtu.be

OrgTbl
markdownで表を書いてると、自分で合わせない限り縦棒がずれているので見にくい。 orgの表を書くモードをMarkdownモードでも使えるので、それを使う。 orgとMarkdownで一部書式が違うので、保存時に変換している。
flymd
Markdownのリアルタイムプレビューはいくつかパッケージがあるけど、今はflymdを使っている。画像や数式もちゃんと表示されるのでいい。
misc
研究ノート的なものをmarkdownで書いているので、C-x C-nで常に同じファイルが開くようになっている。 org-modeでいいじゃないかという気もしている。

その他

git-gutter-fringeの見た目をちょっといじって使っている。いい感じ。

f:id:ksknw:20170108204246p:plain

今後の課題

未だに改行をEnterでしか打てないので、C-mとかC-jとか使っていきたい。

参考

2016-12-29

関係データ学習の実装ツイッターデータのスペクトルクラスタリングとSBM

概要

関係データ学習の学習のために，自分で実装して理解する．ツイッターのフォローフォロワー関係を使って，グラフラプラシアンを求めスペクトルクラスタリングを行った．その結果，なんとなくクラスタリングできた．また，確率的ブロックモデルによる非対称データクラスタリングをStanによって実装しようとした．これはうまくいっていない．

はじめに

関係データ学習という本を買って読んでいる．

www.kspub.co.jp

本の内容は前半と後半に分かれていて，前半は関係データをスペクトルクラスタリングしたり，確率的ブロックモデルでクラスタリングしたりする話．後半は行列分解やテンソル分解の話になっている．まだ前半の途中までしか読めていないが，予想していたよりも数式が簡単だったこともあり，実際のデータに適用してみたくなった．数年前に書いたツイッターのフォローフォロワー関係をダウンロードするスクリプトを見つけたので，これを使ってデータを作り，自分のフォロー中の人たちをクラスタリングしてみることにした． (カット最小化とスペクトルクラスタリングの関係とか，まだよく理解できてないので，理論的な説明はしない．)

ツイッターデータのダウンロード

tweepyを使って関係データをダウンロードする．スクリプトは何年も前に書いたものであり，なんとなく恥ずかしいのでここでは省略するが，

自分がフォローしている人のフォローしている人を全てダウンロード．
ダウンロードした中から自分がフォローしている人たちの関係を抜き出す．

ようになっている．ツイッターのapi制限があるので，100人ちょっとだけど，それなりに時間がかかる．

ダウンロードしたのは以下のようなデータ。つながっている部分が青、ない部分が白。鍵垢のデータはとれてない部分もある。ツイッターではフォローしているがされていないなどの関係もあるので，非対称データになる．

import pandas as pd

data = pd.read_csv("./combinationTable.csv")
uname = data[:1].get_values()[0]

%matplotlib inline
import pylab as plt
import seaborn
import numpy as np

X = (data.get_values()[1:]=="True")

def plot_matrix(matrix, clusters=None):
    f, ax = plt.subplots(figsize=(10, 10))
    plt.pcolor(matrix, cmap=plt.cm.Blues )

    if not clusters is None:
        clusters_diff = np.r_[[0], np.diff(clusters)]
        print clusters_diff
        for i,c in enumerate(clusters_diff):
            if c==1:
                ax.axhline(i, c="grey") #, linewidth=1)
                ax.axvline(i, c="grey")#, linewidth=1)
    plt.savefig("../../temp.svg")
    plt.show()

plot_matrix(X)

f:id:ksknw:20161225112150p:plain

すでにある程度のクラスタができているように見える．これは自分が今まで所属した(ツイッター的意味での)クラスタを表している．自分がフォローした順にアカウントが並んでいるので，右上からだいたいポケモン，大学，大学研究室，研究者みたいな感じ．研究者クラスタ以外は知り合いをフォローしていることが多いので，数が少ない．

はじめはこの関係データを対称データに変換し，2章を参考にしつつスペクトルクラスタリングを行う．以下のように対称データに変換した．

symmetry_X = (X.T + X)
plot_matrix(symmetry_X)

f:id:ksknw:20161225112224p:plain

非正規化グラフラプラシアンに基づくスペクトルクラスタリング

対称データ $X$ とクラスタ数 $K$ が与えられたとき、非正規化グラフラプラシアン $L$ は以下のように定義される。

$L=D-X$

ここで $D$ は次数行列で、オブジェクトから伸びるリンクの総数。 $L$ をプロットすると以下。

D = np.eye(symmetry_X.shape[0]) * symmetry_X.sum(axis=0)
K = 10
L = D - X
plot_matrix(L)

f:id:ksknw:20161225125154p:plain

先ほどまでとはグラフの色が異なるので注意．対角成分だけ1以上の値をとるのでこんな感じの図になる。

スペクトルクラスタリングでは、以下の手順でクラスタを抽出する。

グラフラプラシアンを固有値分解して固有ベクトルを求める。
固有ベクトルをK-meansなどでクラスタリングする。

クラスタ数は適当に10とした。本来なら色々チューニングしないといけないが、面倒なので固定する。

import scipy.linalg

hi = K-1
lo = 0
eigen_value,eigen_vector = scipy.linalg.eigh(L,eigvals=(lo,hi))

def plot_eigen_vec(eigen_vector):
    plt.figure(figsize=(30,10))
    plt.plot(np.abs(eigen_vector))
    plt.show()

from sklearn.cluster import KMeans
def kmeans(eigen_vector, unames, n_clusters):
    kmean = KMeans(n_clusters=n_clusters)
    clusters = kmean.fit_predict(eigen_vector)
    cluster_uname = zip(clusters, uname)
    cluster_uname.sort(key=lambda x:x[0])
    print cluster_uname
    return clusters

def sort_by_cluster(matrix, clusters):
    sorted_mat = zip(clusters, matrix)
    sorted_mat.sort(key=lambda x:x[0])
    _,sorted_mat = zip(*sorted_mat)
    
    sorted_mat = zip(clusters, np.array(sorted_mat).T)
    sorted_mat.sort(key=lambda x:x[0])
    sorted_clusters,sorted_mat = zip(*sorted_mat)
    return np.array(sorted_mat).T, sorted_clusters

plot_eigen_vec(eigen_vector)
clusters = kmeans(eigen_vector, uname, K)
plot_matrix(*sort_by_cluster(symmetry_X, clusters))
sorted_mat, sorted_clusters = sort_by_cluster(symmetry_X, clusters)

クラスタリング結果は以下。クラスタの番号順に並び替えて、境界に線を引いた。

f:id:ksknw:20161225125552p:plain

でかいクラスタが1つと1アカウントしかないクラスタがいくつかできた。正規化しない単純なグラフラプラシアンだとこういう結果になるらしい。

対称正規化グラフラプラシアンに基づくスペクトルクラスタリング

もう少し望ましいクラスタリングを得るために、グラフラプラシアンの正規化を行う。これは(緩和した)グラフカットの正規化と等価になる(?)らしいが、細かい説明は本(とその参考文献)に譲る。読まねば。

対象正規化グラフラプラシアン $L_\mathrm{sym}$ は、上で定義したグラフラプラシアン $L$ から以下のように求まる。

$L_\mathrm{sym} =D^{-\frac{1}{2}} L D^{-\frac{1}{2}}$

L = D - X
L_sym = np.dot(np.dot(scipy.linalg.sqrtm(np.linalg.inv(D)) , L) , scipy.linalg.sqrtm(np.linalg.inv(D))) # D**(-0.5) * L * D**(-0.5)
plot_matrix(L_sym)
eigen_value,eigen_vector = scipy.linalg.eigh(L_sym,eigvals=(lo,hi))
eigen_vector = eigen_vector/np.linalg.norm(eigen_vector)

対象正規化グラフラプラシアンのプロットは以下。非正規化グラフラプラシアンに比べて対角項が小さくなっている。また、クラスタリングする前に固有ベクトルを正規化している。

f:id:ksknw:20161225144212p:plain

plot_eigen_vec(eigen_vector)
clusters = kmeans(eigen_vector, uname, K)
plot_matrix(*sort_by_cluster(symmetry_X, clusters))

クラスタリング結果は以下。先程に比べればだいぶ良さそうだが、真ん中の巨大なクラスタにほとんど周りとつながっていないものがいくつか含まれている。

f:id:ksknw:20161225144619p:plain

酔歩正規化グラフラプラシアン

別の正規化から導かれるグラフラプラシアンを使ってクラスタリングをする。

$L_\mathrm{rw} =D^{-1} L$

L = D - X
L_rw = np.dot(np.linalg.inv(D) , L)  # D**(-1) * L 

plot_matrix(L_rw)
eigen_value,eigen_vector = scipy.linalg.eigh(L_rw,eigvals=(lo,hi))

plot_eigen_vec(eigen_vector)
clusters = kmeans(eigen_vector, uname, K)
plot_matrix(*sort_by_cluster(symmetry_X, clusters))

f:id:ksknw:20161225145741p:plain

でかいクラスタが少し分割されて、まあこっちのほうがいいかなぁという感じ。

固有ベクトルをPCAで2次元に落としてプロットする。左でポケモンクラスタが孤立している。右に研究室のクラスタがきた。ドクターの人と先生のアカウントは研究者のクラスタとエッジをたくさんもっている。研究室が同じかつ大学の同級生は大学同級生クラスタと研究室クラスタの間に配置された。このような感じで、まあまあ妥当なクラスタリングができた。

from sklearn.decomposition import PCA
# from sklearn.manifold import TSNE
# decomp = TSNE()
decomp = PCA()
decomped = decomp.fit_transform(eigen_vector)

%matplotlib inline
plt.figure(figsize=(20,20))
for i,row in enumerate(X):
    for j,is_connected in enumerate(row):
        if is_connected:
            plt.plot(decomped[[i,j], 0], decomped[[i,j], 1], alpha=0.2,linewidth=0.3, c="gray")
plt.scatter(decomped[:,0], decomped[:,1], c=clusters, linewidths=0, cmap=plt.cm.jet)
#for d,u in zip(decomped, uname):
#    plt.text(d[0], d[1], u)
plt.show()

f:id:ksknw:20161229121731p:plain

確率的ブロックモデルによる非対称データクラスタリング

ここまでは、本来非対称データであるツイッターの関係データを対称データとしてクラスタリングした。ここからは、確率的ブロックモデル(SBM)を使って非対称データのままクラスタリングする。

$i$ 番目のアカウントが $j$ 番目のアカウントをフォローしているかどうかを ${ \displaystyle x_{i,j} }$ と書く。また、 $z_{1,i}$ を $i$ 番目のアカウントが所属するフォローする側のクラスタ、 $z_{2,i}$ をフォローされる側のクラスタとする。フォローする側のクラスタとか意味がわからないかもしれないが、例えば機械学習関連のアカウントをフォローしまくっているクラスタとかそういう感じだと思う。

このときSBMでは以下の生成モデルを仮定する。

$\pi_1 \mid \alpha_1 \sim \mathrm{Dirichlet}(\alpha_1)$
$\pi_2 \mid \alpha_2 \sim \mathrm{Dirichlet}(\alpha_2)$
$z_{1,i} = k\mid \pi_1 \sim \mathrm{Discrete}(\pi_1)$
$z_{2,i} = l\mid \pi_2 \sim \mathrm{Discrete}(\pi_2)$
$\theta_{k,l} \mid a_0,b_0 \sim \mathrm{Beta}(a_0,b_0)$
${ \displaystyle x_{i,j} \mid {\theta_{k,l}} , z_{1,i}, z_{2,j} \sim \mathrm{Bernoulli}(\theta_{z_{1,i}, z_{2,j}}) }$

SBMを使うと、「フォローする(1)」と「フォローされる(2)」にわけてクラスタを考えることになる。例えば「フォロワーするクラスタ」は有名な研究者をフォローするクラスタだが、「フォローされるクラスタ」は研究室の人とかそういうことになる。うっ

SBMでは、クラスタ間の関係の強さ $\theta_{k,l}$ とそれぞれのアカウントがどのクラスタに所属するか $z_{1,i}$ 、 $z_{2,i}$ を同時に推定する。本では周辺化ギブスサンプラーによる推論方法の導出を行っている。ここでは、Stanを用いて、SBMを解くことを試みる。まずは、上の数式をそのままStanに変換すると以下のようになる。

data{
  int N;
  int K;
  int x[N, N];
}

parameters{
  vector<lower=0>[K] alpha1;
  vector<lower=0>[K] alpha2;
  int z1[N];
  int z2[N];
  vector[K] pi1;
  vector[K] pi2;
  real<lower=0> a0;
  real<lower=0> b0;
  matrix[K,K] theta;
}

model{
  pi1 ~ dirichlet(alpha1);
  pi2 ~ dirichlet(alpha2);

  for (i in 1:N){
    z1[i] ~ categorical(pi1);
    z2[i] ~ categorical(pi2);
  }

  for (k in 1:K){
    for(l in 1:K){
      theta[k,l] ~ beta(a0, b0);
    }
  }

  for (i in 1:N){
    for (j in 1:N){
      x[i,j] ~ bernoulli(theta[z1[i], z2[j]]);
    }
  }
}

これには色々問題があるが、特にint型の変数がparametersの中に含まれていることが問題である。 Stanでは離散的なパラメータを扱うことができないので、zで周辺化する必要がある。 zで周辺化したり、遅すぎないようにベクトルにしたりで以下のようになる。~~(たぶんあっていると思う。)~~間違っている。

data{
  int N;
  int K;
  int x[N*N];
}

parameters{
  vector<lower=0>[K] alpha1;
  vector<lower=0>[K] alpha2;
  simplex[K] pi1;
  simplex[K] pi2;
  real<lower=0> a0;
  real<lower=0> b0;
  vector<lower=0,upper=1>[K*K] theta;
}

model{
  real gamma[K*K];

  pi1 ~ dirichlet(alpha1);
  pi2 ~ dirichlet(alpha2);
  theta ~ beta(a0, b0);

  for (k in 1:K){
    for (l in 1:K){
      gamma[k+ K*(l-1)] = log(pi1[k]) + log(pi2[l]) + bernoulli_lpmf(x | theta[k+ K*(l-1)]);
    }
  }
  target += log_sum_exp(gamma);
}

サンプリングする。

import pystan
data_dict= {"N":len(X),
            "K":10,
            "x":X.astype(int).reshape(-1)}
fit = pystan.stan("sbm_vector.stan", data=data_dict)
fit

            mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha1[0]  2.6e7   3.1e7  5.3e7 3429.6 4934.1  3.9e4  1.2e7  1.8e8    3.0    3.2
alpha1[1]  3.2e8   3.5e8  5.0e8  3.0e5  1.4e6  3.7e7  5.4e8  1.4e9    2.0   6.44
alpha1[2]  1.0e7   6.9e6  9.7e6 460.13  1.1e6  6.9e6  1.5e7  3.4e7    2.0    2.7
alpha1[3]  2.9e7   3.6e7  5.0e7 9752.8  2.5e4  7.6e5  4.8e7  1.4e8    2.0   9.99
alpha1[4]  9.3e6   1.2e7  1.6e7   1.29  5.7e4  3.8e5  9.9e6  5.1e7    2.0   5.75
...
theta[97]   0.58    0.14   0.25   0.16    0.3    0.6   0.74   0.98    3.0   3.63
theta[98]   0.29    0.04   0.16   0.02   0.23   0.23   0.35   0.78   17.0   1.12
theta[99]   0.35    0.09   0.16   0.11   0.19   0.34   0.54   0.57    3.0   2.23
lp__       -4015   40.35  57.06  -4105  -4053  -4019  -3975  -3929    2.0  12.14

fit.plot()

f:id:ksknw:20161229165851p:plain

うまく収束しなかった。特に対策をしていないために、ラベルスイッチングが起こっているのが原因かなと思ったので、 $\theta$ の片側をorderに変えた。 0~1の制約も残さないといけないので、すこしごちゃごちゃした。

data{
  int N;
  int K;
  int x[N*N];
}

parameters{
  vector<lower=0>[K] alpha1;
  vector<lower=0>[K] alpha2;
  simplex[K] pi1;
  simplex[K] pi2;
  real mu;
  real<lower=0> sigma;
  ordered[K] theta_l[K];
}

transformed parameters{
  vector<lower=0, upper=1>[K] theta[K];
  for (k in 1:K){
    for (l in 1:K){
      theta[k][l] = inv_logit(theta_l[k][l]);
    }
  }
}

model{
  real gamma[K*K];
  pi1 ~ dirichlet(alpha1);
  pi2 ~ dirichlet(alpha2);
  for (k in 1:K)
    theta_l[k] ~ normal(mu, sigma);

  for (k in 1:K){
    for (l in 1:K){
      gamma[k+ K*(l-1)] = log(pi1[k]) + log(pi2[l]) + bernoulli_lpmf(x | theta[k][l]);
    }
  }
  target += log_sum_exp(gamma);
}

サンプリングする。

fit2 = pystan.stan("sbm2.stan", data=data_dict)
fit2

                   mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
    alpha1[0]     2.3e7   1.3e7  2.3e7 6866.1  1.7e6  1.7e7  3.5e7  7.8e7    3.0   2.35
    alpha1[1]     4.4e7   2.0e7  3.4e7  2.1e4  7.3e6  4.3e7  6.1e7  1.1e8    3.0    2.5
    alpha1[2]     3.8e8   5.0e8  7.0e8  2.4e5  1.4e6  4.0e6  4.0e8  2.2e9    2.0   2.96
    alpha1[3]     3.6e7   3.6e7  5.1e7  5.1e6  7.1e6  9.3e6  3.9e7  1.7e8    2.0   3.64
    alpha1[4]     6.1e7   6.1e7  1.1e8  5.6e5  1.1e6  8.5e6  5.7e7  3.5e8    3.0   4.63
    alpha1[5]     3.9e7   3.6e7  5.1e7 2659.3 9569.2  8.9e6  6.2e7  1.5e8    2.0   4.23
    ...
    theta[8,9]     0.95    0.07   0.13   0.51    1.0    1.0    1.0    1.0    4.0   1.66
    theta[9,9]     0.93    0.08   0.16   0.38    1.0    1.0    1.0    1.0    4.0    1.8
    lp__          -4060    6.59  13.17  -4087  -4070  -4060  -4049  -4039    4.0   1.98

が、これもだめ。

諦めて一端ここで終わり。

追記(2017年1月5日)

上記はたぶん間違っているので、以下のようなモデルを作った。しかし、これもうまくいっていない。

data{
  int N;
  int K;
  int x[N,N];
}

parameters{
  simplex[K] pi1;
  simplex[K] pi2;
  real<lower=0> a0;
  real<lower=0> b0;
  /* vector<lower=0>[K] alpha1; */
  /* vector<lower=0>[K] alpha2; */
  matrix<lower=0, upper=1>[K,K] theta;
}

transformed parameters{
  matrix<upper=0>[K,K] soft_z[N, N];
  for (i in 1:N){
    for (j in 1:N){
      for (k in 1:K){
        for (l in 1:K){
          soft_z[i,j, k,l] = log(pi1[k]) + log(pi2[l]) + bernoulli_lpmf(x[i,j] | theta[k,l]);
        }
      }
    }
  }
}

model{
  /* pi1 ~ dirichlet(alpha1); */
  /* pi2 ~ dirichlet(alpha2); */
  to_vector(theta) ~ beta(a0, b0);

  for (i in 1:N){
    for (j in 1:N){
      target += log_sum_exp(soft_z[i,j]);
    }
  }
}