pytorchで Canonical Correlation Analysis (正準相関分析)の実装

はじめに

  • pytorchの練習も兼ねて,Canonical Correlation Analysis (正準相関分析)をpytorchを使って実装する.
  • 本当は分散共分散行列からなる行列の一般化固有値問題を解くが,今回は勾配法で解を求める.
  • pytorchのプログラムが間違っていないことを確認するためにscikit-learnでもやる.

Canonical Correlation Analysis (正準相関分析)

こちらの資料を参考にプログラムを書く.

主成分分析が多次元の値に対して,分散が大きい方向に射影するアルゴリズムなのに対して,正準相関分析では2つの多次元変数を射影先で相関が大きくなるように射影するアルゴリズムである.

多次元のデータ x\in\mathbb{R}^{T\times K_x}y\in\mathbb{R}^{T\times K_y} 間の正準相関を考える. ここでTはデータ数(系列データのときはデータ長),K_x, K_yx_i,y_iの次元を表す.

x,yの射影ベクトルをそれぞれ a\in\mathbb{R}^{K_x} , b\in\mathbb{R}^{K_y} とする. x,yの平均がそれぞれ0のとき,射影後の分散 S_{xx},S_{yy} と共分散 S_{xy} は, S_{xx} = (a^\top x) (a^\top x)^\top=a^\top V_{xx}a , S_{yy} = (b^\top y) (b^\top y)^\top=b^\top V_{yy}bS_{xy} = (a^\top x) (b^\top y)^\top=a^\top V_{xy}b

よって,射影後の相関\rho\rho = \frac{a^\top V_{xy}b}{\sqrt{a^\top V_{xx}a}\sqrt{b^\top V_{yy}b}}になる. ここで,aとかbを単純にk倍しても相関は変化しない. 最適化のときに解が複数あるのは面倒なので,制約条件a^\top V_{xx}a = b^\top V_{yy}b = 1を加える.

以上から以下を解くことで正準相関が求まる. \max_{a,b} a^\top V_{xy}b \ \mathrm{s.t.}\ a^\top V_{xx}a = b^\top V_{yy}b = 1

これはラグランジュの未定乗数法で解くことができて,一般化固有値問題になるが,今回はこの式のまま勾配法で解く.

生成データでのテスト

以下のように適当に生成したデータを使って,正準相関分析のプログラムをテストする.

以下では x,y,\mathrm{dummy}\in\mathbb{R}^{200 \times 2}の3つの信号を使う. x,yはもともと同じ信号に大きめのノイズが乗ったもの.dummyは適当なノイズだけのもの. x,y間では正準相関の値が大きくなり,x,dummy間では正準相関の値が小さくなってほしい.

import pylab as plt
import numpy as np
import torch 
from torch.autograd import Variable
import logging
import coloredlogs
from tqdm import tqdm

logger = logging.getLogger(__name__)
coloredlogs.install(
    #fmt="%(asctime)s %(levelname)s %(message)s"
    fmt="%(levelname)s %(message)s"
)

%matplotlib inline

logger.setLevel("DEBUG")
if logger.level >= logging.WARNING:
    tqdm = lambda x:x
u1 = np.r_[np.zeros(40), np.ones(30), np.ones(50)*2, np.ones(30), np.zeros(20), np.ones(30) ]
u2 = np.r_[np.ones(25), np.zeros(25), np.ones(25)*2, np.ones(55), np.zeros(35), np.ones(35) ]
u = np.c_[u1, u2]
plt.plot(u)
plt.show()

f:id:ksknw:20171213230955p:plain

sx = np.random.random(len(u)) * 10
sy = np.random.random(len(u)) * 10
sd = np.random.random(len(u)) * 10
x = np.c_[
    u1 + sx + np.random.random(len(u)),
    u2 - sx + np.random.random(len(u)),   
]
plt.plot(x)
plt.show()

f:id:ksknw:20171213231026p:plain

y = np.c_[
    u1 + sy + np.random.random(len(u)),
    u2 - sy + np.random.random(len(u)),   
]
plt.plot(y)
plt.show()

f:id:ksknw:20171213231039p:plain

dummy = np.c_[
    sd + np.random.random(len(u)),
    -sd + np.random.random(len(u)),   
]
plt.plot(dummy)
plt.show()

f:id:ksknw:20171213231050p:plain

以上のように,見た目では全部ぐちゃぐちゃなノイズにしか見えない3つの信号を生成した. 見た目では全部バラバラだが,x,yに関しては"いい感じの方向"に写像できればノイズを除去できる. x,yそれぞれだけでは"いい感じの方向"がわからないのでノイズと信号を区別できないが,両方を使って相関が大きくなる方向をいい感じの方向と定義することで, 相関する信号を抽出できる.

Scikit-Learnによる実装

sciikit-learnには既に実装されているので,以下のように簡単に正準相関と射影後の変数を求めることができる.

from sklearn.cross_decomposition import CCA
cca = CCA(n_components=1)

x,y間の正準相関

t_x, t_y = cca.fit_transform(x,y)
plt.scatter(t_x, t_y)
print(np.corrcoef(t_x.T, t_y.T )[0,1])
plt.show()

正準相関の値

0.906462806701

射影した後のデータの分布

f:id:ksknw:20171213231110p:plain

x,\mathrm{dummy}間の正準相関

t_x, t_d = cca.fit_transform(x,dummy)
plt.scatter(t_x, t_d)
print(np.corrcoef(t_x.T, t_y.T )[0,1])
plt.show()

正準相関の値

0.141422242479

射影した後のデータの分布

f:id:ksknw:20171213231121p:plain

これと同じ結果を得ることを目標にpytorchで実装する.

PyTorchによる実装

pytorchにはtorch.eigという固有値を求める関数が用意されているので,これを頑張って使えば解けそうではある. が,諸々の事情からペナルティ関数法で制約条件を加えた誤差関数を設定し,autogradで微分を求めて勾配法で解く.

元の最適化問題はこれだった. \max_{a,b} a^\top V_{xy}b \ \mathrm{s.t.}\ a^\top V_{xx}a = b^\top V_{yy}b = 1

制約条件をペナルティ関数にして,以下のような誤差関数を最適化する. \mathrm{argmin}_{a,b} \mathcal{L} = -a^\top V_{xy}b + \gamma \left| 1 - a^\top V_{xx}a \right| + \gamma \left| 1 - b^\top V_{yy}b \right|

\gammaは適当に0.5にした.小さい値だと制約があまり満たされず,逆に大きい値とかにするとうまく収束しないことが多い.この辺がよくわからないので,「機械学習のための連続最適化」という本を買ったけどまだ読めていない.

optimizerは何でもいいけど,単純にSGDにした.pytorch等で実装してると,なんとなくでoptimizerを選べるのでいいなと思う.

x -= np.mean(x, axis=0)
y -= np.mean(y, axis=0)
dummy -= np.mean(dummy, axis=0)
N = len(x)

def CCA(x,y,
        nb_steps = 100000,
        lr = 0.001,
        γ = 0.5):
    vxx = np.dot(x.transpose(), x)/N
    vyy = np.dot(y.transpose(), y)/N
    vxy = np.dot(x.transpose(), y)/N
    
    X = Variable(torch.from_numpy(x)).float()
    Y = Variable(torch.from_numpy(y)).float()
    
    a = Variable(torch.randn(2,1), requires_grad=True)
    b = Variable(torch.randn(2,1), requires_grad=True)
    
    Vxx = Variable(torch.from_numpy(vxx)).float()
    Vyy = Variable(torch.from_numpy(vyy)).float()
    Vxy = Variable(torch.from_numpy(vxy)).float()
    
    hist = []
    from torch import optim
    optimizer = optim.SGD([a, b], lr = 0.0001)

    for i in tqdm(range(nb_steps)):
        Loss = - torch.mm(torch.mm(torch.t(a), Vxy), b)  \
             + γ * torch.abs(1 - torch.mm(torch.mm(torch.t(a),Vxx),a))  \
             + γ * torch.abs(1 - torch.mm(torch.mm(torch.t(b),Vyy),b)) 
        hist.append(Loss.data.numpy()[0,0])
        optimizer.zero_grad()
        Loss.backward()
        optimizer.step()
        
    x_t = np.dot(a.data.numpy().transpose(), x.transpose())
    y_t = np.dot(b.data.numpy().transpose(), y.transpose())
    cca_score = torch.mm(torch.mm(torch.t(a), Vxy), b)
    
    logger.info("a Vxx a: %s"%torch.mm(torch.mm(torch.t(b),Vyy),b).data.numpy())
    logger.info("b Vyy b: %s"%torch.mm(torch.mm(torch.t(a),Vxx),a).data.numpy())
    return hist, x_t, y_t, cca_score

x,y間の正準相関

hist, x_t, y_t, cca_score = CCA(x,y)
plt.plot(hist)
plt.show()
100%|██████████| 100000/100000 [00:19<00:00, 5232.01it/s]
INFO a Vxx a: [[ 1.00017715]]
INFO b Vyy b: [[ 1.00021267]]

lossの履歴

f:id:ksknw:20171213231145p:plain

plt.scatter(x_t, y_t)
plt.show()
print(cca_score)

正準相関の値

Variable containing:
 0.9066
[torch.FloatTensor of size 1x1]

射影した後のデータの分布

f:id:ksknw:20171213231156p:plain

hist, x_t, y_t, cca_score = CCA(x,dummy)
plt.plot(hist)
plt.show()
100%|██████████| 100000/100000 [00:19<00:00, 5199.69it/s]
INFO a Vxx a: [[ 1.00273335]]
INFO b Vyy b: [[ 0.99813366]]

lossの履歴

f:id:ksknw:20171213231211p:plain

x,\mathrm{dummy}間の正準相関

plt.scatter(x_t, y_t)
plt.show()
cca_score

正準相関の値

Variable containing:
1.00000e-02 *
  9.6734
[torch.FloatTensor of size 1x1]

射影した後のデータの分布

f:id:ksknw:20171213231220p:plain

無事にpytorchによる実装で,ある程度正しそうな結果を得ることができた.

おわりに

pytorchで勾配法を使ったCCAの実装をやった. 制約付き最適化を勾配法で解く良いやり方がいまいちわかってないので,勉強したい.

参考

EmacsでブラウザのテキストボックスにMarkdownを書く (校正もする)

概要

Atomic ChromeEmacsからブラウザのテキストボタンを操作しつつ,textlintをflycheckからよんで文章を校正する. flymdでMarkdownのリアルタイムプレビューもする.

こんな感じになる. f:id:ksknw:20171104122228p:plain

はじめに

ブラウザのテキストボックスに文字を打ち込むのが面倒. ブラウザで操作するテキストボックスには,例えば,はてなブログGithubのissueなどがある. これらの操作が面倒な原因として

などがある.これらの問題を解決する方法として以下をやる.

  • Atomic ChromeEmacsからchromeのテキストボックスを操作する
  • textlintでEmacsから文法チェックする
  • flymdでEmacsからMarkdownのリアルタイムプレビューする

(Emacs最高)

Atomic ChromeEmacsからchromeのテキストボックスを操作する

chromeにはAtomic chromeというプラグインがあり,これを使うことでAtomからテキストボックスを操作できる. このプラグインEmacsから使用できるパッケージが公開されている

.emacs.d/Cask に以下を追記してインストール.

(depends-on "atomic-chrome")

また,設定ファイルのどこかに以下を書いておくと,chrome内のプラグインのボタンを押すと勝手にEmacsにバッファが開くようになる. デフォルトでは"C-x C-s"をしても保存しないようになっているが,ファイルもほしいことが多いので,これは保存のままにしておく.

(require 'atomic-chrome)
(atomic-chrome-start-server)
(bind-key* "C-x C-s" 'save-buffer) 

Emacsで書けばキーバインドだけでなく,表をOrgTblで書いたりもできるので良い. Emacsで書いたものをはてなのテキストボックスに貼り付けてもいいけど,画像のアップロードとかが面倒なので,少しだけ楽. はてなもそうだけど,sharelatexとかもEmacsから触れるので色々良い.

textlintでEmacsから文法チェックする

textlintMarkdownなどのテキストファイルの校正を行うアプリケーション. 本体にルールを追加していくことで,自分の用途にあった校正ができる.

作者の方が日本人であることもあって(?)日本語用のルールも充実していて, 例えば,ですますが混在しないようにするとか1文中に同じ助詞が複数でないようにするとか色々ある.

ルールを個別にインストールしていくこともできるが,上位の概念としてプリセットというのがあり,例えば日本語技術文書向けルールのプリセットなどがあり,かつその中に含まれるルールを個別にオフしたりもできるので,これを使うのが良さそう.

以下のようにインストールする.

$ npm i -g textlint textlint-rule-preset-japanese textlint-rule-preset-jtf-style textlint-rule-preset-ja-technical-writing textlint-spellcheck-tech-word textlint-rule-prh

~/.textlintrc に

{
  "rules": {
    "preset-japanese": {
      "sentence-length": false
    },
    "spellcheck-tech-word": true,
    "preset-ja-technical-writing": {
      "ja-no-mixed-period": {
        "periodMark": ".",
      },
      "sentence-length": false
    },
    "prh": {
     "rulePaths" :["~/prh.yml"]
    }
  }
}

とか書いておくと色々動く. 普通に使うときは以下のような感じで使える.

$ textlint hoge.md

こんな感じになる. f:id:ksknw:20171104121530p:plain

これをEmacsから使う方法としてflycheckのルールにする方法がある. ルールの設定はEmacsからでも↑のやつを勝手によんでくれるので,以下のような感じにする.

(require 'flycheck)
(flycheck-define-checker textlint
  "A linter for prose."
  :command ("textlint" "--format" "unix" source-inplace)
 :error-patterns
  ((warning line-start (file-name) ":" line ":" column ": "
            (id (one-or-more (not (any " "))))
            (message (one-or-more not-newline)
                     (zero-or-more "\n" (any " ") (one-or-more not-newline)))
            line-end))
  :modes (text-mode markdown-mode))

(add-to-list 'flycheck-checkers 'textlint)
(add-hook 'markdown-mode-hook 'flycheck-mode)

こんな感じで校正してくれる.

f:id:ksknw:20171104112919p:plain

flymdでEmacsからMarkdownのリアルタイムプレビューする

もともとMarkdown書くときは使っているのでFlymdを使っている. (できればEmacs内でなんらかのcssが反映されたプレビューを見たい).

以上全部やると,最初のような感じになる.

既知の問題点

  • 一旦ファイルを保存しないとtextlintが動かない.
  • flymdはfirefoxで動くので,chromefirefoxEmacsを開いていて頭が悪い
  • 特に表記ゆれについて,自分でやりがちなものを登録していくともっといい感じになると思う.

参考

CycleGANで普通の木を満開の桜にする

はじめに

定期的に生成系のタスクで遊びたくなる. 今回はCycle GANを使って、普通の木を満開の桜に変換してみることにした。

Cycle GAN

論文はこれ. 中身についてはたくさん解説記事があるので、そちらを参考。

Cycle GANでは2つのドメインの間の写像を学習する。 普通のGANとは異なり(乱数ではなく)、片方のドメインの画像を入力して、もう1方のドメインの画像を出力するGeneratorがある。 また、普通のGANの誤差関数(discriminatorの判別)だけでなく、1度片方のGeneratorに通した後、もう一方のGeneratorに通して再構成誤差を計算し、誤差関数に加える。

pix2pixなどでは対になる画像を用意しないと学習ができないが、CycleGANではそういうのがいらないという利点がある。

実験

ImageNetから桜の画像3000枚と普通の木の画像2500枚をダウンロードした. 画像をざっと見た感じ,桜は木全体だけでなく花だけアップの写真が多くて少し気になった. また,桜と富士山が写っている画像も結構多い.

著者の実装githubで公開されているので,それにダウンロードしたデータを食わせる. gtx780tiで2日ぐらい学習させた. 古いGPUGPUメモリが全然足りないので,オプションで色々設定した。(そろそろ買いたい)

$ python train.py --dataroot ./datasets/tree2sakura --name tree2sakura --model cycle_gan --no_dropout --pool_size 50 --loadSize 128 --fineSize 128

結果

1エポック後

元画像 生成画像 再構成画像
f:id:ksknw:20170827212704p:plain f:id:ksknw:20170827212731p:plain f:id:ksknw:20170827212752p:plain
f:id:ksknw:20170827212712p:plain f:id:ksknw:20170827212723p:plain f:id:ksknw:20170827212800p:plain

50エポック後

元画像 生成画像 再構成画像
f:id:ksknw:20170827212240p:plain f:id:ksknw:20170827212326p:plain f:id:ksknw:20170827212256p:plain
f:id:ksknw:20170827212233p:plain f:id:ksknw:20170827212315p:plain f:id:ksknw:20170827212248p:plain

100エポック後

元画像 生成画像 再構成画像
f:id:ksknw:20170827202729p:plain f:id:ksknw:20170827202725p:plain f:id:ksknw:20170827212050p:plain
f:id:ksknw:20170827202852p:plain f:id:ksknw:20170827202846p:plain f:id:ksknw:20170827212103p:plain

100エポックのネットワークで普通の木を桜に変換したもの

白と黒を変換しがちなのが少し気になるけど,なかなかうまくできたと思う.

良さそうなもの

 元画像 生成画像
f:id:ksknw:20170828195918p:plain f:id:ksknw:20170828195929p:plain
f:id:ksknw:20170828200024p:plain f:id:ksknw:20170828200034p:plain

よくなさそうなもの

 元画像 生成画像
f:id:ksknw:20170828200231p:plain f:id:ksknw:20170828200238p:plain
f:id:ksknw:20170828200318p:plain f:id:ksknw:20170828200326p:plain

空が白っぽいとその部分を黒で塗りつぶす傾向がある。なぜかはわからない。

おわりに

既存のコードを動かしただけだけど、思ったよりもうまく動いた。 が、空が異常に暗くなるなど、できてないやつも多い。体感では半々ぐらい。

せっかく対応する画像がなくても学習できるので、こちらでやられている犬猫変換のように、それなりに離れたドメイン間でやると面白そう。

生成系の研究はレッドオーシャンな感じで自分でやるのは大変そうだけど, 公開されているもので遊ぶのは楽しいので、ちょいちょいやっていきたい.

参考

pythonでNUTSの実装

はじめに

前回 ハミルトニアンモンテカルロ法の実装をやった.

今回は No U-Turn Sampler (NUTS)の実装をやる. 論文を参考にした.

コードはここにもある

github.com

NUTS

ハミルトニアンモンテカルロ (HMC)はパラメータの勾配を利用して, 効率的にMCMCサンプリングを行うことができる手法だった.

HMCの問題点は2つ.

  • 更新ステップ数 { \displaystyle L}を適切に指定しなければいけない.
  • 更新の大きさ { \displaystyle \epsilon} を適切に指定しなければいけない.

NUTSは更新ステップ数{ \displaystyle L}を自動的に決定する手法である. 論文内では{ \displaystyle \epsilon}はdual-averaging (Nesterov 2009)を用いて決定するが,今回は決め打ちにする.

更新ステップ数L

ハミルトニアンモンテカルロでは,正規分布によって発生させた運動量を与えて, { \displaystyle L}ステップの間,点を動かす. 予め決められた{ \displaystyle L}ステップの間,点を動かすので,例えば谷にハマった時などガタガタして無駄な計算をしてしまう.

NUTSでは「Uターンしたら点を動かすのをやめる」という規則でこの無駄な計算をなくす. ただし,単純にUターンしたときに中断したら詳細釣り合い条件を満たさなくなるので,少し工夫する.

具体的には,だいたい

  1. 運動量{ \displaystyle r_0}をランダムに決める.
  2. 時間の向き{ \displaystyle v}を{-1, 1}からランダムに選ぶ.
  3. 選んだ時間の向きの方向に{ \displaystyle 2^j}回,点を移動させる(点の軌跡は全て記憶しておく)
  4. どこかでUターンしていないかを確認する
  5. Uターンしていたら6へ,それ以外は{ \displaystyle j+1}して2へ戻る
  6. これまでに計算した軌跡からランダムに1点選んで,サンプリング結果に加えて,1に戻る

という感じ.

実装

以下のようなデータの平均と分散をサンプリングする.

import numpy as np
from numpy import exp
from copy import deepcopy
import pylab as plt
import seaborn as sns
from scipy.stats import norm, gamma
from tqdm import tqdm
from numpy import random
from matplotlib.animation import FuncAnimation

sns.set_style("white")

true_μ = 3
true_σ = 1
nb_data = 1000

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

print(x.mean(), x.std())
sns.kdeplot(x)
plt.show()
2.99882646635 1.02180321215

f:id:ksknw:20170806182124p:plain

NUTSのプログラムは以下. 論文のNaive-NUTSの実装をする. Leapfrogで点を移動させるなどはハミルトニアンモンテカルロと同じ.

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

def H(θₜ, p):
    return -norm_lpdf(θₜ[0], θₜ[1]) + 0.5 * np.dot(p, p)

def Leapfrog(x, θ, r, ε):
    θ_d = deepcopy(θ)
    r_d = deepcopy(r)
    r_d -= 0.5 * ε * log_dh(θ_d[0], θ_d[1])
    θ_d[0] = θ_d[0] + ε * r_d[0]
    θ_d[1] = θ_d[1] + ε * r_d[1]
    r_d -= 0.5 * ε * log_dh(θ_d[0], θ_d[1])
    return θ_d, r_d
norm_lpdf = lambda μ, σ: np.sum(norm.logpdf(x, μ, σ))
gamma_lpdf = lambda a: np.sum(gamma.logpdf(x, a))

Δ_max = 1000
ε = 0.05
L = norm_lpdf
M = 100

θ0 = np.array([random.randn(), random.gamma(1)])
list_θₘ = [θ0]

{ \displaystyle 2^j} 回の移動は再帰で実装される.

def BuildTree(θ, r, u, v, j, ε):
    if j == 0:
        θd, rd = Leapfrog(x, θ, r, v * ε)
        if np.log(u) <= (L(*θd) - 0.5 * np.dot(rd, rd)):
            Cd_ = [[θd, rd]]
        else:
            Cd_ = []
        sd = int(np.log(u) < (Δ_max + L(*θd) - 0.5 * np.dot(rd, rd)))
        return θd, rd, θd, rd, Cd_, sd
    else:
        θ_minus, r_minus, θ_plus, r_plus, Cd_, sd = BuildTree(θ, r, u, v, j - 1, ε)
        if v == -1:
            θ_minus, r_minus, _, _, Cdd_, sdd = BuildTree(θ_minus, r_minus, u, v, j - 1, ε)
        else:
            _, _, θ_plus, r_plus, Cdd_, sdd = BuildTree(θ_plus, r_plus, u, v, j - 1, ε)
        sd = sdd * sd * int((np.dot(θ_plus - θ_minus, r_minus) >= 0) and (np.dot(θ_plus - θ_minus, r_plus) >= 0))
        Cd_.extend(Cdd_)

        return θ_minus, r_minus, θ_plus, r_plus, Cd_, sd
hist_L = []
hist_C = []
for m in tqdm(range(M)):
    r0 = random.randn(2)
    u = random.uniform(0, exp(L(*list_θₘ[-1]) - 0.5 * np.dot(r0, r0)))

    θ_minus = deepcopy(list_θₘ[-1])
    θ_plus = deepcopy(list_θₘ[-1])
    r_minus = deepcopy(r0)
    r_plus = deepcopy(r0)
    j = 0
    C = [[deepcopy(list_θₘ[-1]), deepcopy(r0)]]
    s = 1

    while s == 1:
        v = random.choice([-1, 1])
        if v == -1:
            θ_minus, r_minus, _, _, Cd, sd = BuildTree(θ_minus, r_minus, u, v, j, ε)
        else:
            _, _, θ_plus, r_plus, Cd, sd = BuildTree(θ_plus, r_plus, u, v, j, ε)

        if sd == 1:
            C.extend(Cd)
        s = sd * int((np.dot(θ_plus - θ_minus, r_minus) >= 0) and (np.dot(θ_plus - θ_minus, r_plus) >= 0))
        j += 1

    index = random.choice(list(range(len(C))))
    list_θₘ.append(C[index][0])

    hist_L.append(L(C[index][0][0], C[index][0][1]))
    hist_C.append(C)
100%|██████████| 100/100 [00:00<00:00, 760.56it/s]
def plot(list_θₘ, hist_C):
    fig = plt.figure()
    ax = fig.add_subplot(111)

    def update(i):
        fig.canvas.draw()
        ax.cla()
        j = i // 3
        if (i % 3) == 0:
            ax.scatter(np.array(hist_C[j])[:, 0, 0], np.array(hist_C[j])[:, 0, 1], linewidth=0, marker=".")
            ax.plot(list_θₘ[:(j), 0], list_θₘ[:(j), 1], c="gray", linewidth=0.3, alpha=0.4)
            ax.scatter(list_θₘ[:(j), 0], list_θₘ[:(j), 1], c="g", linewidth=0, marker=".")

        elif (i % 3) == 1:
            ax.scatter(np.array(hist_C[j])[:, 0, 0], np.array(hist_C[j])[:, 0, 1], linewidth=0, marker=".")
            ax.plot(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="gray", linewidth=0.3, alpha=0.4)
            ax.scatter(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="g", linewidth=0, marker=".")
        else:
            ax.scatter(np.array(hist_C[j])[:, 0, 0], np.array(hist_C[j])[:, 0, 1], linewidth=0, marker=".", c="w")
            ax.plot(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="gray", linewidth=0.3, alpha=0.4)
            ax.scatter(list_θₘ[:(j + 1), 0], list_θₘ[:(j + 1), 1], c="g", linewidth=0, marker=".")

        plot_lim = 30
        if(j) > plot_lim:
            temp_list = np.array([[np.min(np.array(C)[:, 0], axis=0),
                                   np.max(np.array(C)[:, 0], axis=0)] for C in hist_C[j - plot_lim: j]])

            temp_xlim = [np.min(temp_list[:, :, 0]), np.max(temp_list[:, :, 0])]
            xlim_range = temp_xlim[1] - temp_xlim[0]
            temp_ylim = [np.min(temp_list[:, :, 1]), np.max(temp_list[:, :, 1])]
            ylim_range = temp_ylim[1] - temp_ylim[0]
            plt.xlim(temp_xlim[0] - xlim_range * 0.1, temp_xlim[1] + xlim_range * 0.1)
            plt.ylim(temp_ylim[0] - ylim_range * 0.1, temp_ylim[1] + ylim_range * 0.1)

    ani = FuncAnimation(fig, update, frames=len(hist_C) * 3 - 2)
    ani.save("temp.gif", writer="imagemagick", fps=3)
    
plot(list_θₘ, hist_C)

たぶんできていると思う きれいに書き直すつもりだったけど、面倒だったからやめた

「2017年3月14日をもって動画のアップロード機能は終了いたしました」 じゃないんだよなぁ

参考

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

はじめに

今までなんとなく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はまとめられそうだけど(というか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

おわりに

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

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

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

参考