ハミルトニアンモンテカルロ法の実装

はじめに

今までなんとなくStanなどを使ってMCMCをやっていた。 ギブスサンプリングぐらいなら昔勉強したけど、ハミルトニアンモンテカルロや、ましてやNUTSなどは何をやっているのかあまり理解していなかった。 基礎からのベイズ統計学という本を読んで、ハミルトニアンモンテカルロまではなんとなくわかったので、プログラムを書いて理解する。

NUTSはまだわかってないので、

を実装する。

理論は他の方のブログや本を参考にしてほしい。 自分の理解で適当に書いているので、間違ってたらごめんなさい。

データ

なんでもいいけど、平均3標準偏差1の1次元の正規分布に従うデータから、その平均と分散を推定することにする。

import numpy as np
from numpy import random
from copy import deepcopy
import pylab as plt
import seaborn as sns
from matplotlib.animation import FuncAnimation
from scipy.stats import norm, gamma
from tqdm import tqdm
sns.set_style("white")

true_μ = 3
true_σ = 1
nb_data = 1000

data = np.random.normal(true_μ, true_σ, nb_data)

print(data.mean(), data.std())

sns.kdeplot(data)
plt.show()
    3.02683958884 0.998470831061

    /home/kei/anaconda3/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
      y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j

f:id:ksknw:20170709222104p:plain

メトロポリスヘイスティングス

メトロポリスヘイスティングス法(MH法)では、

  1. { \displaystyle
a \sim q(a|\theta_t)
} からパラメータの候補aを生成する
  2. { \displaystyle r = \frac{q(\theta_t | a)f(a)}{q(a|\theta_t)f(\theta_t)} } を計算する。
  3. 確率 { \displaystyle \min(1,r) }でaを受容し、{ \displaystyle \theta_{t+1}=a }とする。

を繰り返してパラメータを更新していく。 よって、現在のパラメータから新しいパラメータを生成するときの確率{ \displaystyle q(a|\theta_t)}、および、尤度{ \displaystyle f}を計算するときの分布の形状(提案分布)を決めると、パラメータを推定することができる。

今回は提案分布を正規分布として、その平均と標準偏差を推定する。

独立メトロポリスヘイスティングス

まずは{ \displaystyle q(a|\theta_t) = q(a)}、つまり次のパラメータを前のパラメータに依存せず独立にサンプリングする。 今回は平均と標準偏差の推定なので、それぞれ正規分布およびガンマ分布からパラメータをサンプリングする。

sampling_μ_mean = 0
sampling_μ_std = 10
sampling_σ_a = 2

norm_lpdf = lambda x, μ, σ: np.sum(norm.logpdf(x, μ, σ))
gamma_lpdf = lambda x, a: np.sum(gamma.logpdf(x, a))

logf = lambda θₜ: norm_lpdf(data, θₜ["μ"], θₜ["σ"])
logq = lambda θₜ: norm_lpdf(θₜ["μ"], sampling_μ_mean, sampling_μ_std) \
    + gamma_lpdf(θₜ["σ"], sampling_σ_a)
logr = lambda a, θₜ: logq(θₜ) + logf(a) - logq(a) - logf(θₜ)

nb_samples = 1000
θₜ = {}
a = {}  # 候補点

hist = []

θₜ["μ"] = random.normal(sampling_μ_mean, sampling_μ_std)
θₜ["σ"] = random.gamma(sampling_σ_a)

for t in tqdm(range(nb_samples)):
    a["μ"] = random.normal(sampling_μ_mean, sampling_μ_std)
    a["σ"] = random.gamma(sampling_σ_a)

    if logr(a, θₜ) > 0:
        θₜ = deepcopy(a)
        hist.append([a["μ"], a["σ"], 0])
    else:
        if np.exp(logr(a, θₜ)) > random.uniform(0, 1):
            θₜ = deepcopy(a)
            hist.append([a["μ"], a["σ"], 1])
        else:
            hist.append([a["μ"], a["σ"], 2])
hist = np.array(hist)

print(θₜ)
    100%|██████████| 1000/1000 [00:03<00:00, 323.83it/s]

    {'σ': 1.1803787969064548, 'μ': 2.801370701839039}

サンプリングがどのように進んだかをプロットしてみる。

def update(i, hist):
    if hist[i, 2] == 0:
        plt.scatter(hist[i, 0], hist[i, 1], marker=".", c="b", alpha=0.8)
    elif hist[i, 2] == 1:
        plt.scatter(hist[i, 0], hist[i, 1], marker=".", c="g", alpha=0.8)
    elif hist[i, 2] == 2:
        plt.scatter(hist[i, 0], hist[i, 1], marker="x", c="r", alpha=0.4)

fig = plt.figure()
plt.xlabel(r"$\mu$")
plt.ylabel(r"$\sigma$")
ani = FuncAnimation(fig, update, fargs=(hist,), frames=nb_samples)
ani.save("independent.gif", writer="imagemagick")

f:id:ksknw:20170709221856g:plain

青い点が受容されたパラメータ、赤いバツが棄却された候補を表している。 殆どの候補が棄却されて、むちゃくちゃ効率が悪いことがわかる。

ランダムウォークMH法

{ \displaystyle a = \theta_t + e}として、新しいパラメータ候補{ \displaystyle a}を生成すると、これまで見つかった良さそうなパラメータの周辺を探索することになるので、受容率が上がって効率のいいサンプリングが行える。

今回は{ \displaystyle e}が平均0、標準偏差0.05の正規分布に従うとしてサンプリングを行った。 このとき、{ \displaystyle q(a|\theta_t) = \mathcal{N}(\theta_t|0.05)} になる。 { \displaystyle q(\theta_t|a)}も同様。

sampling_mean = 0
sampling_std = 0.05

logf = lambda θₜ: norm_lpdf(data, θₜ["μ"], θₜ["σ"])
logq = lambda a, θₜ: norm_lpdf(θₜ["μ"], sampling_mean+a["μ"], sampling_std) \
    + norm_lpdf(θₜ["σ"], sampling_mean+a["σ"], sampling_std)
logr = lambda a, θₜ: logq(θₜ,a) + logf(a) - logq(θₜ,a) - logf(θₜ)

nb_samples = 1000
θₜ = {}
a = {}  # 候補点

hist = []
hist_2sigma = []
nb_i = 40
sample_2sigma = [np.array([2*sampling_std*np.sin(i*2*np.pi/nb_i), 
                           2*sampling_std*np.cos(i*2*np.pi/nb_i)]) for i in range(nb_i)]

θₜ["μ"] = random.normal(1, 1)
θₜ["σ"] = random.gamma(3)

for t in tqdm(range(nb_samples)):
    a["μ"] = θₜ["μ"] + random.normal(sampling_mean, sampling_std)
    a["σ"] = θₜ["σ"] + random.normal(sampling_mean, sampling_std)

    if logr(a, θₜ) > 0:
        θₜ = deepcopy(a)
        hist.append([a["μ"], a["σ"], 0])
    else:
        if np.exp(logr(a, θₜ)) > random.uniform(0, 1):
            θₜ = deepcopy(a)
            hist.append([a["μ"], a["σ"], 1])
        else:
            hist.append([a["μ"], a["σ"], 2])
            
    hist_2sigma.append([np.array([θₜ["μ"],θₜ["σ"]]) + sample_2sigma_i 
                        for sample_2sigma_i in sample_2sigma])

hist = np.array(hist)
hist_2sigma = np.array(hist_2sigma)
print(θₜ)
    100%|██████████| 1000/1000 [00:01<00:00, 592.47it/s]

    {'σ': 0.9944398277566455, 'μ': 3.0750991483862506}
def plot(hist, hist_2sigmam, filename, T=1):
    lim_range = 0.2
    fig = plt.figure()
    ax = fig.add_subplot(111)
    fill, = ax.fill(hist_2sigma[0, :, 0], hist_2sigma[0, :, 1], alpha=0.2, c="gray")
    plt.xlabel(r"$\mu$")
    plt.ylabel(r"$\sigma$")

    def update(i, hist, hist_2sigma, line):
        fill.set_xy(hist_2sigma[i])
        fig.canvas.draw()
        
        if hist[i, 2] == 0:
            ax.scatter(hist[i, 0], hist[i, 1], marker=".", c="b", alpha=0.8, linewidths=0)
        elif hist[i, 2] == 1:
            ax.scatter(hist[i, 0], hist[i, 1], marker=".", c="g", alpha=0.8, linewidths=0)
        else:  # hist[i, 2] == 2:
            ax.scatter(hist[i, 0], hist[i, 1], marker="x", c="r", alpha=0.4)
        
        plot_lim = 10
        if i > plot_lim:
            plt.xlim(min(min(hist[i - plot_lim:i, 0]), 3-lim_range), 
                     max(max(hist[i - plot_lim:i, 0]), 3+lim_range))
            plt.ylim(min(min(hist[i - plot_lim:i, 1]), 1-lim_range), 
                     max(max(hist[i - plot_lim:i, 1]), 1+lim_range))

    ani = FuncAnimation(fig, update, fargs=(hist, hist_2sigma, line), frames=nb_samples//T)
    ani.save(filename, writer="imagemagick")
plot(hist, hist_2sigma, "random_walk.gif", 3)

f:id:ksknw:20170709220916g:plain

青い点および赤いバツは上と同じ、確率rで受容されたパラメータ。 96%ぐらいのパラメータの候補が灰色の領域から選ばれる。 先ほどと比べるとだいぶ多くのパラメータが受容されている。

プロットの範囲がいまいちなのは、どうすればいい感じになるのかよくわからなかったからです…

ハミルトニアンモンテカルロ法

ランダムウォークMH法は{ \displaystyle e}を適切に設定しないと特にパラメータの次元が高くなるとうまくいかない(らしい)。 パラメータによる勾配が利用できるのであれば、それを使ってパラメータを更新すると良さそうである。

ハミルトニアンモンテカルロ法では慣性付きの勾配法っぽい感じでパラメータを更新する。

  1. モーメンタム{ \displaystyle p}を標準正規乱数から発生させる。
  2. ポテンシャル{ \displaystyle h(\theta)}を対数事後分布として、慣性付きの勾配法を進めて、{ \displaystyle L}ステップ後の値を{ \displaystyle \theta_a}とする。
  3. { \displaystyle r = exp(H(\theta_t, p_t) - H(\theta_a, p_a))}を求める。
  4. 確率 { \displaystyle \min(1,r)}でパラメータ{ \displaystyle \theta_a}を受容する。

ここで{ \displaystyle H}ハミルトニアン{ \displaystyle H(\theta, p) = h(\theta) + \frac{1}{2} p \cdot p}

普通の勾配法と違うのは、

  • { \displaystyle L}ステップ更新するごとに、慣性項を標準正規乱数でリセットする。
  • まれに更新を棄却する
  • 粘性項がない

といったところかと思う。
(普通の確率勾配法で更新したらいかんのかと思ってしまうけど、そういう研究もあるんだろうな)

def log_dh(x, μ, σ):
    return np.array([-np.sum(x - μ) / σ**2,
                     len(x) / (2 * σ**2) - np.sum((x - μ)**2) / (2 * σ**4)])

def H(x, θₜ, p):
    return -norm_lpdf(x, θₜ["μ"], θₜ["σ"]) + 0.5 * np.dot(p, p)

nb_samples = 1000
θₜ = {}
θₐ = {}

ε = 0.01
L = 100

hist = []
hist_2sigma = []
θₜ["μ"] = random.normal(1, 1)
θₜ["σ"] = random.gamma(3)
nb_i = 40
p_2sigma = [np.array([2*np.sin(i*2*np.pi/nb_i), 
                      2*np.cos(i*2*np.pi/nb_i)]) for i in range(nb_i)]

def update(pₜ, θₜ):
    pₐ = deepcopy(pₜ)
    θₐ = deepcopy(θₜ)
    for t in range(L):
        pₐ -= ε / 2 * log_dh(data, θₐ["μ"], θₐ["σ"])
        dμ, dσ = ε * pₐ
        θₐ["μ"] += dμ
        θₐ["σ"] += dσ
        θₐ["σ"] = np.abs(θₐ["σ"])
        pₐ -= ε / 2 * log_dh(data, θₐ["μ"], θₐ["σ"])    
    return  pₐ, θₐ 

for i in tqdm(range(nb_samples)):
    pₜ = np.random.randn(2)
    pₐ, θₐ = update(pₜ, θₜ)
    if np.exp(H(data, θₜ, pₜ) - H(data, θₐ, pₐ))>1:
        θₜ = deepcopy(θₐ)
        hist.append([θₐ["μ"], θₐ["σ"], 0])
    elif np.exp(H(data, θₜ, pₜ) - H(data, θₐ, pₐ)) > random.uniform(0,1):
        θₜ = deepcopy(θₐ)
        hist.append([θₐ["μ"], θₐ["σ"], 1])
    else:
        hist.append([θₐ["μ"], θₐ["σ"], 2])
    list_θ = [update(p_2sigma_i, θₜ)[1] for p_2sigma_i in p_2sigma]
    hist_2sigma.append(np.array([[θ["μ"],θ["σ"]] for θ in list_θ]))

hist = np.array(hist)
hist_2sigma = np.array(hist_2sigma)
print(θₜ)
    100%|██████████| 1000/1000 [05:34<00:00,  2.24it/s]

    {'σ': 1.0048874933483771, 'μ': 3.0559467768730531}
plot(hist, hist_2sigma, "hamiltonian.gif", 3)

f:id:ksknw:20170709220450g:plain

だいたい今までと同じだけど、灰色の部分はpを発生させるときに96%ぐらいにしているので、勾配法で更新されたあとでこの範囲にサンプリングが収まっているかは保証されないと思う。 matplotlibのfill関数の使い方が間違っていなければ、最初のほうでサンプリングの範囲がぐちゃぐちゃになっている部分があるが、これは関数がぐちゃぐちゃだから(??)

ランダムウォークと比べても、棄却されるサンプルが圧倒的に少なく効率的にサンプリングできている。

おわりに

自分の勉強のために

を実装した。

勾配法との関係がいまいちわかっていなくて、例えばAdamとか使っても詳細釣り合い条件が満たされるのかよくわからないし、そもそも勾配法だけではいけないのか。 次はNUTSを理解したい。

参考

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に続くかも。

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

概要

関係データ学習を見ながらpythonSBMの実装をした。
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の事後確率

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

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

ここで、

はガンマ関数で、 はパラメータ

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

以上をそのままpythonのプログラムにした。 z1とz2はまとめられそうだけど、ややこしくなるのが嫌だったので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

おわりに

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

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

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

参考

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

概要

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

はじめに

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

%matplotlib inline
import numpy as np
import pylab as plt
import seaborn
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:20170326233209p: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のほうがだいぶ近そうという結果が得られた。

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

import matplotlib.gridspec as gridspec

def plot_path(m, 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)
    
    path = []
    path.append(m[-1][-1][1])
    while True:
        path.append(m[path[-1][0]][path[-1][1]][1])
        if path[-1]==(0,0):
            break
    path = np.array(path)
    ax2.plot(path[:,1], path[:,0], c="r")
    
    ax1.plot(A, range(len(A)))
    ax1.invert_xaxis()
    ax4.plot(B, c="g")
    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)
plot_path(m, A, B)

f:id:ksknw:20170326233224p:plain

f:id:ksknw:20170326233229p:plain

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

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

f:id:ksknw:20170326233240p:plain

f:id:ksknw:20170326233252p: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)
plot_path(m, A, D)

f:id:ksknw:20170326233307p:plain

f:id:ksknw:20170326233317p: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:20170326233331p:plain

f:id:ksknw:20170326233338p:plain

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

実データ(国内1週間の気温データ)

正弦波ばっかりいじっててもあれなので、気象庁のサイトから適当な地点の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

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

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

おわりに

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

参考

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

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

以下では,ベクトル{ \displaystyle
\bf{x}}と行列{ \displaystyle
A}を使う.
{ \displaystyle
{\bf x}  \in \mathbb{R}^2}
{ \displaystyle
A  \in \mathbb{R}^{2x2}}

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])と書くと列ベクトルっぽくなるが,後の処理が面倒になるのでしない.

トレース,行列式逆行列 { \displaystyle \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. ]])

ベクトルの内積 { \displaystyle 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]])

行列とベクトルの積 { \displaystyle 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]])

行列の積 { \displaystyle 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]])

要素ごとの積 (アダマール積){\displaystyle  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]])

ベクトルから行列 (テンソル積 ?) { \displaystyle 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で書くことになるので,慣れてないとバグを生む可能性が高い気がする. 皆さんどうしてるんですか.

関係データ学習の実装 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

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

まとめ

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

参考

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

はじめに

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

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

環境

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とか使っていきたい。

参考