開かない瓶の蓋をバケットクラッシャーで破砕する

機械学習・バイオインフォマティクス関連の気になった内容を記事にしていきます．

Mixtureモデルの周辺尤度をモンテカルロ積分で近似してみる

はじめに

この記事では、Mixtureモデルの周辺尤度をモンテカルロ積分で近似するアルゴリズムChib's methodを試します。 Mixtureモデルの周辺尤度は簡単には求まらない（数値的にも）という印象だったので、サンプリングを使って計算する手法がどこまでうまく行くのか検証しようと思いました。

経緯

Non-local prior というモデル選択のための事前分布の紹介qiitaを読んで、興味を持ちNon-local priorの論文も読みました。この論文はクラスタのダブりや空のクラスタが現れる事前確率を0にすることで、周辺尤度を使ったモデル選択（この場合はクラスタ数選択）の性能を上げるというものでした。この論文は面白かったのですが、周辺尤度の計算には疑問が残りました。見てみると周辺尤度の計算にはこの論文（Marin & Robert, 2008）の方法を使っており、この方法はこの論文内でChib's methodと呼ばれる方法をベースとしています。そこで、Chib's methodからMarin & Robertの方法からNon-local priorまですべて計算機実験でためそう！となったわけです。

そして、Chib's methodの検証までで飽きたのでそこまでを記事にすることにしました。

想定するモデル

想定するモデルは以下のようなものです。

  • Mixtureモデル
  • サンプルはi.i.d.
  • パラメタの事前分布はconjugate

Chib's method

まず、Chib's methodを説明します。ベイズの定理より、

p(X) = \frac{p(X|\theta)p(\theta)}{p(\theta|X)}

が成り立ちます。ここで、 Xはデータ、 \thetaはパラメタを指します。今回欲しいのは周辺尤度 p(X)なので、事後分布 p(\theta|X)さえ計算できればOKです。しかも左辺は \thetaに依存しないので、 \thetaとして任意の値 \theta^\starを代入した確率密度が計算できればよいです。簡単簡単。と、思いきや、肝心なことを忘れています。そうですクラスタを示す潜在変数です。 \thetaに含めて考えることもできますが、潜在変数とパラメタの同時分布の確率密度を得る必要が出てきて苦しいです。

そこで、潜在変数 Zを周辺化して上式を得ることを考えます。 p(X|\theta)と p(\theta)は簡単です。前者は Zの事前分布の割合で按分するだけで、後者はそもそも必要ないです。問題は p(\theta|X)です。簡単には求まりません。

解決法は単純です。MCMCによって Zのサンプルを得てそれを用いて近似して評価します。すなわち、

p(\theta|X) = \sum_Z p(\theta, Z|X) = \sum_Z p(\theta|Z, X) p(Z|X) \approx \frac{1}{S} \sum_s^S p(\theta | X, z^{(s)})

と近似します。ここで、 z^{(s)}はMCMCによる s番目の Zの事後サンプルを示し、 Sはサンプリング回数を示します。 これによって周辺尤度 p(X)が近似できます。この方法がChib's methodです。 また、Marin & Robertの方法ではMCMCにおけるlabel switchingの問題に対処するため \theta^\starを各クラスタについてパーミュテーションした値全通りで計算して平均を取っています。 今回は \theta^\starとして各クラスタで同じ値を使っています。いい値を選ばないと精度が悪くなるとのことですが、今回は気にしていません。

検証実験

ここからは検証実験の内容になります。

モデル

用いるモデルは超簡単なGaussian mixture modelです。具体的な生成過程は、

\mu_k \sim \mathrm{Normal}(\mathbf{0}, \mathbf{I})
z_i \sim \mathrm{Categorical}(\mathbf{1}/K)
x_i \sim \mathrm{Normal}(\mu_{z_i}, \mathbf{I})

です( \thetaが \muになってます)。厳密には記号の定義をするべきですが、しません。ハートで理解してください。

正解のための解析的な計算

ここからが問題です。Chib's methodを検証するためには正確な周辺尤度を計算しないといけません。しかし今回対象にできるモデルはmixtureモデルのみなのでそれができたら苦労しない...。なので、サンプルサイズを小さくし、 Zによる積分を厳密に計算できるデータで検証を行いました。 \sum_Zが計算できるなら正確な周辺尤度を出せるはず...！つまりは、上の周辺事後分布の

\sum_Z p(\mu|Z, X) p(Z|X)

の部分を計算します。 まず p(\mu|Z, X)です。こちらは、単峰分布の場合の事後分布の計算とほぼ変わりません。すなわち、

p(\mu_k|Z, X) \propto \left(\prod_{i: z_i = k} p(x_i | \mu_k)\right)p(\mu_k)

のように計算できます。右辺はクラスタkに属するデータについてだけ（クラスタの分け方すなわち Zは与えられているので可能）単峰の場合の事後分布を計算することになるので、 p(\mu_k)が共役事前分布であることを使って正規化し左辺を求められます。 今回の場合では、

p(\mu_k|Z, X) = \mathrm{Normal}\left(\frac{\sum_{i: z_i=k} x_i}{N_k+1}, \frac{\mathbf{I}}{N_k+1}\right)

となります。ただし、 N_kはクラスタkに属するデータ数です。 p(\mu|Z, X)を計算するコードは以下になります。

def calc_ln_alpha(x, z, mu_star, D, K):
    ln_alpha = 0
    for k in range(K):
        x_k = x[z==k]
        N_k = x_k.shape[0]
        if N_k == 0:
            ln_alpha += multivariate_normal.logpdf(
                mu_star, mean=np.zeros(D), cov=np.identity(D))
            continue
        post_mu_ex = x_k.sum(0) / (N_k + 1)
        post_mu_cov = np.identity(D) / (N_k + 1)
        ln_alpha += multivariate_normal.logpdf(
            mu_star, mean=post_mu_ex, cov=post_mu_cov)
    return ln_alpha

次に p(Z|X)です。ベイズの定理から、

p(Z|X) \propto p(X|Z)p(Z)

が成り立ちます。ここで、

p(X|Z)=\int p(X, \mu| Z) d\mu=(\prod_i p(x_i|\mu,z_i))(\prod_k p(\mu_k))=\prod_k \prod_{i: z_i=k} p(x_i|\mu_k) p(\mu_k)

となり、kに所属するxについてだけ注目したときの、周辺尤度p(X)を計算すればよいことがわかります。これは単峰の場合の周辺尤度なので計算できます。具体的には、単峰を考えて Zを省略するとして、

p(X) = \frac{p(X|\mu)p(\mu)}{p(\mu|X)}
です。これは、対数を取って、
\left(\sum_i \log \left( \mathrm{Normal}(x_i | \mu_k, \mathbf{I}) \right) \right) + \log \left( \mathrm{Normal}(\mu_k | \mathbf{0}, \mathbf{I}) \right) - \log \left( \mathrm{Normal}\left(\mu_k \middle| \frac{\sum_i x_i}{N_k+1}, \frac{\mathrm{I}}{N_k + 1}\right) \right)
で計算できます。最終的には、
\log p(X) = - \frac{DN}{2} \log 2 \pi -\frac{D}{2} \log (N_k + 1) -\frac{1}{2} \sum_i x_i^{\mathrm{T}} x_i - \frac{1}{2} \frac{1}{N_k +1} \left( \sum_i x_i \right)^{\mathrm{T}}\left( \sum_i x_i \right)
と計算できます（実は何度か計算ミスしましたが、多分これであってます）。 すべての Zの取り方（KN通り）についてこれと p(Z)をかけ合わせたものを計算して、logsumexpで正規化すれば p(Z|X)が計算できます。 p(Z|X)を計算するコードは以下になります。 
def calc_ln_betas(x, N, D, K):
    ln_betas = []
    for i in range(pow(K, N)):
        this_z = []
        cluster_str = None
        if K == 1:
            cluster_str = ''.zfill(N)
        else:
            cluster_str = np.base_repr(i, K).zfill(N)
        for s in cluster_str:
            this_z.append(int(s))
        this_z = np.array(this_z)
        ln_beta = 0
        for k in range(K):
            x_k = x[this_z==k]
            N_k = x_k.shape[0]
            if N_k == 0:
                continue
            ln_betak = - D * N_k * np.log(2 * np.pi) / 2
            ln_betak -= D * np.log(N_k + 1) / 2
            x_term = 0
            for i in range(N_k):
                x_term += x_k[i] @ x_k[i]
            x_term -= x_k.sum(0) @ x_k.sum(0) / (N_k+1)
            x_term /= -2
            ln_betak += x_term
            ln_beta += ln_betak
        ln_betas.append(ln_beta)
    ln_betas -= logsumexp(ln_betas)
    return ln_betas
最後に、掛け合わせてKN通りの和を取って、 \sum_Z p(\mu|Z, X) p(Z|X)を得ます。 コードは以下です。 \mu^{\star}としていくつかの値を計算して同じ値になったので実装に問題はない...と思ってます。 
def calc_true_log_posterior(x, mu_star, N, D, K):
    ln_etazs = []
    ln_betas = calc_ln_betas(x, N, D, K)
    for i in range(pow(K, N)):
        this_z = []
        cluster_str = None
        if K == 1:
            cluster_str = ''.zfill(N)
        else:
            cluster_str = np.base_repr(i, K).zfill(N)
        for s in cluster_str:
            this_z.append(int(s))
        this_z = np.array(this_z)
        ln_etaz = calc_ln_alpha(x, this_z, mu_star, D, K) + ln_betas[i]
        ln_etazs.append(ln_etaz)
    return logsumexp(ln_etazs)

Chib's methodのためのギブスサンプリング実装

こちらは普通にギブスサンプリングです。初期値から、 \muを使った尤度と事前確率の値を正規化した分布から Zをサンプリング、 Zをgivenとして共役事前分布を利用して計算した事後分布から\muをサンプリングします。まあ定式化はしなくていいでしょう。ギブスサンプリングとそれを使って対数事後確率を近似するコードは以下です。 

def GMM_gibbs_sample(x, K, T):
    N, D = x.shape
    mu_sample = np.random.normal(0, 1, size=(K, D))
    z_samples = []
    mu_samples = []
    mu_samples.append(mu_sample)
    for t in tqdm(range(T)):
        z_sample = np.zeros((N), dtype=int)
        for i in range(N):
            ln_posterior = np.zeros(K)
            for k in range(K):
                ll = multivariate_normal.logpdf(x[i],
                                                mean=mu_sample[k],
                                                cov=np.identity(D))
                ln_prior = np.log(1 / K)
                ln_posterior[k] = ll + ln_prior
            ln_posterior -= logsumexp(ln_posterior)
            z_sample[i] = np.random.choice(range(K), p=np.exp(ln_posterior))
        z_samples.append(z_sample)
        for k in range(K):
            x_k = x[z_sample==k]
            N_k = x_k.shape[0]
            mu_sample[k] = np.random.multivariate_normal(
                x_k.sum(0) / (N_k + 1), np.identity(D) / (N_k + 1))
        mu_samples.append(mu_sample.copy())
    return np.array(z_samples), np.array(mu_samples)



def calc_GS_approx_log_posterior(x, mu_star, D, K):
    T = 1000
    z_samples, _ = GMM_gibbs_sample(x, K, T)
    n_burn_in = 100
    ln_alphas = []
    for z_sample in z_samples[n_burn_in:]:
        ln_alphas.append(calc_ln_alpha(x, z_sample, mu_star, D, K))
    lpos_approx = logsumexp(ln_alphas) - np.log(T-n_burn_in)
    return lpos_approx

カーネル密度推定

加えて、共役事前分布を利用しないケースも考えてみます。共役事前分布でないと p(\mu|X)を解析的に求められないため、MCMCの \muのサンプルからカーネル密度推定で確率密度関数の値を近似することにします。 この場合のMCMCはギブスサンプリングを使ってしまうとズルなので（上記のギブスサンプリングは共役事前分布を利用している）、Stanでサンプリングすることにします。 Stanコード（を格納する文字列変数を定義するpythonコード）は以下です。 

stan_code = """
data {
    int<lower=0> N;
    int<lower=0> D;
    int<lower=0> K;
    vector[D] X[N];
}
parameters {
    vector[D] mu[K];
}
transformed parameters{
    vector[K] ps[N];
    for (i in 1:N){
        for (k in 1:K){
            ps[i, k] = normal_lpdf(X[i] | mu[k], 1) + log(1.0/K);
        }
    }
}
model {
    for (k in 1:K){
        mu[k] ~ normal(0, 1);
    }
    for (i in 1:N){
        target += log_sum_exp(ps[i]);
    }
}
このStanコードを使って近似した p(\mu^{\star})を得るコードは以下です。 
def calc_KDE_approx_log_posterior(x, mu_star, N, D, K):
    stan_data = {"N": N,
                 "D": D,
                 "K": K,
                 "X": x,
                }
    random_seed = np.random.randint(0, 10e6)
    posterior = stan.build(stan_code, data=stan_data, random_seed=random_seed)
    fit = posterior.sample(num_chains=4, num_samples=10000)
    # クラスタを区別しないmuの事後サンプル
    mu_nuts = fit["mu"]
    mu_post_samples = mu_nuts.transpose((2, 0, 1)).reshape(-1, D)
    kde = KernelDensity(kernel='gaussian', bandwidth='silverman').fit(mu_post_samples)
    log_density = kde.score(mu_star.reshape(1, -1))
    return log_density

対数周辺尤度の計算

これで、三者の方法による事後分布の確率密度関数の値の計算ができました。あとは対数周辺尤度の計算です。 事前分布と尤度は共通です。以下のコードで計算できます（ \mu^\starとして各クラスタで同じ値を用いる前提のコードです）。 

def calc_log_prior(mu_star, D, K):
    lpri = 0
    for k in range(K):
        lpri += multivariate_normal.logpdf(mu_star, np.zeros(D), np.identity(D))
    return lpri



def calc_log_likelihood(x, mu_star, D):
    lls = multivariate_normal.logpdf(x, mu_star, np.identity(D))
    return lls.sum()
対数周辺尤度は、事後分布を計算する関数を引数として、以下のように書きました。 
def calc_lml(x, mu_star, N, D, K, calc_log_posterior):
    ll = calc_log_likelihood(x, mu_star, D)
    lpri = calc_log_prior(mu_star, D, K)
    lpos = calc_log_posterior(x, mu_star, N, D, K)
    return ll + lpri - lpos

実験結果

では、いよいよ三者の比較を行います。確認です。 True: 正解周辺尤度を厳密に計算。共役事前分布が必要で、Nは30が限度。 GS: Chib's methodを使った近似計算。共役事前分布が必要。 KDE: KDEを使ったChib's method。どんなモデルでも使える。KDEのためのハイパーパラメタは必要。

データは以下のコードで生成しました。まあ正解パラメタを当てたいわけではないのでサンプルサイズと次元さえ決められればどんな作り方でも良いです。 

def generate_data_and_param(N, D, K):
    mu = np.random.multivariate_normal(np.zeros(D), np.identity(D), size=K)
    z = np.random.choice(range(K), size=N)
    x = []
    for i in range(N):
        x.append(np.random.multivariate_normal(mu[z[i]], np.identity(D)))
    x = np.array(x)
    return x, z, mu
サンプルサイズは3, 5, 10、サンプルの次元は1, 3, 5、クラスタ数は1, 2の全ての組み合わせで作成したデータで実験を行いました。 比較実験は以下のコードで行いました。各条件で生成したデータ一つに対して、各近似手法で100回ずつ計算を行って近似のブレも調べました。 
T = 100
settings = [(N, D, K) for N in [3, 5, 10] for D in [1, 3, 5] for K in [1, 2]]
for N, D, K in settings:
    output_filename = f'N{N}D{D}K{K}.pickle'
    print(output_filename)
    if os.path.exists(output_filename):
        continue
    x, z, mu = generate_data_and_param(N, D, K)
    mu_star = x.mean(0)
    true_lml = calc_lml(x, mu_star, N, D, K,
                        calc_true_log_posterior)
    GS_lmls = []
    KDE_lmls = []
    for t in tqdm(range(T)):
        GS_lml = calc_lml(x, mu_star, N, D, K,
                          calc_GS_approx_log_posterior)
        KDE_lml = calc_lml(x, mu_star, N, D, K,
                           calc_KDE_approx_log_posterior)
        GS_lmls.append(GS_lml)
        KDE_lmls.append(KDE_lml)
    result = [true_lml, GS_lmls, KDE_lmls]
    with open(output_filename, mode='wb') as f:
        pickle.dump(result, f)
そして可視化はこちらのコードで行いました。 
fig, axes = plt.subplots(3, 6, figsize=(30, 15))
fontsize = 20
settings = [(N, D, K) for N in [3, 5, 10] for D in [1, 3, 5] for K in [1, 2]]
for t, (N, D, K) in enumerate(settings):
    i = t // 6
    j = t % 6
    ax = axes[i][j]
    condition = f'N{N}D{D}K{K}'
    output_filename = f'{condition}.pickle'
    with open(output_filename, mode='rb') as f:
        true_lml, GS_lmls, KDE_lmls = pickle.load(f)
    approx_lml_df = pd.DataFrame([GS_lmls, KDE_lmls], index=['GS', 'KDE']).T
    sns.violinplot(x='variable', y='value', data=pd.melt(approx_lml_df), scale='width', ax=ax)
    plt.sca(ax)
    plt.hlines(true_lml, -0.5, 1.5, color='blue')
    plt.title(condition, fontsize=fontsize)
    plt.xticks(fontsize=fontsize)
    plt.yticks(fontsize=fontsize)
    plt.xlabel('Method', fontsize=fontsize)
    plt.ylabel('Log marginal likelihood', fontsize=fontsize)
plt.plot()
plt.tight_layout()
結果は以下になります。
各実験設定における対数周辺尤度。横軸は近似手法、縦軸は対数周辺尤度の値を示す。青色の直線は真の値を示す。それぞれの図はそれぞれの実験設定に対応する。例えば、「N10D3K2」はサンプルサイズ10、サンプルの次元3、クラスタ数2のデータとモデルでの結果を示す。

まず、肝心のGS（Chib's method）ではそこそこ当たっていそうな感じです。あと、サンプルサイズや次元を大きくすると急激に精度が悪化するようなこともなさそうです（K=2までしか調べてないですし、各条件でデータセットも一つしか試していないので雰囲気です）。KDEだとバイアスが大きそうですね。やり方をちゃんと吟味すればうまく行くかもしれません。 周辺尤度の近似精度としてはそんな感じです。

やり残したこと

近似精度の解釈

ただ、周辺尤度を計算したいときというのはモデル選択を行いたいときが多いはずで、他のモデルでの周辺尤度との差も考慮して誤差を評価したいところです。例えば、近似値のほとんどが真の値の±0.1くらいに収まったとしても、比較したいモデルの真の値に0.001くらいの差しかなければ近似手法の有用性は低そうです（ノイズに差が埋もれてしまう）。

渡辺ベイズ理論における周辺尤度の近似

渡辺ベイズ（ベイズ統計の理論と方法、5.1.4 自由エネルギーの近似）では周辺尤度は逆温度 \betaを少しずつ変えた複数回のMCMCにより計算できる一方で、ナイーブなモンテカルロ積分ではうまくいかないと書かれています。 WBICも含めて、渡辺ベイズ理論の手法とも比較したいところでした。畜生、根気さえあれば...！

\theta^\starの吟味

元論文では最尤推定などで事後分布の大きいところを占める \thetaの値を \theta^\starとして設定することが推奨されていました。今回は各クラスタで同じ値を使ってしまったので、理論の上では同じことだとしても、性能に差が出る可能性があります。

高次元データでの性能

今回は低次元データのみでの検証でした。高次元データではどうなるのかわかりません...ので、軽くやってみたところ... ん...？あれ...？若干分布が怪しいですね。中央値を見ればそこそこ近そうですが... さらなる調査が望まれる。

まとめ

いかがでしたか？ タイトルの内容を試して、理論値（と言っていいのか？正解の値？）と比較しました。 やり残したことはありますが、個人的には満足です。 では。

ガウス過程の導関数を求める

はじめに

ガウス過程回帰による予測値の導関数をGPyのpredictive_gradientsメソッドによって簡単に得られることがわかったので記しておこうと思います．

ガウス過程の導関数

まずガウス過程によって推定された関数の導関数が求められることを確認します．こちらのページを参考にしました．

まず，ガウス過程回帰における予測値(事後期待値)は

E[y^* | x^*, D] = k_{*}^{T} K^{-1} y

で計算できます．ここで，

y^*は新規観測値の目的変数， x^*は新規観測値の説明変数， Dはデータ， k_{*}はデータの説明変数と新規観測値の目的変数のカーネル値(データは複数あるのでベクトル)， Kはデータの説明変数のカーネル行列， yはデータの目的変数を示します．そしてこれらを x^*で微分すると，

\frac{d}{dx^*}E[y^*|x^*, D] = \frac{dk_{*}^{T}}{dx^*} K^{-1} y

これにより，ガウス過程回帰の事後期待値の導関数はカーネル値のベクトルの導関数を用いて求められることが分かります．

GPyによる導関数値の取得

まずはデータを作って可視化します．サインカーブにでもしてみましょう．

import GPy
import numpy as np
import matplotlib.pyplot as plt

N=100
np.random.seed(0)
x = np.random.normal(0, 1, size=N)
y = np.sin(x) + np.random.normal(0, 0.1, size=N)

kernel = GPy.kern.RBF(1)

model = GPy.models.GPRegression(x.reshape(-1, 1), y.reshape(-1, 1), kernel=kernel)
model.optimize()
model.plot()
plt.show()

f:id:oaietnbvgq30:20200801221327p:plain

GPyではガウス過程回帰の導関数がpredictive_gradientsメソッドで簡単に求められます．

このデータの導関数値を数値解析的に求めます．

x_new = np.array([[0.0]])
epsilon = 1.0e-5
x0_minus_0 = model.predict_quantiles(x_new - epsilon, quantiles=(5, 50, 95))[1][0][0]
x0_plus_0 = model.predict_quantiles(x_new + epsilon, quantiles=(5, 50, 95))[1][0][0]
print((x0_plus_0 - x0_minus_0) / (2 * epsilon))
# 0.9633798891961475

ここでは小さい値だけ前後にずらした2点の傾きを求めています． 今度はpredictive_gradientsメソッドを使ってみます．

print(model.predictive_gradients(x_new)[0][0][0][0]) 
# 0.9633798901678498

ちゃんと求められていますね(数値解析的に求めたものを答えとするのってどうなんだろう…)．GPyのgithubからも分かるように予測値の分散の導関数(導関数の分散ではない)も一つ目の添字を1にすることで得られます．

まとめ

今回はガウス過程回帰の予測値の導関数を求められることの確認とGPyライブラリによる計算方法の紹介でした．

pystanでglassoを実装してみる

はじめに

変数間のinteractionを発見するのに使われるglassoの拡張モデルを作る準備としてstanでglassoを実装してみます．ただglassoを使うだけなら

sklearn.covariance.graphical_lasso — scikit-learn 0.23.1 documentation

を使うのが一番手っ取り早いです．

stanファイルの実装

stanでパラメタ推定を行う場合はモデルを記述するstanファイルを用意する必要があります． 詳細は省きますが，glassoに置けるL1正則化項は正規分布の精度行列パラメタの事前分布と見なせます． この事前分布は精度行列の各要素を対角成分と非対角成分に分けて，それぞれ独立同分布とみなすものです． ただし，この事前分布は精度行列としての条件(逆行列＝分散共分散行列が半正定値行列)を満たすことが保証されないのでstanでの実装においてはcov_matrix型を定義することで対処します． このcov_matrix型で宣言されたパラメタは分散共分散行列としての条件を満たすことが保証されます． そして実装したのが以下になります．

// glasso.stan
data {
    int<lower=0> n;
    int<lower=0> p;
    vector[p] x[n];
    real<lower=0> lambda;
    // real<lower=0> prior_parameter_for_lambda;
}

transformed data {
}

parameters {
    cov_matrix[p] invsigma;
    vector[p] mu;
    // real<lower=0> lambda;
}

transformed parameters {
    real<lower=0> invlambda;
    real<lower=0> lambda_half;
    invlambda = 1/lambda;
    lambda_half = lambda/2;
}

model {
    for (i in 1:n) {
        x[i] ~ multi_normal_prec(mu, invsigma);
    }
    for (j1 in 1:p) {
        for (j2 in 1:j1) {
            if (j1==j2) {
                invsigma[j1,j2] ~ exponential(lambda_half);
            } else {
                invsigma[j1,j2] ~ double_exponential(0,invlambda);
            }
        }
    }
    //Prior for lambda
    // lambda ~ exponential(prior_parameter_for_lambda);
}

generated quantities {
}

こちらの論文&コードを参考にしました(なんでバイオインフォやねんというのはたまたま思い浮かんだのがこの論文でして...)

人工データによる実験

生成過程から人工データを作成して真のパラメタを当てられるか検証します． まず，importとデータの次元を決めます．

import pandas as pd
import numpy as np
import sys
import seaborn as sns
sns.set()
import matplotlib.pyplot as plt
import pystan

N = 1000
p = 5

精度行列を作ります．ただし，精度行列をglassoの生成過程通りに作るとスパースにならないので，categorical分布で0か否かを生成してから値を生成します．

prec_matrix = np.zeros([p,p])

for i in range(p):
    for j in range(i+1):
        if i == j:
            prec_matrix[i, j] = np.random.gamma(100,1)
        else: 
            prec_matrix[i, j] = np.random.choice(2, p=[0.8, 0.2])
            prec_matrix[i, j] = prec_matrix[i, j] * np.random.normal(0,100)
            prec_matrix[j, i] = prec_matrix[i, j]

固有値の最小値を出力し，負でないことを確認します(負であれば再生成します)．

np.linalg.eigvals(np.linalg.inv(prec_matrix)).min()
# 0.005406751950232075

精度行列を可視化します．

sns.heatmap(prec_matrix, cmap='coolwarm', center=0)

f:id:oaietnbvgq30:20200528100608p:plain

データを生成します．平均は0にします．

x = []
for i in range(N):
    x.append(np.random.multivariate_normal(np.zeros(5), np.linalg.inv(prec_matrix)))
x = pd.DataFrame(x)

罰則項の係数は1にしてみます．

stan_data = {'n': N, 'p': p, 'x': x, 'lambda': 1}
fit = pystan.stan(file='./glasso.stan', data=stan_data, iter=1000, chains=1)
print('Sampling finished.')

一応MAP推定をしたいので，サンプル列からposteriorが最大のindexを取ってきて，該当するサンプルを可視化します．

la  = fit.extract(permuted=True)
names = fit.model_pars
map_index = pd.Series(fit.get_logposterior()[0]).idxmax()
la  = fit.extract(permuted=True)
invsigma_sample = pd.DataFrame(la[names[0]].reshape(-1, la[names[0]][0].size))
# burn_in期間の500を引いたindexを使用
estimated_prec_matrix = pd.DataFrame(np.array(invsigma_sample.iloc[map_index - 500]).reshape(p, p))
sns.heatmap(estimated_prec_matrix, cmap='coolwarm', center=0, annot=True)

f:id:oaietnbvgq30:20200528103020p:plain

うまく推定できていますね． ただ，posteriorが最大になるサンプルを取ってきているだけなので0に見える部分は0ではなかったりします．

おわりに

glassoはposteriorを最大化した場合にスパースになりやすいという手法なのでMCMCで行うのは邪道かもしれませんね．精度行列の非対角成分にラプラス分布の代わりに正規分布をおいてもたいして変わらないと思います．しかし，MCMCで近似された事後分布の形はおおよそ正しい訳ですから，考案したモデルをスクラッチから実装する前の確認としては十分有用なのでは無いかと思います．

特異値分解(singular value decomposition)でベクトルを直交化する

はじめに

特異値分解を使った直交化をする機会があったので書き留めておきます．数式で説明するか文章で説明するかはフィーリングです．

モチベーション

行列Mが与えられたとします．Mの列ベクトルを，変換後の空間のスケール(厳密な用語が分かりません...)を変えないように直交したベクトルに変換することを行いたいという話です．なぜかについては説明しません(面倒ですので・・・)．

説明

間違ってたらご一報頂けると幸いです．

特異値分解

行列分解の一種．行列MをM=U \Sigma Vのように分解します．ここで， U, Vは直交行列， \Sigmaは対角行列です．

直交行列の性質

行列が直交行列である，とは以下のことと同値です．

  1. 各行・列ベクトルが直交しておりノルムが1
  2. 逆行列は転置した行列
  3. 変換が等長写像(ノルムを変えない写像)になる

今回の直交化においてこの性質が大事になってきます．

特異値分解(SVD)による直交化

SVDによる直交化は M=U \Sigma V を利用して M'=MV ^T=U \Sigma のように直交化します．ここで，性質2より，Vの逆行列をかける代わりに転置行列をかけています．

特異値分解(SVD)による直交化の性質

SVDによる直交化では，被変換行列Mによる変換後の単位円を変えずに各ベクトルを直交にすることができます．以下，説明です．全てのノルムが1のベクトルを考えたとき，座標とみなしたベクトルの集合によって円ができます．これを単位円と言います．Mによってその各座標を変換すると単位円の形が変わります(楕円もしくは円に)．SVDによる直交化では，この楕円の形を変えないベクトルに直交化できるのです．もちろん，各ベクトルについて見てみれば直交化する前と後では値が変わってしまいます．ただ，ノルムが一定のベクトルによってできる，集合レベルでは変わらないようになっているということですね． この話の直感的な図がWikipediaにあります．

f:id:oaietnbvgq30:20191209170006p:plain
赤と黄色のベクトルは変換対象のベクトルです．Mによる変換(左上⇨右上)とU \Sigmaによる変換(左下⇨右上)では単位円の形を変えないことが分かります．
en.wikipedia.org

この直交化を行列Mに対して行うと，

MV ^T=U \Sigma

となり，SVDによって得られるUと\Sigmaで変換後の行列を得られることになります(もちろん元の行列にV ^Tをかけてもよい)． 更にU \Sigmaを見てみると直交行列と対角行列の行列積になっているので，性質1より，変換後の行列は列ベクトルが直交していることが分かります(行ベクトルは直交するとは限らない)．

要するに，V ^Tによる変換をしてあげると，変換後の空間のスケール(厳密な用語がry)を変えずに列ベクトルを直交化させることができる，というわけですね．

実験

実際にSVDを使って得られるベクトルを眺めてみます．

実験条件

pythonで動かしていきます． まずimportします．numpyと可視化関連です(個人的な好みでseabornを使ってますがimportしなくてもこの後のコードは実行できます，figureの見た目が良くなるだけです)．

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

適当な行列を作ります．

M = np.array([[1,4],[3,4]])

SVDしてU \Sigmaを計算することでMを直交化します(numpyの関数だと \Sigmaは対角成分のみのベクトルで返ってくるためnp.diagで行列にしている)．

U, sigma, V = np.linalg.svd(M, full_matrices=True)
M_prime = np.dot(U, np.diag(sigma))

一応直交しているか確認．

np.dot(M_prime.T, M_prime)
# array([[  4.04164878e+01,  -1.21399145e-17],
# [ -1.21399145e-17,   1.58351216e+00]])

大丈夫そう．最後に可視化します．

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.set_xlim([-6,5])
ax.set_ylim([-6,5])
ax.plot([0,M.transpose()[0][0]], [0,M.transpose()[0][1]], color='b')
ax.plot([0,M.transpose()[1][0]], [0,M.transpose()[1][1]], color='b')
ax.plot([0,M_prime.transpose()[0][0]], [0,M_prime.transpose()[0][1]], color='r')
ax.plot([0,M_prime.transpose()[1][0]], [0,M_prime.transpose()[1][1]], color='r')
plt.show()

f:id:oaietnbvgq30:20191205012324p:plain

直交化後のM'の列ベクトルが描く単位円が変わってなさそうなことが分かります(適当)．