Abnormally Distributed

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

説明変数のある時系列予測(状態空間モデル)

概要

時系列データは実務で取り扱う機会が多い。その際、予測対象の時系列データ(目的変数)と似たような変動をする他の時系列データ(説明変数)を手がかりに予測したい場合がある。

説明変数を用いて目的変数を予測する最も単純な方法として、線形回帰がある。ただ、線形回帰では時系列データの特徴を捉えきれないという問題があり、高い精度は期待できない。

この記事では時系列分析の1手法である状態空間モデルによる予測を行い、線形回帰に比べて何が優れているかを説明する。

データ

ここでは、下図のように変動するX, Yという2つの変数を考える。 Yが予測対象の変数、Xは説明変数である。 Xを利用して、Yをリアルタイムに予測したい。つまり、t時点のYを予測する際には、t-1時点までのYの値とt時点までのXの値が利用可能とする。

f:id:kibya:20200104224855p:plain

XとYの散布図は下図のようになる。

f:id:kibya:20200104225925p:plain

XとYにははっきりと相関があり、XはYの予測に役立ちそうである。

線形回帰

ここで単純な予測モデル、すなわち線形回帰することを考える。
これはt-1時点までのX, Yを用いて回帰直線を作成し(切片αと傾きβを求め)、t時点のXの値を代入することでYの予測値を算出する。

 {
\begin{eqnarray} 
y_t &=& \alpha+ \beta x_t + \varepsilon_t \\
\varepsilon_t &\sim& N(0, \sigma^2_\varepsilon) \\
\end{eqnarray} 
}

下図は、t=10までのデータ(グレー)を用いて、回帰直線を作成し、t=11のYを予測している。

f:id:kibya:20200104230756p:plain

全時点の予測結果をプロットすると下図の通り。

f:id:kibya:20200104231634p:plain

そこまで悪い予測ではないが、途中に連続して予測を大きく外してしまうケースが見られる。一期先予測誤差(MAE)は2.38程度となった。 これは線形回帰では、過去のデータ点をすべて同等に扱っており、直近のデータは類似しているというデータからはっきり読み取れる傾向を無視しているからである。
この直近のデータは類似している、という性質こそ時系列データの特徴であり、予測の際に考慮すべき情報なのだ。

状態空間モデル

時系列予測の1手法として、状態空間モデルがある。
特に今回扱うのは、一般的に「回帰成分を含むローカルレベルモデル」、と呼ばれるモデルで、式で表すと下記の通り。
線形回帰と違うのは、回帰直線の傾きと切片が、1つ前の時点の傾き、切片にノイズが加わって得られる、としている点である。

 {
\begin{eqnarray} 
y_t &=& \alpha_{t} + \beta_{t} x_t + \varepsilon_t \\
\alpha_{t} &=& \alpha_{t-1} + \eta_{\alpha,t} \\
\beta_{t} &=& \beta_{t-1} + \eta_{\beta,t} \\
\varepsilon_t &\sim& N(0, \sigma^2_\varepsilon) \\
\eta_{\alpha,t} &\sim& N(0, \sigma^2_{\eta_{\alpha}}) \\
\eta_{\beta,t} &\sim& N(0, \sigma^2_{\eta_{\beta}}) \\
\end{eqnarray} 
}

回帰直線の切片と傾きが「状態」と呼ばれる潜在変数である。各時点の状態は1時点前の状態によって決まる。
つまり予測の際は直近のデータの影響が大きくなる。

状態の推定にはカルマンフィルタと呼ばれる手法が用いられる。詳しくはこちらの記事を参照。 logics-of-blue.com

カルマンフィルタの要点としては、観測値が得られるたびに状態の推定値が補正され、補正された状態を用いて次の時点における状態と観測値の予測が行われる、ということになる。

状態空間モデルによる一期先予測の結果は下図の通り。

f:id:kibya:20200107075456p:plain

一期先予測誤差(MAE)は1.81程度となり、線形回帰よりはっきり改善している。

予測結果の比較

2つのモデルの違いを見るため、t=47, 48における予測結果を比べてみる。

まず、線形回帰の場合。この場合は過去の全ての観測値を同等に扱って予測を行う。結果、t=47, 48のいずれでも実際の観測値より大幅に小さい値を予測してしまっている。

f:id:kibya:20200107080906p:plainf:id:kibya:20200107080925p:plain

次に、状態空間モデルの場合。t=47では予測誤差が大きいのは同じだが、t=48の予測誤差は比較的小さくなっている。 これは、t=47の観測値を用いた状態の補正が行われているためである。つまり、直近の観測値は予測値よりだいぶ大きかった、という情報をモデルに与えることで、次の時点の予測値を大きくするように補正が働く

f:id:kibya:20200107081003p:plainf:id:kibya:20200107081013p:plain

アニメーションも載せておく。

f:id:kibya:20200111084459g:plainf:id:kibya:20200111084530g: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})}

下図は、3つの二変量正規分布からなる混合分布の例。それぞれの分布の平均を✖️で、共分散行列を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')

参考文献

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

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

はじめてのパターン認識

はじめてのパターン認識

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分布から得られたサンプルを精度行列とする2次元正規分布について、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. 正定値対称行列の制約があるため、全ての成分を少しずつ動かしたサンプルが有効となる確率はほぼ0となる。最初に制約のない空間に変換した上でサンプリングし、制約のある空間に戻した方が好ましい。
  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法まで)

「基礎からのベイズ統計学」で紹介された例を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次独立なものの数」である。
ベクトルの11次独立というのは、他のベクトルを定数倍したり足し合わせたりして表すことができないベクトルの組を意味する。


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による写像ではすべての点が1直線上に乗ってしまうので、元の2次元空間の情報は失われてしまう。
つまり逆写像は存在せず、逆行列も存在しないということになる。

このように、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転職活動まとめ①〜エージェント〜

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

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

 

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

 

エージェント

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

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

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

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

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

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

 

それぞれの寸評。

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

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

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

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

 

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

 

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

 

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

 

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

 

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

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

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

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

 

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

コンサル(総合系・IT):A、P、D、K、I

メーカー:H、D、E

IT:L、D

保険:S