「Pythonではじめるベイズ機械学習入門」を読む　～第4章潜在変数モデル②(ベイジアン主成分分析)～

はじめに

Pythonによるベイズモデルの実装をきちんと学ぼうと思い，森賀新・木田悠歩・須山敦志著「Pythonではじめるベイズ機械学習入門」を読むことにした。

本記事は，第4章「潜在変数モデル」のうちベイジアン主成分分析に関する読書メモである。

本書の紹介ページ

www.kspub.co.jp

4.2 行列分解モデル

行列分解モデル(Matrix Factorization Model)は，観測値を表す行列を，より次元の小さい行列の積に分解することである。
ノイズを除外した本質的なデータの特徴を得られることや，欠損値を補完できることがその利点である。

4.2.1 モデル概要 : 主成分分析

主成分分析(Principal Component Analysis)は，多次元データの座標軸を，データの分散が最大になる方向に変換する手法である。
高次元データの次元圧縮などの目的で利用される。

4.2.2 モデル概要 : ベイジアン主成分分析

主成分分析は，「高次元のデータを低次元に圧縮する」という発想であったが，ベイジアン主成分分析(Bayesian PCA)はこれとは逆に， $K$ 次元の潜在変数 $z$ の線形変換にノイズが句分かり， $D$ 次元の観測変数 $x$ が得られると考える(ここで， $K \lt D$ である)。
すなわち，観測されるデータは高次元([tex : D]次元)であるが，実質的な情報は $K$ 次元の潜在変数により表現できると考える。

$N$ 個のデータ点からなる観測変数の集合を $X = \{ x_1, \cdots, x_n \}$ ，対応する潜在変数の集合を $Z = \{ z_1, \cdots, z_N \}$ とする。ただし， $x_n \in \mathbb{R}^D$ ， $z_n \in \mathbb{R}^K$ とする。

潜在変数 $z_n$ は平均が0，分散共分散行列が単位行列であるガウス分布に従うとする。

$\begin{align} z_n \sim \mathcal{N}(0, I) \\ \end{align}$

次に，潜在変数 $z_n$ が得られた下での観測変数 $x_n$ について，

$\begin{align} P(x_n | z_n) &= \mathcal{N}(x_n | W z_n + \mu, \sigma^2 I) \\ \\ W &\in \mathrm{R}^{D \times K} \\ \\ \mu &\in \mathrm{R}^D \\ \end{align}$

と表す。

ベイジアン主成分分析ではさらに， $W, \mu, \sigma$ の事前分布を設定する。ただしここでは簡単のため， $\mu = 0$ に固定し， $W, \sigma$ の事前分布のみを設定する。

本書のP152では， $W$ の列ベクトルを $w_d$ ，事前分布 $p(W)$ を，

$\begin{align} p(W) = \prod_{d=1}^{D} \mathcal{N} (w_d | 0, I) \\ \end{align}$

と説明しているが， $x_n, W, z_n$ の関係は下図のようになり， $w_d \in \mathbb{R}^D, d=1,\cdots,K$ となるはずなので，実際は，事前分布 $p(W)$ は

$\begin{align} p(W) = \prod_{d=1}^{K} \mathcal{N} (w_d | 0, I) \\ \end{align}$

となる。

$\sigma$ は正の値を取るため，事前分布は半コーシー分布となる。

モデル式をまとめると，以下のようになる。

$\begin{align} z_n &\sim \mathcal{N}(0, I) \\ \\ p(W) &= \prod_{d=1}^{K} \mathcal{N} (w_d | 0, I) \\ \\ p(\sigma) &= \mathrm{HalfCauchy}(\sigma | \beta) \\ \\ P(x_n | z_n) &= \mathcal{N}(x_n | W z_n + \mu, \sigma^2 I) \\ \end{align}$

4.2.3 実装

サンプルコードを動かしながら，挙動を確認した。
github.com

データの準備

以下のような観測データを対象とする。

モデルの定義

4.2.2節で示したようなモデルを定義する。
ただし，観測変数の次元は $D=2$ ，潜在変数の次元は $K=1$ ，サンプル数は $N=500$ ，標準偏差の真値は $\sigma=0.8$ としている。
また変換行列 $W$ は2行・1列の行列となる。

推定対象は，変換行列 $W$ の列ベクトル $w_d \in \mathbb{R}^1, d=1,2$ の値と，標準偏差 $\sigma$ である。

MCMCによる事後分布からのサンプリング

MCMCによるサンプリングを行なう。2024年12月現在，PyMC3ではなくPyMCを使うので，ライブラリのインポート設定は以下のようになる。

#import pymc3 as pm
import pymc as pm

サンプルの分布とトレースプロットは下図の通りである。

主成分分析における固有値・固有ベクトルの計算において，符号が反転することがある。トレースプロットを確認しても，チェーンによっては値の正負が入れ替わることがあるので，MCMCにおいても同様の現象が確認できた。

予測分布

サンプルコードはそのままでは動かなかったので，以下のように修正した。

サンプルコード(修正前)

with model:
    x_post = pm.sample_posterior_predictive(trace, 1, random_seed=42)

サンプルコード(修正後)

with model:
    x_post = pm.sample_posterior_predictive(trace, var_names=['x'], random_seed=42)

修正後のコードでは，取り出したい変数を”var_names”引数によって指定する必要がある。

予測分布の生成のコード変更により，描画用コードも変わってくる。

サンプルコード(修正前)

plt.scatter(x_train[0, :], x_train[1, :], color='blue', s=30, alpha=0.3, label='observed')
plt.scatter(x_post['x'][0, 0, :], x_post['x'][0, 1, :], color='red', marker='x', s=30, alpha=0.3, label='simulated')
plt.legend()
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')

サンプルコード(修正後)

plt.scatter(x_train[0, :], x_train[1, :], color='blue', s=30, alpha=0.3, label='observed')
plt.scatter(x_post.posterior_predictive['x'].sel(chain=[0], draw=[0], x_dim_2=[0]), 
            x_post.posterior_predictive['x'].sel(chain=[0], draw=[0], x_dim_2=[1]), 
            color='red', marker='x', s=30, alpha=0.3, label='simulated') 
plt.legend()
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')