Abnormally Distributed

統計解析担当のライフサイエンス研究者 -> データサイエンティスト@コンサル

Variational AutoEncoderのお気持ち

深層学習 変分推論

VAEの元論文を読んだ際のほぼ自分用まとめ記事。 VAE自体の分かりやすい解説記事は下記リンクを参照。

qiita.com

www.slideshare.net

1. Auto-Encoding Variational Bayes (Diederik P Kingma, Max Welling)

https://arxiv.org/abs/1312.6114
連続潜在変数を含んだ確率モデルに対する、スケーラブルな変分推論の方法の提案。
reparameterization trickの導入により、確率的勾配降下法で最適化可能な変分下限を導出した。これにより、潜在変数モデルを効率的に近似推論できるようになった。

MCMCは計算量が多く、平均場近似による変分推論では解析的に事後分布を導出することができないことが多い。 一方、Stochastic Gradient Variational Bayesは連続変数を含むほとんどどのようなモデルにも適用可能。

潜在変数zの事前分布を設定。xのzによる条件付き分布からxが生成。 いずれもパラメータθで規定された確率分布であり、pdfが微分可能とする。

f:id:kibya:20200518215705p:plain
VAEのグラフィカルモデル

推論の目的は下記の3つ。 1. パラメータθの推定 パラメータ自体に興味がある場合。また、人工データを生成したい場合。 2. 潜在変数zの事後分布の推定 暗号化やデータ表現のタスクに有用。 3. xの周辺尤度の推論 画像のノイズ除去など。

潜在変数zの事後分布の近似分布q(z|x)を考える。確率分布のパラメータφとする。 q(z|x)をエンコーダ、p(x|z)をデコーダと呼ぶ。

変分下限を変形すると、第１項のKLダイバージェンスはq(z|x)の分布に関する制約である正規化項。第２項はEncoder, Decoderを表す項であり、Encoderで推論したzの確率分布を用いて、対数尤度logp(X|z)の期待値を取っている。

f:id:kibya:20200518220306p:plain

2. Semi-Supervised Learning with Deep Generative Models (Diederik P. Kingma, Danilo J. Rezende, Shakir Mohamed, Max Welling)

https://arxiv.org/abs/1406.5298
生成モデルを用いた半教師あり学習。 ラベル割り当てのコストが高い、または不可能であり、一部のデータにしかラベルがないケースがある。ラベルがないデータも活用して、どのように分類精度を上げるかが問題。

3. Tutorial on Variational Autoencoders (Carl Doersch)

[1606.05908] Tutorial on Variational Autoencoders

1. Introduction

生成モデルは確率分布p(X)を取り扱うモデル。Xは高次元空間におけるdatapoint。
例えば画像生成の場合、datapoint Xは画像であり、Xは多数のpixelからなる高次元データ。
生成モデルの用途としては、既存のデータに似たデータを生成したいケースがある。
ある確率分布にしたがって生成されたサンプルから、なるべく元の確率分布に近いモデルを学習することで、入力と類似したサンプルを生成できる。

潜在変数モデル

0-9の数字の画像を生成する問題を考える。モデルから生成する際は、まず0-9の数字から生成対象とする数字zをランダムに決定し、その数字の画像になるように各ピクセルの値を生成する。zを潜在変数と呼ぶ。 簡単にサンプリング可能な確率分布p(z)に従う潜在変数zと決定論的な関数f(z;theta)を用いて観測変数xのサンプル（データ点）が得られる。 f(z;theta)が高確率で実際のデータに近くなるように、パラメータthetaを最適化したい。

f:id:kibya:20200531153722p:plain

f(z;theta)をP(x|z;theta)に修正すると、(1)式のP(X)を最大化する問題となる。これは最尤推定である。 VAEではP(x|z;theta)は下記のように平均f(z;theta)の正規分布とすることが多いが、例えばXが二値の場合はベルヌーイ分布など、他の分布でも良い。

f:id:kibya:20200531211011p:plain

Variational Autoencoders

(1)式を解く上では、潜在変数zをどう決めるか、zに関する積分をどう取り扱うか、の２点を決める必要がある。

画像を生成する際に、潜在変数zの各次元が表す内容や次元間の相関を明示的に決めるのは難しい。
そこで、VAEでは潜在変数zは単純な正規分布（共分散行列が単位行列）からサンプル可能とする（P(z)=N(0, I))。単純な分布からのサンプルを、複雑な関数で変換することで、モデルが必要とする潜在変数に変換できる、と仮定している。

f(z;theta)が多層ニューラルネットワークとすれば、最初の何層かで正規分布に従うzを潜在変数の値（文字の太さ、角度など）にマッピングし、残りの層で潜在変数を観測される数字の画像に変換することを想定できる。（この点については、理論的な保証はなさそう？）
実際に上記のような潜在的な構造があるかは心配する必要がなく、もしそのような構造が数字の画像の生成に役に立つならばいずれかの層で学習されるだろう。

残った問題は、zが標準正規分布に従う時、(1)式のP(x)をどうやって最大化するかである。 機械学習の分野で一般的なように、P(X)の勾配をとって、SGDで最適化できる。 そのためには、P(x)の計算式が必要になる。概念的には、zのサンプルを多数取得し、P(x|z)の期待値を取れば良いが（単純モンテカルロ法）、高次元空間においてP(X)を正確に計算するには無数のサンプルが必要になるため難しい。

実際にところ、ほとんどのzについて、P(X|z)は0に近い値になる。 VAEのキーとなる考え方は、Xを生成したであろうzのみを用いて、P(x)を計算することである。 そこで、新しい関数Q(z|X)を導入する。QはXを受け取り、Xを生成しやすいようなｚの分布を与える。 Qの元でもっともらしいzの空間は、事前分布p(z)のもとでのzの空間よりも小さいことが期待される。それにより、例えば E _ {z \sim Q} P(X \mid z)の計算が容易になる。

次に、 E _ {z \sim Q} P(X \mid z)とP(X)と関係が、変分ベイズ法のポイントとなる。

簡単な式変形により、logP(X)の下限（変分下限、ELBO）は下式で表されることが分かる。

f:id:kibya:20200531215113p:plain

左辺の第2項のKLダイバージェンスを最小化することが、ELBOの最大化と同値となる。P(z|X)は解析的に求めることはできないが、Q(z|X)を任意の精度を持つモデルで表現すると、P(z|X)に一致することが期待され、KLダイバージェンスは０となる。そのため、右辺の最大化はP(X)の最大化と同値となる。

ELBOの最大化の際は、上式の右辺をSGDにより最大化する。 右辺はencode, decoderを表す第１項と、zの近似事後分布と事前分布のKLダイバージェンスの第２項に分解できる。
zの分布系を正規分布とする場合、第２項は正規分布どうしのKLダイバージェンスは解析的に計算可能。

第１項の期待値を正確に計算するためにはzのサンプルを多数得る必要があるが、計算コストが高い。
そこで、SGDで標準的な通り、zのサンプルを1つ得て、そのzに対するP(X|z)を E _ {z \sim Q} P(X \mid z)の期待値として取り扱う。

実際は、データセットDのサンプルであるXを用いて、下式を最大化する事になる。 f:id:kibya:20200531220839p:plain

この式の勾配をとると、勾配の計算は期待値の中に移動することができる。そのため、まずXを1つサンプルし、Q(z|X)からzを1つサンプルし、下式の勾配を求めれば良い。

f:id:kibya:20200531221636p:plain

ただし、この式で表現されるネットワークは下式に左図の場合にあたる。 この場合、順伝播は問題ないが、バックプロパゲーションの際、Q(Z|X)からzをサンプルする層が非連続的となり、勾配が計算できない。
そこでreparametrization trickが必要になる。

f:id:kibya:20200531222553p:plain
　

reparametrization trickは、 N( \mu(X), \Sigma(X))を計算する際に、最初に \varepsilon \sim N(0, I)を得て、次に z = \mu(X) + \Sigma ^ {1/2}(X) * \varepsilonを計算する。 この結果、勾配を求める式は下記のようになる。上図の右側の場合に相当する。

f:id:kibya:20200531223406p:plain

EMアルゴリズムを可視化で理解

最適化

EMアルゴリズムは潜在変数を含む確率モデルのパラメータ推定に用いられる最適化手法。 有名な手法だが、あまりちゃんと理解していなかったので、PRMLで勉強してみた。

ここでは数式を書くのは辛いので、可視化でイメージを掴んでもらえればと思う。 潜在変数を含むモデルとして、混合正規分布を題材とする。

混合正規分布とは？

複数の正規分布の重み付き和が混合正規分布となる。ここで重み\pi_kは混合係数と呼ばれ、 \sum ^ {K} _ {k=1}\pi _ {k}=1である。

\displaystyle{p(\boldsymbol x) =\sum ^ {K} _ {k=1}\pi _ {k} \mathcal{N}( \boldsymbol x | \boldsymbol \mu _ {k}, \boldsymbol \Sigma _{k})}

下図は、３つの二変量正規分布からなる混合分布の例。それぞれの分布の平均を✖️で、共分散行列を95%信頼区間を示す楕円で表現している。

f:id:kibya:20190910214621p:plain
混合正規分布

EMアルゴリズムはこの混合分布のパラメータ、すなわち平均、共分散行列、混合係数をデータから推定するための方法である。
つまり、下図のようなラベルのないデータが与えられた際に、上述のパラメータを求める際に使う。なお、混合した分布の数は既知とする。

f:id:kibya:20190910214833p:plain
混合正規分布（ラベルなし）

パラメータの最尤推定を行う際は、対数尤度関数をパラメータで偏微分して0とおいた式を解くのが基本的なアプローチだが、潜在変数を含むモデルでは、解析的にこの式を解くことができない。そこでEMアルゴリズムの出番となる。

EMアルゴリズムの流れ

  1. 平均、共分散行列、混合係数をランダムに初期化。
  2. Eステップ：負担率(responsibility)を算出。
  3. Mステップ：潜在変数を固定した上で、最尤法により平均、共分散行列、混合係数を新たに算出。
  4. 対数尤度が収束するまで上記のE Step, M Stepを繰り返す。

順を追って説明する。

1. パラメータの初期化

ここでは平均は一様分布からサンプリング、共分散行列は0.2I、混合係数は\alpha=1000のディリクレ分布からサンプリングした。

f:id:kibya:20190910220353p:plain
パラメータの初期化

2. E step

E Stepでは潜在変数の事後確率、すなわち負担率(responsibility) \gamma(z _ {i,k})を算出する。 ここで\gamma(z _ {i,k}) はデータi がクラスターkから得られた確率を表す。

下図は\gamma(z _ {i,0}), \gamma(z _ {i,1}), \gamma(z _ {i,2}) をRGBとして各データ点に色をつけている。
つまり赤色の点は赤い✖️のクラスターに属する確率が高く、紫の点は赤・青のクラスターに属する確率が半々くらい、と言うことになる。

負担率の計算には初期化したパラメータを用いる。2回目以降のE StepではM Stepで更新したパラメータを用いる。

f:id:kibya:20190910220413p:plain
E Step

3. M step

M Stepでは最尤法により平均、共分散行列、混合係数の各パラメータを算出する。
詳細な説明は省くが、混合正規分布の対数尤度を各パラメータで微分して0とおくと、各パラメータの最尤推定値は負担率を用いて表せることが分かる。
つまりE Stepで得られた負担率を用いて、より良いパラメータを再計算できることになる。

その様子を表したのが下図で、平均、共分散行列が修正されたことが分かる。

f:id:kibya:20190910220435p:plain
M Step

4. E Step, M Stepの繰り返し

上記のE Step, M Stepを繰り返すと、対数尤度が常に増加することを示せる。
実用上は、対数尤度の増加があるしきい値より小さくなった時点で、アルゴリズムが収束したと判定する。

そもそもなぜE Step, M Stepに分ける必要があるかと言うと、E Stepでは負担率の計算にパラメータが必要、M Stepではパラメータの計算に負担率が必要だからである。
言い換えると、「パラメータを固定して負担率を計算 → 負担率を固定してパラメータを計算」を繰り返して最適化していく方がEMアルゴリズムだと言える。



EMアルゴリズムが収束する様子をアニメーションで示す。この場合は最初の数回でほぼ最適解が得られていることが分かる。

f:id:kibya:20190910230921g:plain
EMアルゴリズム　アニメーション

今回の場合、対数尤度の変化量のしきい値を0.001とすると64 iterationで収束した。
対数尤度は最初のステップで大幅に増加し、以降は小刻みに増え続けている。

f:id:kibya:20190910224010p:plain
対数尤度



結び

EMアルゴリズムは古典的で重要な手法だが、最尤推定には特異性の問題など重大な限界がある。
ベイズ的なアプローチである変分推論では最尤推定が抱える問題を解決できるとのことなので、次回はこちらを解説したい。

今回作成したコードはこちら。

import numpy as np
from numpy import random
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter
import matplotlib.animation as ani
from matplotlib.patches import Ellipse
%matplotlib inline


# draw 95% confidence interval of multivariate normal distribution
def draw_ellipse(mu, cov, ax, **kwargs):
    var, U = np.linalg.eig(cov)
    if U[1, 0] == 0:
        angle = 0
    else:
        angle = 180. / np.pi * np.arctan(U[0, 0]/ U[1, 0])

    e = Ellipse(mu, 2 * np.sqrt(5.991 * var[0]), 2 * np.sqrt(5.991 * var[1]), 
                angle=angle, facecolor='none', **kwargs)
    ax.add_artist(e)
    return e


def make_plot(data, mu, cov, scatter_color, ax):
    artists = []
    im1 = ax.scatter(data[:,0], data[:,1], s=20, c=scatter_color, alpha=0.5)
    artists.append(im1)
    for i in range(n_cluster):
        im2 = ax.scatter(mu[i, 0], mu[i, 1], marker="x", s=50, c=colors[i])
        artists.append(im2)
        im3 = draw_ellipse(mu[i], cov[i], ax, edgecolor=colors[i])
        artists.append(im3)
    return artists


def calc_L(data, mu, cov, pi):
    L = np.zeros((len(data), len(pi)))
    for i in range(len(pi)):
        L[:,i] = pi[i] * stats.multivariate_normal(mean=mu[i], cov=cov[i]).pdf(data)
    return L


#＃# サンプルデータ作成
n_cluster = 3
n_sample = 500
z_list = np.array([0.2, 0.5, 0.3])
colors = ["r", "g", "b"]

random.seed(18)

mu_list = stats.uniform(loc=-3, scale=6).rvs(size=(3, 2))
cov_list = stats.invwishart(df=3, scale=np.eye(2)).rvs(3)

# 各サンプルをどの分布から発生させるか
lab = random.choice(n_cluster, n_sample, p = z_list)
class_count = Counter(lab)

data = []
for i in range(n_cluster):
    y = stats.multivariate_normal(mean=mu_list[i], cov=cov_list[i]).rvs(class_count[i])
    data.append(y)
    
D = np.vstack(data)  
    

#＃# EMアルゴリズム
threshold = 1e-4
fig, ax = plt.subplots()
ax.set_xlim(-3, 6)
ax.set_ylim(-5, 4)
ims = []

## mu, Sigma, piの初期化
(xmin, ymin), (xmax, ymax) = D.min(axis=0), D.max(axis=0)

random.seed(0)
mu_x = random.uniform(xmin, xmax, n_cluster)
mu_y = random.uniform(ymin, ymax, n_cluster)
mu = np.c_[mu_x, mu_y]

cov = np.tile(0.1 * np.eye(2), (n_cluster, 1, 1))
pi = random.dirichlet(1000 * np.ones(n_cluster))

## E Step 
L = calc_L(D, mu, cov, pi)

# log likelihood
loglike = np.log(L.sum(axis=1)).sum()


for step in range(50):
    loglike_prev = loglike
    
    ## E Step 
    # responsibility
    gamma = L / L.sum(axis=1, keepdims=True)

    ## plot
    artists = make_plot(D, mu, cov, gamma, ax)
    ims.append(artists) 
    
    ## M Step
    # 各クラスターの実質的なデータ数
    n_k = gamma.sum(axis=0)

    # 負担率による重み付き和
    mu = np.zeros((n_cluster, 2))
    for k in range(n_cluster):
        mu[k] = (gamma[:,k, np.newaxis] * D).sum(axis=0) / n_k[k]

    # 共分散行列
    cov = np.zeros((n_cluster, 2, 2))
    for k in range(n_cluster):
        diff = D - mu[k]
        sigma = np.dot(diff.T, gamma[:, k, np.newaxis] * diff) / n_k[k]
        cov[k] = sigma

    # 混合係数
    pi = n_k / n_sample
    
    ## calculate log likelihood
    L = calc_L(D, mu, cov, pi)
    loglike = np.log(L.sum(axis=1)).sum()

    loglike_diff = loglike - loglike_prev
    
    ## plot
    artists = make_plot(D, mu, cov, gamma, ax)
    ims.append(artists) 
    
    print("decrease of log likelihood: {0:.4f}".format(loglike_diff))
    
    if loglike_diff < threshold:
        print(f"log likelihood converged after {step} iteration")
        break
        
# アニメーションを作成する。
plt.tight_layout()
anim = ani.ArtistAnimation(fig, ims, interval=100, repeat_delay=1000)
# anim.save('animation.gif', writer='imagemagick')

参考文献

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

はじめてのパターン認識

はじめてのパターン認識

  • 作者:平井 有三
  • 発売日: 2012/07/31
  • メディア: 単行本（ソフトカバー）

Wishart分布の可視化

確率分布

ベイジアンモデリングで多次元正規分布の共分散行列の事前分布を設定することがある。

その際、よく取り上げられるのがWishart分布である。Wishart分布は多次元正規分布の精度行列（分散共分散行列の逆行列）に対する共役事前分布になることが知られている。分散共分散行列の事前分布としては、InverseWishart分布を用いれば良い。

WIshart分布の確率密度関数は下記の通り。パラメータ\nuは自由度であり、\nu > D-1である必要がある。 また、パラメータ\Sigmaは正定値対称行列であり、scipyの実装ではscale matrixと呼ばれている。

{ \displaystyle \text{Wishart}(S|\nu,\Sigma) =C_w(\nu, \Sigma) \left| S \right|^{(\nu - K - 1)/2} \ \exp \left(- \frac{1}{2} \ \text{tr}\left( \Sigma^{-1} S \right) \right)}

Wishart分布の期待値は下記の通り。

E(S)=\nu\Sigma

(Inverse)Wishart分布は行列を生成する確率分布なので、直感的にイメージしづらい。 そこで、Wishart分布から得られたサンプルを精度行列とする２次元正規分布について、95%信頼区間を図示して可視化を試みる。

50個の精度行列をサンプリングして、それぞれについて95%信頼区間を赤線で示している。 また、青線は精度行列の期待値を用いた場合の95%信頼区間である。  

まずは\Sigmaを単位行列として、自由度を変化させた場合。

f:id:kibya:20190816160824p:plain
自由度nuを変化させた場合のWishart分布
 

次に、自由度を固定して\Sigmaを定数倍して行った場合。

 

f:id:kibya:20190816161042p:plain
scale matrix (sigma)を変化させた際のWishart分布

\Sigmaを非対角行列とすることで、相関がある共分散行列を得られやすくすることもできる。

f:id:kibya:20190816161436p:plain
期待値が相関係数が0.8の共分散行列となるWishart分布

f:id:kibya:20190816161221p:plain
期待値が相関係数が-0.4の共分散行列となるWishart分布

ただし、PyMC3やStanでは共分散行列の事前分布としては、LKJ分布の利用が勧められている。

PyMC3のissueにも挙げられている通り、

  1. 正定値対称行列の制約があるため、全ての成分を少しずつ動かしたサンプルが有効となる確率はほぼ０となる。最初に制約のない空間に変換した上でサンプリングし、制約のある空間に戻した方が好ましい。
  2. Wishart分布はかなり裾が重い分布であり、サンプリング効率が悪くなってしまう。LKJ分布はWishart分布ほど裾が重くない。
  3. LKJ分布の方がWishart分布よりも直感的に解釈しやすい。

とのことらしい。 LKJ分布についてはまた記事を書く。

プロット出力に用いたPythonコードは下記の通り。

from matplotlib.patches import Ellipse
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np
from scipy import stats


# draw 95% confidence interval of multivariate normal distribution
def draw_ellipse(cov, ax, **kwargs):
    var, U = np.linalg.eig(cov)
    if U[1, 0]:
        angle = 180. / np.pi * np.arctan(U[0, 0]/ U[1, 0])
    else:
        angle = 0

    e = Ellipse(np.zeros(2), 2 * np.sqrt(5.991 * var[0]), 2 * np.sqrt(5.991 * var[1]), 
                angle=angle, facecolor='none', **kwargs)

    ax.add_artist(e)
    ax.set_xlim(-8, 8)
    ax.set_ylim(-8, 8)

        
# visualize samples from Wishart distribution
def visualize_Wishart(nu, sigma, n_sample, ax):
    delta_samples = stats.wishart(df=nu, scale=sigma).rvs(n_sample)
    
    for delta in delta_samples:
        cov = np.linalg.inv(delta)
        draw_ellipse(cov, ax, edgecolor='r', linewidth=0.5, alpha=0.4)
    # 期待値
    draw_ellipse(np.linalg.inv(nu * sigma), ax, edgecolor='b', linewidth=2.0, linestyle='--', alpha=0.7)
    ax.set_title(f'df={nu}, scale={sigma.round(2)}')


## 様々なパラメータで可視化
# nuを変化させる
nu_list = [2, 3, 4]
sigma = np.eye(2)

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

np.random.seed(2)
for i, nu in enumerate(nu_list):
    visualize_Wishart(nu, sigma, 50, axes[i])
plt.tight_layout()
plt.savefig('Wishart_sample_nu.png')


# sigmaを変化させる
scale_list = [2, 3, 4]
sigma = np.eye(2)
nu = 3

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

np.random.seed(2)
for i, scale in enumerate(scale_list):
    visualize_Wishart(nu, scale * sigma, 50, axes[i])
plt.tight_layout()
plt.savefig('Wishart_sample_sigma.png')


# 相関係数0.8
nu_list = [2, 3, 4]
sigma = np.linalg.inv(np.array([[1.0, 0.8], [0.8, 1.0]]))

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

np.random.seed(2)
for i, nu in enumerate(nu_list):
    visualize_Wishart(nu, sigma, 50, axes[i])
plt.tight_layout()
plt.savefig('Wishart_sample_nu_corr0.8.png')


# 相関係数-0.4
scale_list = [2, 3, 4]
nu = 3
sigma = np.linalg.inv(np.array([[1.0, -0.4], [-0.4, 1.0]]))

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

np.random.seed(2)
for i, scale in enumerate(scale_list):
    visualize_Wishart(nu, scale * sigma, 50, axes[i])
plt.tight_layout()
plt.savefig('Wishart_sample_sigma_corr-0.4.png')

「基礎からのベイズ統計学」の例を実装しながら理解する（①MH法まで）

ベイジアンモデリング MCMC

「基礎からのベイズ統計学」で紹介された例をPythonで実装してみる。この記事は基本的には備忘録であり、書籍で詳しく解説してある部分は省略し、ポイントと実装して確認した結果のみ記載する。

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門


MCMCを使う理由

何はともあれ、ベイズの定理。

\displaystyle{f(\theta|x)=\frac {f(x|\theta)f(\theta)} {\int f(x|\theta)f(\theta)d\theta}}

事後分布の計算には分母の積分が必要になるが、解析的に解けないことがほとんど。
そこで、分子の（事前分布）x （尤度）の部分（本書でいう事後分布のカーネル）を利用して、母数の事後分布に従う乱数を発生させることで、母数に関する推論を行うことを考える。そのための手法がMCMC（マルコフ連鎖モンテカルロ法）である。

 

詳細釣り合い条件

\displaystyle{f(\theta|\theta')f(\theta')=f(\theta'|\theta) f(\theta)}

MCMCが定常分布に収束するための十分条件が詳細釣り合い条件である。これが、あらゆる2点間に成り立っている必要がある。
上式は、\thetaと\theta'の確率の比が1:pならば、\thetaから\theta'への移動は\theta'から\thetaへの移動のp倍起こりやすい、ということを意味する。

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

MCMCのうち、一番単純な手法。適当な提案分布をマルコフ連鎖の遷移核（f(\theta|\theta')）として利用する。詳細釣り合い条件を満たすように遷移核を補正していくことで、事後分布に従うサンプルを得られるようになる。

独立MH法

本書の通り、「波平釣果問題」を題材とする。

波平さんの過去10回の釣果は0匹が6回、1匹が3回、2匹が1回だった。波平さんが次に釣りに行くときに、釣果が0匹となる確率を求めよ。ただし、釣果はポアソン分布に従っているとする。また、事前分布は\alpha=6, \beta=3のガンマ分布とする（平均2、標準偏差0.8程度）。

この場合、ポアソン分布の共役事前分布であるガンマ分布を使っているので、事後分布もガンマ分布となり、解析的に確率を求めることが可能である。ここでは、サンプリング結果の正しさを確認するために、敢えてこの問題を題材としている。


独立MH法では、提案分布を無条件分布とする。つまり、提案される候補値は互いに独立となる。
本書では平均1.0 、分散0.5の正規分布を提案分布としている。この場合、負の値がサンプリングされることがあるが、ポアソン分布の母数は正の値である必要があるため、不適である。今回の実装では提案された値が負の場合は、事後分布のカーネルが0となるようにしている。こうすることで、負の値が提案された場合は必ず受容されないことになる。

独立MH法の欠点として、適切な提案分布を選ばないと収束に時間がかかる、または現実的な時間では収束しないことが挙げられている。
提案分布を正規分布として、平均は1.0に固定、標準偏差を0.5, 1.5, 0.1とした場合でサンプリングを試みる。

f:id:kibya:20190331110425p:plain
目標分布と提案分布


それぞれについて、6000個のサンプリングを行い、最初の1000個をバーンイン期間として破棄して、残りの5000個でヒストグラムを作成した結果を下に示す。
受容率はそれぞれ、53.7%, 20.1%, 39.7%となった。

提案分布の幅が広すぎる場合（標準偏差：1.5）は、受容率が低くなってしまう。また、幅が狭すぎる場合（標準偏差：0.1）は目標分布からの満遍ないサンプルを得ることが難しく、収束させることは難しい。

f:id:kibya:20190331110504p:plain
独立MH法によるサンプリング結果

ランダムウォークMH法

 独立MH法では適切な提案分布の設定が難しいという問題があった。ランダムウォークMH法では候補の提案にランダムウォークを利用することで、この問題を解決している。提案分布に対象な分布を選ぶと、受容確率（補正係数）の式から提案分布の部分が消え、事後カーネルの比となる。

提案分布を \mu=0の正規分布として、 3通りの\sigmaでサンプリングを行なった結果が下記。6000個のサンプリングを行い、最初の1000個はバーンイン期間として破棄している。受容率はそれぞれ、49.3%, 92.2%, 14.6%となっている。提案分布の分散 \sigma_\epsilonが適切な値であれば、独立MH法よりも受容率がかなり高くなることが分かる。
 

f:id:kibya:20190401092145p:plain
ランダムウォークMH法

次に、トレースプロット（各ステップで採択されたパラメータを順番にプロットしていったもの）を見て、それぞれの\sigma_\epsilonの値の場合で何が起こっているかを確認する。なお、提案された値のうち、破棄されたものを赤点で示している。

f:id:kibya:20190401092219p:plain
ランダムウォークMH法　トレースプロット

 \sigma_\epsilon=0.5は適切な値であり、少し離れた所にあった初期値から速やかに0.8付近に近づいていることが分かる。  一方、 \sigma_\epsilon=0.05と小さすぎる場合は、受容率こそ高いものの、サンプルされる値が少しずつしか変化しないため、収束するには時間がかかる。また、 \sigma_\epsilon=2.0と大きすぎる場合は、受容率が低くなってしまい、やはりサンプリングの効率が悪くなってしまう。

今回のようにパラメータが1つの場合は\sigma_\epsilonの調整は難しくないが、パラメータ数が多くなると\sigma_\epsilonの調整が難しくなり、受容率が低下するという問題がある。
これを解決するための手法が、ハミルトニアンモンテカルロ法（HMC法）である。次回は、こちらを実装してく。


Pythonコードはこちら。Jupyter notebookを前提に書いている。

%matplotlib inline
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import cm
import seaborn as sns
from tqdm import tqdm_notebook
from scipy.stats import norm, poisson, gamma, uniform
sns.set(style="white")

np.random.seed(1)

x = np.linspace(-1, 5, 100)
plt.plot(x, norm.pdf(x, loc=1.0, scale=0.5), label="proposal (mu=1.0, sigma=0.5)")
plt.plot(x, norm.pdf(x, loc=1.0, scale=1.5), label="proposal(mu=1.0, sigma=1.5)")
plt.plot(x, norm.pdf(x, loc=1.0, scale=0.1), label="proposal(mu=1.0, sigma=0.1)")
plt.plot(x, gamma.pdf(x, a=11, scale=1./13), linestyle="--", label="posterior")
plt.xlabel("theta")
plt.ylim(-0.1, 2)
plt.legend()
plt.tight_layout()
plt.savefig("proposal_dists.png")


lam = 2
data = [0, 0, 0, 0, 0, 0, 1, 1, 1, 2]


def posterior_kernel(lam):
    if lam < 0:
        return 0
    likelihood = poisson.pmf(data, mu=lam).prod()
    prior = gamma.pdf(lam, a=6, scale=1./3)
    return prior * likelihood


# 独立MH法
def independent_MH(mu_proposal, scale_proposal, ax, theta_init=1.0, n_sample=6000, burn=1000):
    theta_current = theta_init
    posterior = [theta_current]
    count = 0 # count accepted proposals
    
    for i in tqdm_notebook(range(n_sample)):

        proposal_dist = norm(loc=mu_proposal, scale=scale_proposal)
        theta_proposal = proposal_dist.rvs()

        # Lieklihood x prior
        p_current = posterior_kernel(theta_current)
        p_proposal = posterior_kernel(theta_proposal)

        # probability to accept proposal
        p_accept = proposal_dist.pdf(x=theta_current) * p_proposal / \
        (proposal_dist.pdf(x=theta_proposal) * p_current)

        if np.random.rand() < p_accept:
            count += 1
            theta_current = theta_proposal

        posterior.append(theta_current)
    
    g = sns.distplot(posterior[burn:], ax=ax, label="posterior samples")
    g.plot(x, gamma.pdf(x, a=11, scale=1./13), label="theoretical posterior")
    g.legend()
    g.set_title(f"proposal mu: {mu_proposal}, sigma: {scale_proposal}")
    g.set_xlim(0, 2.5)
    print("accept ratio: {:.2%}".format(count / n_sample))


fig, axes = plt.subplots(1, 3, figsize=(14, 4))
independent_MH(1.0, 0.5, axes[0])
independent_MH(1.0, 1.5, axes[1])
independent_MH(1.0, 0.1, axes[2])
plt.savefig("independent_MH.png")


# ランダムウォークMH法
def randomwalk_MH(scale_proposal=0.5, theta_init=4.0, n_sample=6000, burn=1000):
    theta_current = theta_init
    posterior = [theta_current]
    rejected = []
    count = 0 # count accepted proposals
    
    for i in tqdm_notebook(range(n_sample)):

        # 前のサンプルに基づき、次の候補が提案される
        theta_proposal = norm.rvs(loc=theta_current, scale=scale_proposal)

        # Lieklihood x prior
        p_current = posterior_kernel(theta_current)
        p_proposal = posterior_kernel(theta_proposal)

        # probability to accept proposal
        p_accept = p_proposal / p_current

        if np.random.rand() < p_accept:
            count += 1
            theta_current = theta_proposal
        else:
            rejected.append((i, theta_proposal))

        posterior.append(theta_current)
        
    print("accept ratio: {:.2%}".format(count / n_sample))
    return np.array(posterior), np.array(rejected)


def plot_posterior(posterior, burn=1000):
    g = sns.distplot(posterior[burn:], label="posterior samples")
    g.plot(x, gamma.pdf(x, a=11, scale=1./13), label="theoretical posterior")
    g.legend()
    g.set_xlim(0, 2.5)
    

def plot_trace(posterior, rejected):
    plt.plot(posterior, label="trace")
    plt.scatter(rejected[:, 0], rejected[:, 1], marker="o", color="r", s=4, label="rejected")
    plt.legend()


sigma1, sigma2, sigma3 = 0.5, 0.05, 2.0
posterior1, rejected1 = randomwalk_MH(sigma1)
posterior2, rejected2 = randomwalk_MH(sigma2)
posterior3, rejected3 = randomwalk_MH(sigma3)

fig, axes = plt.subplots(figsize=(14, 4))
plt.subplot(131)
plot_posterior(posterior1)
plt.title(f"proposal sigma: {sigma1}")

plt.subplot(132)
plot_posterior(posterior2)
plt.title(f"proposal sigma: {sigma2}")

plt.subplot(133)
plot_posterior(posterior3)
plt.title(f"proposal sigma: {sigma3}")

plt.savefig("randomwalk_MH.png")


fig, axes = plt.subplots(figsize=(10, 8))
plt.subplot(311)
plot_trace(posterior1, rejected1)
plt.title(f"proposal sigma: {sigma1}")

plt.subplot(312)
plot_trace(posterior2, rejected2)
plt.title(f"proposal sigma: {sigma2}")

plt.subplot(313)
plot_trace(posterior3, rejected3)
plt.title(f"proposal sigma: {sigma3}")

plt.tight_layout()
plt.savefig("randomwalk_MH_trace.png")

行列のランクが意味すること（Pythonで可視化）

線形代数

最近、こちらの本を参考に、線形代数を復習している。

プログラミングのための線形代数

プログラミングのための線形代数


行列は「写像」であるという観点で、行列式、固有値などが空間上で何を意味するかを説明していて、個人的には目から鱗、という感じだった。線形代数は大学の教養過程や大学院の授業でも習ったけど、こういう観点の解説があれば良かったなと思う。

 
また、プログラミングのためのと銘打っているだけに、逆行列や固有値などをコンピュータで効率よく計算するための手法に関する解説がある点も、普通の数学の教科書とは差別化されている。

 
今回は、こちら本を読んで、行列のランクが意味することをよく理解できたので、備忘録を書いておく。
行列のランクは重回帰分析や主成分分析のような基本的なデータ解析手法とも関わってくる大事な概念である。

ここでは、視覚的に把握するため、2x2正方行列Aによる写像を考え、ランクが何を表すかを説明する。

A=\begin{pmatrix} 1 & 2 \\ 3 & 1 \\ \end{pmatrix}
 
この行列Aによる写像は、点(1, 0) が (1, 2) に、点 (0, 1) が (3, 1) に映るような写像となる。
ここで、 下図左のような原点を中心とした正方形上のデータ点を考えると、行列Aによる写像は、中央の図のようになる。回転や伸縮の結果、歪んだ四角形となることがわかる。

f:id:kibya:20190203102348p:plain
行列Aによる写像


また、歪んだ四角形上の点を元の正方形上の点に戻す写像は、Aの逆行列による写像である。
PAx = x \Leftrightarrow PA = I \Leftrightarrow P = A^{-1}
 

さて、ここから本題の行列のランクについて説明する。
行列Aのランク（rank(A)）とは「行列Aの行ベクトル（または列ベクトル）のうち、1次独立なものの数」である。
ベクトルの１1次独立というのは、他のベクトルを定数倍したり足し合わせたりして表すことができないベクトルの組を意味する。


2次の正方行列の場合、行は2つしかないので、当然ランクは最大で2となる。
実は上の行列Aはランクが2の場合である。
n次正方行列のランクがnのとき、この行列は正則行列と呼ばれ、逆行列が存在する。

では、2次の正方行列でランクが1の場合はどうなるか。
例えば、下記の行列Bによる写像を考える。

B=\begin{pmatrix} 1 & 2 \\ 0.5 & 1 \\ \end{pmatrix}
 
行列Bの2列目は1列目の2倍として表すことができる。つまり、行列Bの列のうち、1次独立なものは1つしか存在しないので、rank(B)は1である。
行列Aの場合と同じように、正方形上の点のBによる写像を確認すると、すべての点が1つの直線上に乗ってしまう。
実はこの正方形上に限らず、平面状のどの点もこの直線上に写像される。

f:id:kibya:20190203104217p:plain
行列Bによる写像


Bによる写像ではすべての点が１直線上に乗ってしまうので、元の２次元空間の情報は失われてしまう。
つまり逆写像は存在せず、逆行列も存在しないということになる。

このように、n次正方行列Aについて、rank(A) < nとなることをランク落ちという。
 ランク落ちした行列による写像は、元のn次元空間を「つぶして」しまい、rank(A)次元に落とすような写像となる。
このとき、Aの逆行列は存在しない。


まとめると、n次正方行列Aについて、下記の5つの表現は同じことを意味する。

 

Pythonコードはこちら。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
%matplotlib inline

# 正方形上の点
x = np.linspace(-1, 1, 9)
tmp1 = np.array([(i, j) for i in [-1, 1] for j in x]) 
tmp2 = np.array([(j, i) for i in [-1, 1] for j in x]) 
data = np.vstack([tmp1, tmp2])

lab = np.abs(data.sum(axis=1)) == 1

plt.figure(figsize=(14, 5))
plt.subplot(131, aspect=1.0)
plt.scatter(data[:,0], data[:,1], s=10, color="gray")
plt.scatter(data[lab, 0], data[lab, 1], marker="s", s=20, c=["red", "blue", "green", "orange"])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

# 行列Aによる写像
A = np.array([[1, 2], [3, 1]])
data2 = np.dot(data, A)

plt.subplot(132, aspect=1.0)
# plt.axes().set_aspect('equal')
plt.scatter(data2[:,0], data2[:,1], s=10, color="gray")
plt.scatter(data2[lab, 0], data2[lab, 1], marker="s", s=20, c=["red", "blue", "green", "orange"])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

# 逆行列=逆写像
A_inv = np.linalg.inv(A)
data_inv = np.dot(data2, A_inv)

plt.subplot(133, aspect=1.0)
plt.scatter(data_inv[:,0], data_inv[:,1], s=10, color="gray")
plt.scatter(data_inv[lab, 0], data_inv[lab, 1], marker="s", s=20, c=["red", "blue", "green", "orange"])
plt.xlim(-5, 5)
plt.ylim(-5, 5)


B = np.array([[1, 2], [0.5, 1]])
data3 = np.dot(data, B)

plt.figure(figsize=(14, 5))
plt.subplot(131, aspect=1.0)
plt.scatter(data[:,0], data[:,1], s=10, color="gray")
plt.scatter(data[lab, 0], data[lab, 1], marker="s", s=20, c=["red", "blue", "green", "orange"])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

plt.subplot(132, aspect=1.0)
plt.scatter(data3[:,0], data3[:,1], s=10)
plt.scatter(data3[lab, 0], data3[lab, 1], marker="s", s=20, c=["red", "blue", "green", "orange"])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

 ちなみに、行列のランクはnumpy.linalg.matrix_rankで算出することができる。

2017転職活動まとめ①〜エージェント〜

初めての転職活動がほぼ終わったので、まとめてみる。

現職のメーカーＴには2016年4月に入社した。入社前から元々、ちょっと違うかもなーとは思っていて、新人導入研修でこの会社無理そうだと悟った。笑

 

転職活動を始める、となるとまずはエージェントの登録を考えるかと思う。そこでまず、自分が使ったエージェントのまとめ。

 

エージェント

結局５社に登録したが、こんなに必要はなかった。

登録した順に、DODA、JAC、ムービン、マイナビ、リクルート。

担当者との相性もあると思うので、何社か登録するべきだとは思う。自分の希望条件を話しても、先方の受け止められ方が違ったり、チョイスしてくる求人も異なるので、まずは色々なエージェントから紹介してもらうのが良いと思う。

ただ、その後応募する際にはなるべく１社にまとめた方が良い。そのほうが複数の会社を同時に進める場合の面接の日程調整がスムーズ。そして何より、複数社から内定をもらって、どこに入ろうか迷う段階になった際、各社の応募に使ったエージェントが異なると公平な意見が得られない、というのが大きなデメリット。自分は完全にこれに陥った。

なお、自分は職場も住居も東京からやや離れていることもあり、キャリアアドバイザーとはほぼ電話でしか話さなかった。キャリアアドバイザーは朝は遅いが夜は大抵21時過ぎまでいる。なので仕事終わりとかに普通に電話できる。昼休みに電話することもあった。

リクルートエージェントだけはだいぶ休日が多かったが、他社は休日が少なく、なかなか大変そうだった。なお、ムービン以外は女性のキャリアアドバイザーだった。

 

それぞれの寸評。

DODA：最初に登録し、そのままメインで使うことになった。登録したきっかけは覚えていないが、新卒のときマイナビのキャリアフェアやリクナビのエージェントがイマイチすぎてネガティブなイメージを持っていたからだと思う。

担当キャリアアドバイザーがつき、面談で話した希望条件に合わせて求人を紹介してくれる。これとは別に採用プロジェクト担当という、特定の企業の担当者がメールを送りつけてくる。これが１日に複数届くくらいの多さで、結局は迷惑メールフィルタで対処した。キャリアアドバイザーがついている期間はこのメールが止められないシステムらしい。

また、ここだけ先方のオフィスを訪問して、面接対策を受けた。面接対策の専門家のおばさんで、「転職先で〜を学びたいではなく、自分は〜ができると伝える」など、当たり前だが面接で大事そうなことを改めて確認できる、結構良い機会だったと思う。

ちなみにDODAとは別にエグゼクティブサービス？というのを同じ会社がやっているようで、製薬担当のリクルーターから直接連絡があった。良い人そうだったが、だいぶ忙しそうで結局連絡は途絶えたまま。

 

JAC：外資に強いというイメージで登録。JAC internationalというのも同時に登録した。ここは担当のキャリアアドバイザーは特につかず、最初の面談の情報を共有して企業担当者がこちらにアプローチしてくるシステムらしい。たまに微妙な会社を紹介されるくらいで、結局ほとんど使わなかった。もう少し職務経験がある人向けだったのかもしれない。

 

ムービン：コンサルを考えていたこともあり、コンサル業界No.1との触れ込みを見て登録した。こちらは唯一、登録時にキャリアアドバイザーと対面で面談した。エージェントは東大卒の方で、週末の深夜2時など結構滅茶苦茶な時間にメールを送ってくることもあり、だいぶ忙しそうな方だった。さすがにコンサルに特化しているだけあり、受けた企業の面接担当者などしっかり把握していた。

 

マイナビ：担当エージェントのＴさんがとても良い方だった。応募を進めている段階で登録したので、結局こちら経由ではメーカーＨしか受けられなかったのが残念。良い方だが、紹介してくる求人はちょっとずれているものも多かった。あと、登録しても自分で求人検索ができないのが残念なところ。

 

リクルート：キャリアアドバイザーとアシスタントの二人体制だった。他社よりアドバイザーの休日が多いのか、担当者がたまたまそういう勤務体系だったのか。マイページから自分で求人検索ができたり、求人票に募集人数が書いてあったり、お気に入り登録している求人の応募状況がたまにメールで送られてきたり、エージェントレポートなる企業担当者からの紹介文書が一部求人についていたりなど、なかなか参考にさせて頂いた。こちらも選考が進んでからの登録で、ITのLとメーカーDしか受けなかった。さすがに最大手なだけはあり、とりあえず登録するならここが良さそう。

 

また、外資に行きたい気持ちもあり、元々作ってはいたLinkedInのプロフィールをちゃんと書き、転職シグナルを発信する設定にした。すると、４ヶ月で10人くらいのリクルーターからメッセージが来た。

大体は小さい外資系の人材紹介会社。ただ、日系化学メーカーNからはリクルーターから直接メッセージが来た。AIの活用を進める組織を新設したので、その中心的な人材として来てほしいとのこと。美女だったので少し迷ったが、結局返事はしなかった。現職の競合他社で、雰囲気も似たようなもので自分には合わないと考えた。

また、LinkedIn経由で外国人のリクルーターと一回だけ、恵比寿のオフィスで直接面接したがイマイチだった。ヘルスケア関連のデータ解析の仕事がしたいと言っているのに、微妙そうな外資の技術営業のようなポジションばかり紹介された。

また、別のリクルーターから、日本に進出しようとしている画像解析系のベンチャーのポジションも紹介され、電話で話した。英語で電話はなかなか厳しい。こちらも、会社自体は面白そうだったが、Sales engineerのポジションだったこと、そもそもベンチャーに行くのは周囲の強い反対にあうことも予想されていたので、断ることにした。

 

別記事に書くが、今回の活動で受けたのは全11社。うち、メーカーEとITのL以外は幸いなことに内定を頂いた。

コンサル（総合系・IT）：A、P、D、K、I

メーカー：H、D、E

IT：L、D

保険：S