読者です 読者をやめる 読者になる 読者になる

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

参考

関係データ学習の実装 ツイッターデータのスペクトルクラスタリングと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

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

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

{ \displaystyle
L=D-X}

ここで{ \displaystyle
D}は次数行列で、オブジェクトから伸びるリンクの総数。 { \displaystyle
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以上の値をとるのでこんな感じの図になる。

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

クラスタ数は適当に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アカウントしかないクラスタがいくつかできた。 正規化しない単純なグラフラプラシアンだとこういう結果になるらしい。

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

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

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

{ \displaystyle
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

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

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

{ \displaystyle
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)を使って非対称データのままクラスタリングする。

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

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

{ \displaystyle
\pi_1 \mid \alpha_1 \sim \mathrm{Dirichlet}(\alpha_1) }
{ \displaystyle
\pi_2 \mid \alpha_2 \sim \mathrm{Dirichlet}(\alpha_2)}
{ \displaystyle
z_{1,i} = k\mid \pi_1 \sim \mathrm{Discrete}(\pi_1)}
{ \displaystyle
z_{2,i} = l\mid \pi_2 \sim \mathrm{Discrete}(\pi_2)}
{ \displaystyle
\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

が、これもだめ。

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

おわりに

関係データ学習の学習のために,対称データのスペクトルクラスタリングと非対称データのSBMを実装した。 スペクトルクラスタリングの結果、なんとなくクラスタを抽出できた。 SBMをStanで実装しようとしたが、いまいち収束しなかった。

とりあえず、最後まで読んで、頑張れたら頑張ります。

参考

[http://www.kspub.co.jp/book/detail/1529212.html:title]

追記(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]);
    }
  }
}

Stanで気温変動のフィッティング

はじめに

この記事は Stan Advent Calendar 2016n日目の記事ではありません。 21日目の記事です。

話題のアヒル本(言われるまで何がアヒルなのかわからなかった)を読んだ。

statmodeling.hatenablog.com

Stanの文法を知るために買ったけれど、全体的に統計モデリング、解析について書いてあって、勉強になった。 本を読んで勉強になったことはたくさんあるけど、特にグラフを死ぬほど書かないといけないという認識を持った。 自分ではそれなりにグラフを書いているつもりだったけれど、pariplotとかviolinplotとかちゃんと書いたことがなかった。

せっかく本を読んだので、適当にデータをとってきて解析してAdvent Calendarなるものに挑戦しようと思っていたけど、もう埋まっていた。悲しいので一人で勝手にやることにする。 よく見たら21日が空いていたので、フライングということにさせてください。

今回は特に12章を見ながら時系列データをいじくることにした。

知識がないので、間違ったことをしていたら教えていただけると嬉しいです。

データ

時系列データとして天気に関するデータを使うことにした。 気象庁のサイト から天気に関するいろんなデータがダウンロードできるので利用する。

今回は愛知県の気温とか降水量とかを過去5年分ぐらいダウンロードした。 日本語のファイルがやや面倒だった。

import pandas as pd
aichi = pd.read_csv("./weather_aichi.csv", encoding='SHIFT-JIS')
aichi["年"] = aichi[u"年月日"].apply(lambda x: x.split("/")[0])
aichi["月_s"] = aichi[u"年月日"].apply(lambda x: "%02d"%int(x.split("/")[1]))
aichi["日_s"] = aichi[u"年月日"].apply(lambda x: "%02d"%int(x.split("/")[2]))
aichi["m"] = aichi[u"年月日"].apply(lambda x: int(x.split("/")[1]))
aichi["d"] = aichi[u"年月日"].apply(lambda x: int(x.split("/")[2]))
aichi["is_rain"]  = aichi[u"降水量の合計(mm)"]>0
aichi = aichi.set_index(u"年月日")
month_groups = aichi.groupby("月_s")
aichi
Unnamed: 0 平均気温(℃) 降水量の合計(mm) 月_s 日_s m d is_rain
年月日
2011/12/3 1 12.2 8.5 2011 12 03 12 3 True
2011/12/4 2 11.9 0.0 2011 12 04 12 4 False
2011/12/5 3 10.4 0.0 2011 12 05 12 5 False
2011/12/6 4 8.0 0.5 2011 12 06 12 6 True
... ... ... ... ... ... ... ... ... ...
2016/11/30 1825 8.2 0.0 2016 11 30 11 30 False
2016/12/1 1826 12.0 0.0 2016 12 01 12 1 False
2016/12/2 1827 11.0 0.0 2016 12 02 12 2 False
2016/12/3 1828 9.8 0.0 2016 12 03 12 3 False

1828 rows × 9 columns

今回は、平均気温をモデリングしてみることにした。 とりあえずグラフを書く。

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

ax = aichi[[ u"平均気温(℃)"]].plot(figsize=(20,10))

f:id:ksknw:20161211134350p:plain

当たり前だけど春夏秋冬がある。 雨のデータも後で使うのでプロットする。

ax = aichi[ u"平均気温(℃)"]["2011/12/3": "2012/12/3"].plot(figsize=(20,10))
aichi[ u"降水量の合計(mm)"]["2011/12/3": "2012/12/3"].plot( kind="bar", ax=ax)
plt.show()

f:id:ksknw:20161211134429p:plain

その他色々プロットを書く。

sns.pairplot(aichi[[u"平均気温(℃)",   u"降水量の合計(mm)", "m"]], hue="m", size=4)

f:id:ksknw:20161211134443p:plain

アヒル本では対角線の上と下で別のプロットを書いていたけど、正直手でやるのはちょっと面倒。 変数の型とか見て自動でやるやつをそのうち作りたい。

fig, ax = plt.subplots(figsize=(20,20))
sns.violinplot(data=aichi, x="m", y =u"平均気温(℃)")
fig, ax = plt.subplots(figsize=(20,20))
sns.violinplot(data=aichi, x="m", y =u"平均気温(℃)", hue="is_rain", split=True)
plt.show()

f:id:ksknw:20161211134455p:plain

f:id:ksknw:20161211134516p:plain

雨が降ると夏は気温が下がるけど、冬は逆に気温が高い。そりゃそうだな。

temperatures = aichi[u"平均気温(℃)"].get_values()
data = {"Y":temperatures,
       "N":len(aichi),
       "is_rain":aichi["is_rain"].astype(int).get_values()+1}

状態空間モデルでのフィッティング

アヒル本12章を見ながら、いくつかのモデルを使って気温の変動をモデル化する。

状態空間モデル(1階差分)

観測された気温({ \displaystyle
Y
})は 真の気温({ \displaystyle
\mu
})を平均とする正規分布に従い、 真の気温は、前日の真の気温を平均とする正規分布に従うというモデルを仮定する。 式的には以下。

{ \displaystyle
 \mu_n \sim \mathcal{N}(\mu_{n-1}, s_{\mu})
}

{ \displaystyle
 Y_n \sim \mathcal{N}(\mu_{n}, s_{Y})
}

以下のようにStanでモデルを書ける。

data{
  int N;
  vector[N] Y;
}

parameters{
  vector[N] mu;
  real<lower=0> s_mu;
  real<lower=0> s_Y;
}

model{
  mu[2:N] ~ normal(mu[1:(N-1)], s_mu);
  Y ~ normal(mu, s_Y);
}

python からStanをよんでサンプリングして、色々とプロットする。

import pystan
import numpy as np
def fit(model_filename):
    fit = pystan.stan(model_filename, data=data)
    ms = fit.extract()
    fit.plot()
    return fit, ms

def fit_plot(fit, ms, list_max_x=[-1,500,100,50], ylim=None):
    mu_mean = ms["mu"].mean(axis=0)
    mu_std  = ms["mu"].std(axis=0)

    def plot_line(max_x):
        fig = plt.figure(figsize=(10,10))
        plt.plot(mu_mean[:max_x], linewidth=0.5)
        
        plt.scatter(range(len(temperatures[:max_x])), temperatures[:max_x], marker=".",# linestyle="None", 
                    c=aichi["is_rain"].get_values()[:max_x].astype(float))
    
        plt.fill_between(range(len(mu_mean[:max_x])), (mu_mean-mu_std)[:max_x], (mu_mean+mu_std)[:max_x], alpha=0.3)
        plt.fill_between(range(len(mu_mean[:max_x])), (mu_mean-mu_std*2)[:max_x], (mu_mean+mu_std*2)[:max_x], alpha=0.2)

        plt.xlim(0, len(temperatures[:max_x]))
        if not ylim==None:
            plt.ylim(*ylim)
        plt.show()
   
    for max_x in list_max_x:
        plot_line(max_x)
fit1 = fit("./model.stan") 
fit_plot(*fit1, list_max_x=[-1,100])

f:id:ksknw:20161211140401p:plain

1000のあたりとか、見るからにサンプリングが安定していない感じがある。 サンプリング結果を見てみる。

fit1[0]
           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu[0]     12.18  6.1e-3   0.39  11.34  11.99  12.19  12.36  12.98 4000.0    1.0
mu[1]     11.84  6.0e-3   0.38  10.98  11.66  11.87  12.04  12.63 4000.0    1.0
mu[2]     10.36  6.1e-3   0.38   9.49  10.16  10.39  10.54  11.15 4000.0    1.0
mu[3]      8.21  6.5e-3   0.41   7.47   7.96   8.12   8.46   9.13 4000.0   1.05
mu[4]      9.99  6.2e-3   0.39    9.1   9.77  10.04   10.2  10.74 4000.0   1.02

mu[1825]  11.77  6.5e-3   0.41  10.85  11.53  11.87  12.04  12.48 4000.0   1.07
mu[1826]  10.97  6.1e-3   0.39  10.14  10.78  10.99  11.18  11.76 4000.0    1.0
mu[1827]   9.86  6.3e-3    0.4   9.08   9.67   9.83  10.06  10.75 4000.0    1.0
s_mu       1.72    0.04   0.07   1.59   1.67   1.72   1.78   1.85    3.0   1.71
s_Y        0.37    0.12   0.18   0.07    0.2   0.42   0.49   0.63    2.0   3.15
lp__     -622.4  903.47 1277.7  -1830  -1483  -1218 121.56 1970.3    2.0   5.28

アヒル本によると、

n_effはStanが自己相関等から判断した実効的なMCMCサンプル数である。

とあるので、今回のサンプリングではs_muやs_Yはまともにサンプリングされていないことがわかる。

また、Rhatはパラメータが収束したかどうかを表す値であり、

「chain数が3以上ですべてのパラメータでRhat < 1.1 となること」を「収束した」と見なすことにする。

とある。 ということで、このモデルではうまく収束しなかった。 サンプリング数を増やしたりしてもいいのかもしれないが、今回は別のモデルを作ることにする。

状態空間モデル(2階差分)

真の気温の変化はその前の日の真の気温の変化を平均とする正規分布に従うというモデル。 式的には以下。

{ \displaystyle
 \mu_n-\mu_{n-1} \sim \mathcal{N}(\mu_{n-1}-\mu_{n-2}, s_{\mu})
}

{ \displaystyle
 Y_n \sim \mathcal{N}(\mu_{n}, s_{Y})
}

Stanのコードは以下。

data{
  int N;
  vector[N] Y;
}

parameters{
  vector[N] mu;
  real<lower=0> s_mu;
  real<lower=0> s_Y;
}

model{
  mu[3:N] ~ normal(2 * mu[2:(N-1)] - mu[1:(N-2)], s_mu);
  Y ~ normal(mu, s_Y);
}

サンプリングする。

fit2 = fit("./model2.stan") 
fit_plot(*fit2, list_max_x=[-1,100])

f:id:ksknw:20161211140436p:plain

サンプリング結果

fit2[0]
           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu[0]     12.12    0.02   1.06  10.07  11.42  12.13  12.83  14.17 4000.0    1.0
mu[1]     11.26    0.01   0.76   9.76  10.75  11.25  11.78  12.75 4000.0    1.0
mu[2]     10.41    0.01   0.67   9.12   9.95  10.41  10.86  11.73 4000.0    1.0
mu[3]      9.64    0.01   0.68   8.28   9.19   9.65  10.09  10.97 4000.0    1.0
mu[4]      9.04    0.01   0.68    7.7   8.59   9.04   9.49  10.38 4000.0    1.0

mu[1825]  10.28    0.01   0.67   8.99   9.82  10.28  10.73  11.59 4000.0    1.0
mu[1826]   10.4    0.01   0.75   8.94    9.9  10.38   10.9  11.89 4000.0    1.0
mu[1827]  10.42    0.02   1.04   8.42   9.72  10.42  11.14  12.39 4000.0    1.0
s_mu       0.49  9.5e-3   0.06   0.35   0.45   0.49   0.53   0.61   42.0    1.1
s_Y        1.43  7.5e-3   0.06   1.33   1.39   1.43   1.46   1.55   55.0   1.07
lp__      -1158   28.55 180.58  -1459  -1281  -1171  -1060 -711.0   40.0   1.09

Rhatの値を見てみる1.1ギリギリという感じ。 n_effも50とかしかないので、厳しそう。

一応{ \displaystyle \mu}ベイズ信頼区間のグラフを見る。 白丸が晴れの日、黒丸は雨の日を表す。 濃い青は68%ぐらい、薄い青は95%ぐらい。 青線は平均値を表す。

f:id:ksknw:20161211181516p:plain

なんとなく雨の日が上手くあっていないような気がする。 雨が降ると前日との温度差が大きくなることが予想されるので、そのへんのモデルをちゃんと作ったほうが良さそう。

状態空間モデル(2階差分 降水量を考慮するモデル)

はじめの方に描いたグラフなどから天気を考慮する必要がありそう。 そこで、この題材に対するモチベーションとモデルの当てはまりから以下のようなモデルを仮定する。

{ \displaystyle
 \mu_n-\mu_{n-1} \sim \mathcal{N}(\mu_{n-1}-\mu_{n-2}, s_{\mu})
}

雨のとき

{ \displaystyle
      Y_n \sim \mathcal{N}(\mu_{n}, s_{Y|\mathrm{rain}})
}

雨がふらなかったとき

{ \displaystyle
      Y_n \sim \mathcal{N}(\mu_{n}, s_{Y|\mathrm{no\ rain}})
}

雨が振った日と振らなかった日では、真の気温から観測値へ乱数の分散が変わるというモデル 本当は分散だけでなく平均も変わったほうがいい気がするが、その変わり方が季節ごとに違うので面倒でやめた。

fit3 = fit("./model3.stan") 
fit_plot(*fit3, list_max_x=[-1,100])

f:id:ksknw:20161211181548p:plain

fit3[0]
       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu[0]     11.41    0.04   2.55   6.51   9.72  11.38  13.13   16.4 4000.0    1.0
mu[1]     10.83    0.03   1.59   7.72   9.74  10.84   11.9  13.97 4000.0    1.0
mu[2]     10.26    0.01   0.93   8.46   9.64  10.25  10.89  12.09 4000.0    1.0
mu[3]      9.71    0.01   0.72   8.32   9.22    9.7  10.18  11.11 4000.0    1.0

mu[1826]  10.52  9.9e-3   0.63   9.31  10.08  10.52  10.96  11.77 4000.0    1.0
mu[1827]  10.38    0.01   0.88   8.65   9.78  10.39  10.99  12.07 4000.0    1.0
s_mu       0.69  7.5e-3   0.07   0.55   0.64   0.69   0.74   0.82   88.0   1.03
s_Y[0]     1.06  5.5e-3   0.06   0.95   1.01   1.05    1.1   1.18  123.0   1.02
s_Y[1]      1.7  2.9e-3   0.07   1.57   1.65    1.7   1.75   1.84  583.0   1.01
lp__      -1508    13.2 122.38  -1712  -1595  -1520  -1433  -1247   86.0   1.03

前のモデルよりちゃんと収束してる感はある。 s_Y[1]が雨が降った日なので、雨の日は晴れの日よりも、真の気温からの分散が大きくなる。

前のモデルとの違いがぱっとわからないので、重ねてグラフをかく。

fig = plt.figure(figsize=(10,10))
max_x = 100
colors = ["b", "g"]
for i,ms in enumerate([fit2[1], fit3[1]]): 
    mu_mean = ms["mu"].mean(axis=0)
    mu_std  = ms["mu"].std(axis=0)

    plt.plot(mu_mean[:max_x], linewidth=0.5)    
    plt.scatter(range(len(temperatures[:max_x])), temperatures[:max_x], marker=".",# linestyle="None", 
                c=aichi["is_rain"].get_values()[:max_x].astype(float))
    
    #plt.fill_between(range(len(mu_mean[:max_x])), (mu_mean-mu_std)[:max_x], (mu_mean+mu_std)[:max_x], alpha=0.3)
    plt.fill_between(range(len(mu_mean[:max_x])), (mu_mean-mu_std*2)[:max_x], (mu_mean+mu_std*2)[:max_x], alpha=0.2
                    ,color=colors[i])

plt.xlim(0, len(temperatures[:max_x]))
plt.ylim([-5, 20])
plt.show()

f:id:ksknw:20161211141729p:plain

青が単純な2階差分、緑が雨で条件わけしたモデル。 正直良くなっているのかどうかよくわからない。

まとめと今後

アヒル本を参考にして、気温変動のモデリングを行った。 Stanの使い方はなんとなくわかってきたけど、 統計的に正しいのかとか、いまいちよくわかっていない。

今回は年ごとの関係を考慮しないでモデルを作った。 sinとかでモデル化できそうな気がする。 また、雨の日と晴れの日でsinの強さが変わるようにすると、もっといい感じにモデル化できるように思う。

Stanでまだわかってないことの1つにディリクレ混合過程がある。 アヒル本には載っていなかったけど、このへんを見ながらやってみたい。

fit.plot()のfig_sizeを大きくしたいんだけど、どうやらpython側でなくstan側からプロットされている(?)ようで、いまいちやり方がわからなかった。

(この記事には、グラフを死ぬほどかかないといけないという認識があまり反映されていないのではないか。)

参考

jupyter notebookで触れるプロットを描く

昔、Emacsから使えないからjupyter notebookは使わないという旨の記事を書いた。 実際まだ、Emacsからは使えていないけど、最近はだいたいのコードをEmacs+pythonで書いて、実行やグラフを描く部分をjupyter notebookでやるということをよくやっている。

jupyter notebook内にグラフを描く方法として

%matplotlib inline

という文を実行しておくという方法がある

大体の場合にはこれで間に合うが、3Dグラフを回したり、アニメーションを動かしたりすることができなかった。 調べてみると他にも色々できることがわかったのでまとめる。 この記事では以下のようなプロットを描く。

  • 回転できる3D
  • アニメーション
  • スライドバー付き
  • もっと綺麗なグラフを描く

markdownで出力すると全く触れなくなったので動画を作った。

www.youtube.com

回転できる3Dプロット

%matplotlib inline

でプロットを作成すると、ノートブックの中に画像として出力されてしまうので回せない。

%matplotlib notebook

を使うとノートブックの外に出力したときと同じように、インタラクティブなグラフが描ける。 こちらの3次元プロットを例として描く。

%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

x = np.arange(-3, 3, 0.25)
y = np.arange(-3, 3, 0.25)
X, Y = np.meshgrid(x, y)
Z = np.sin(X)+ np.cos(Y)

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_wireframe(X,Y,Z)

plt.show()

f:id:ksknw:20160804225927p:plain

スライドバーつき

ipywidgetsというやつを使えばできる。 jupyter notebookをインストールしただけで使おうとすると以下のコマンドを実行するように表示されたので、おとなしく従うと動くようになった。

sudo jupyter nbextension enable --py --sys-prefix widgetsnbextension

以下を実行すると、スライドバーが表示される。 このバーをいじると変数の値が動的に変わるので、コールバックにプロットを登録しておくと、スライドバーをいじるだけで勝手にプロットが変わる。

自分の環境では

%matplotlib notebook

とするとカクカクでいまいちだったが、

%matplotlib inline

にするといい感じ。

%matplotlib inline
from ipywidgets import interact
import numpy as np

def scatter(num_data):
    x = range(num_data)
    y = [np.sin(t/5.0) for t in x]
    plt.plot(x, y)
    plt.show()

interact(scatter, num_data=(1,200, 1), value=2)

f:id:ksknw:20160804225709p:plain

ipywidgetsには他にもボタンとかプログレスバーとか色々あって、ノートを簡単なアプリにできる。

%matplotlib inline

import pylab as plt
import time
from ipywidgets import FloatProgress
from IPython import display

prg = FloatProgress(min=0, max=99, value=1)
display.display(prg)

for i in range(100):
    prg.value = i
    time.sleep(0.01)

もっと綺麗なグラフを描く

matplotlibではなくて、plotlyというライブラリを使えばできる。 matplotlibと書き方が結構違うので、改めて覚えてまでやるほどかと言われると微妙だが、グラフは綺麗でかっこいい。 コードやグラフ、データをオンラインで共有したりもできるらしい。

デモのサイトを見ていると可能性を感じる。

from plotly.offline import iplot, init_notebook_mode
import plotly.plotly as py
from plotly.graph_objs import Scatter, Data


init_notebook_mode(connected=True) 

trace0 = Scatter(
    x=[1, 2, 3, 4],
    y=[10, 15, 13, 17]
)

trace1 = Scatter(
    x=[1, 2, 3, 4],
    y=[16, 5, 11, 9]
)

data = Data([trace0, trace1])
unique_url = iplot(data, filename = 'basic-line')

f:id:ksknw:20160804230515p:plain

参考

DCGANで湖岸道路ぽい画像を生成する

f:id:ksknw:20160625232128g:plain:w300

はじめに

先日、小倉であった某学会に参加した。 面白そうな話はいくつかあったけど、中でもDNNを使った画像の生成について興味を持った。 深層学習ウェイ系の某先生もマルチモーダルとか生成とか言ってた気がするし、 判別するより生成するほうが見た目に楽しそうなので、 こちらの発表でも使われていたDCGANを使って、画像を生成してみることにした。

Deep Convolutional Generative Adversarial Networks

論文はこちらとてもわかりやすいこちらのブログがわかりやすい。

GANはAdversarialという名の通り、2つのネットワークを競合させて学習を行うアルゴリズム。 GANでは普通のAutoEncoderなんかで画像を生成する場合と違って、Discriminatorというやつを作る。 Discriminatorは入力された画像が、Generatorの生成した画像か元画像かを判別する。 GeneratorはDiscriminatorを騙せるように画像の生成を学習し、Discriminatorはそれを判別できるように学習を行う。 これによって、最終的に単に入力画像だけからGeneratorを学習するよりも良くなるということらしい。

具体的には以下のような式を最適化する。

f:id:ksknw:20160626130749p:plain

右辺第1項は入力画像に対してDが出力を正しく出したときに大きくなり、第2項はDがGの画像に反応しないときに大きくなる。 Gはこの評価関数を小さくする、つまりDが騙される方に学習し、Dは逆に正しく見分けられるように学習する。

面白いのは、画像をGeneratorに直接学習させて、例えば入力画像との誤差を最小化するのではなく、 Discriminatorを騙せさえすればいいというところ(と思う)。 実際には入力画像にはなくても、Dが騙される「ぽい」画像を作れば評価関数の値を下げることができる。

DCGANはDiscriminatorにConvNet、GaneratorにDeconvNetを使う

mattyaさんによるchainerの実装

とてもおもしろい研究だし、arxivを見る感じ、最初にDCGANの論文がでたのが2015年の11月と少し古いのもあって、 既に色々なライブラリを使った実装が公開されていて、自分でコードを書く必要はなさそう。 chainerがお気に入りなので、その中から、mattyaさんによる実装を使わせてもらう。

git clone して動かそうとすると、残念ながら2GしかないGPUではメモリ足りなくて動かなかった。 コードを見ていると順次メモリを解放してやればギリギリいけそうだったので、210行目あたりに

del z
del x
del x2
del yl
del yl2
del L_dis
del L_gen

と書いてやると何とか動くようになった。

入力データ

入力データとして琵琶湖(マザーレイク)周辺の道路画像を使う。 動画を撮影しに行くのは大変なので、これを使ってGoogleMapのストリートビューから画像を取ってくる。

とってきた画像は以下のようなもの。 f:id:ksknw:20160612114759g:plain

道路に対して前だけでなく、横を向いている画像もある。 画像データは合計で600枚ぐらい。 この画像を237x119にリサイズしたあと、96x96の画像をランダムに切り出して4560枚の画像を作った。 入力画像からランダムに10枚選んだものを以下に示す。 思ったよりちゃんと道が写っているものが少ない。

f:id:ksknw:20160626005852p:plainf:id:ksknw:20160626005853p:plainf:id:ksknw:20160626005854p:plainf:id:ksknw:20160626005855p:plainf:id:ksknw:20160626005916p:plainf:id:ksknw:20160626005917p:plainf:id:ksknw:20160626005918p:plainf:id:ksknw:20160626005919p:plainf:id:ksknw:20160626010125p:plainf:id:ksknw:20160626010126p:plain

画像の中にはかなり類似したものもあるだろうし、そもそもDNNの学習に数千枚の画像というのは少ないように思うけど、JSAIの発表でもそれぐらいの枚数だったように思うので、まずはやってみる。

結果

GTX750Tiで学習。うるさいパソコンとともに数日過ごした。 10万ぐらいのグラボがほしい。

学習途中で生成された画像は以下。 上半分(50枚)の画像はエポック間で同じzベクトルで生成している画像。下半分は毎回ランダムにサンプリングしなおしたzベクトルで生成している。

最初 f:id:ksknw:20160626000823p:plain

1エポック後 f:id:ksknw:20160626000902p:plain

5エポック後 f:id:ksknw:20160626000936p:plain

10エポック後 f:id:ksknw:20160626001032p:plain

61エポック後 f:id:ksknw:20160625230705p:plain

62エポック後 f:id:ksknw:20160626003455p:plain

63エポック後 f:id:ksknw:20160626003511p:plain

大体10エポックぐらいから遠目にはそんなに変わってないように思う。 あと、上半分は同じzベクトルから生成しているはずなのに、61から63エポックで毎回結構大きく変わっているのも気になる。学習が収束している感じがあまりしない。画像が少ないとかそういうことが関係しているのだろうか。

なんか赤いなにかを作ったりしているし、まだもうちょっと学習してほしい感もある。あと1週間ぐらい回したほうがいいのかもしれない。

zベクトルをいじって色々画像を作る

DCGANでは学習時に乱数のベクトル(zベクトル)をDeconvNetworkに入力して画像を生成する。 生成した全ての画像において、Discriminatorを騙せるように学習するので、 理想的には学習後はどのzベクトルを選択しても、それっぽい画像を生成できるようになるはずである。 しかも、論文を読んでるとzベクトルを動かしていくと連続的に画像が変化していくらしい。楽しそう。

ということで適当に画像を生成した後、そこから3枚選んでそれらの間にある画像を生成してみる。 以下のように適当にプログラムを書く。

def interpolate(index1, index2, img_i):
    import pandas as pd

    epoch = 61

    all_z = pd.read_csv("z.csv").get_values()
    print all_z.shape
    sub = all_z[index2] - all_z[index1]
    gen = Generator().to_gpu()
    dis = Discriminator().to_gpu()

    serializers.load_hdf5("%s/dcgan_model_dis_%d.h5" % (out_model_dir, epoch), dis)
    serializers.load_hdf5("%s/dcgan_model_gen_%d.h5" % (out_model_dir, epoch), gen)
    gen.to_gpu()
    dis.to_gpu()
    z = xp.random.uniform(-1, 1, (batchsize, nz), dtype=np.float32)
    for i in range(100):
        z[i, :] = xp.array(all_z[index1] + sub / 100.0 * i)
    pylab.rcParams['figure.figsize'] = (0.96, 0.96)

    print z
    z = Variable(z)
    x = gen(z, test=True)
    x = x.data.get()

    for i_ in range(100):
        tmp = ((np.vectorize(clip_img)(x[i_, :, :, :]) + 1) / 2).transpose(1, 2, 0)
        pylab.imshow(tmp)
        pylab.axis('off')
        pylab.savefig('interpolate_imgs/%09d.png' % (img_i))
        img_i += 1

        pylab.clf()
    f.close()

interpolate(4, 10, 0)  # index+1がtest_imgsの画像の番号と対応する
interpolate(10, 23, 100)
interpolate(23, 4, 200)

深い理由はあまりないけど、なんとなく色々な画像を作っていて良さそうなので、61エポックのときの学習モデルを使うことにした。 以下のような3枚を選んで、それらの間のzベクトルを線形につないで、それぞれの間ごとに100枚ずつ画像を生成する。 f:id:ksknw:20160625235309p:plain:w300

生成された画像をつなげて動かしてみると以下のようになった。 f:id:ksknw:20160625232128g:plain:w300

なんか、ぽいものが生成されている。

最初の何枚かを静止画で見る。

f:id:ksknw:20160626012429p:plainf:id:ksknw:20160626012430p:plainf:id:ksknw:20160626012431p:plainf:id:ksknw:20160626012432p:plainf:id:ksknw:20160626012433p:plainf:id:ksknw:20160626012434p:plainf:id:ksknw:20160626012435p:plainf:id:ksknw:20160626012436p:plainf:id:ksknw:20160626012437p:plainf:id:ksknw:20160626012438p:plainf:id:ksknw:20160626012439p:plainf:id:ksknw:20160626012440p:plainf:id:ksknw:20160626012441p:plainf:id:ksknw:20160626012442p:plainf:id:ksknw:20160626012443p:plainf:id:ksknw:20160626012444p:plainf:id:ksknw:20160626012445p:plainf:id:ksknw:20160626012446p:plainf:id:ksknw:20160626012447p:plainf:id:ksknw:20160626012448p:plainf:id:ksknw:20160626012449p:plainf:id:ksknw:20160626012450p:plainf:id:ksknw:20160626012451p:plainf:id:ksknw:20160626012452p:plain f:id:ksknw:20160626012453p:plainf:id:ksknw:20160626012454p:plainf:id:ksknw:20160626012455p:plainf:id:ksknw:20160626012456p:plainf:id:ksknw:20160626012457p:plainf:id:ksknw:20160626012458p:plainf:id:ksknw:20160626012459p:plainf:id:ksknw:20160626012500p:plainf:id:ksknw:20160626012501p:plainf:id:ksknw:20160626012502p:plainf:id:ksknw:20160626012503p:plainf:id:ksknw:20160626012504p:plainf:id:ksknw:20160626012505p:plainf:id:ksknw:20160626012506p:plainf:id:ksknw:20160626012507p:plainf:id:ksknw:20160626012508p:plainf:id:ksknw:20160626012509p:plainf:id:ksknw:20160626012510p:plainf:id:ksknw:20160626012511p:plainf:id:ksknw:20160626012512p:plainf:id:ksknw:20160626012513p:plainf:id:ksknw:20160626012514p:plainf:id:ksknw:20160626012515p:plainf:id:ksknw:20160626012516p:plainf:id:ksknw:20160626012517p:plainf:id:ksknw:20160626012518p:plainf:id:ksknw:20160626012519p:plainf:id:ksknw:20160626012520p:plainf:id:ksknw:20160626012521p:plainf:id:ksknw:20160626012522p:plainf:id:ksknw:20160626012523p:plainf:id:ksknw:20160626012524p:plainf:id:ksknw:20160626012525p:plainf:id:ksknw:20160626012526p:plainf:id:ksknw:20160626012527p:plainf:id:ksknw:20160626012528p:plainf:id:ksknw:20160626012529p:plainf:id:ksknw:20160626012530p:plainf:id:ksknw:20160626012531p:plainf:id:ksknw:20160626012532p:plainf:id:ksknw:20160626012533p:plainf:id:ksknw:20160626012534p:plainf:id:ksknw:20160626012535p:plainf:id:ksknw:20160626012536p:plainf:id:ksknw:20160626012537p:plainf:id:ksknw:20160626012538p:plain

木が建物っぽい何かになりつつ、道の向きも気づいたら変わっている(気がする)。

まとめ

DCGANを使って、琵琶湖周辺の道路っぽい画像を生成した。 zを連続的に変化させると画像も連続的に変化していった。 ぱっと見はそれなりにできているように思うけれど、1枚1枚ちゃんと見ると、全体的にまだ線がぼやっとしていたり、自然な画像はまだ少し遠い。

誰か早くLSTM-GAN発表してコードをくれ。どうせ半年後までには誰かやるんだろ、わかってんだぞ

参考