東京大学・株式会社Nospareの菅澤です．
今回から多次元データの従属構造に対するベイズ分析の方法を2回に分けて紹介します．前半となる今回は，順位尤度を用いた事後分布とMCMCアルゴリズムについて紹介します．

従属構造のモデリングとコピュラ

多次元データが与えられたとき，変数間の関係性（従属構造）を把握することはデータの構造を理解する上で非常に重要です．多次元データの確率的な構造を考えるためには，データが従う多次元分布について考える必要があります．

多次元分布の分解

多次元分布は，ある条件のもとで各変数ごとの周辺分布と従属性を表す構造に一意に分解できることがスクラー（Sklar）の定理から知られています．これは多次元分布を考える際に，各変数ごとの確率的なバラつきと，従属構造は分けて考えられることを意味しています．このとき，従属構造の部分は一般的に コピュラ（copula） と呼ばれています．

以下では，各変数の周辺分布には興味がなく，従属構造（コピュラ）のみに興味がある状況を考え，コピュラのみを直接ベイズ推定することを考えます．

また，現実の応用例では連続値や離散値（カウント値や二値）が混ざった多次元データを扱うことが多々ありますが，そのような連続離散混合の多次元データから従属構造を推定する一般的な設定を考えます．

順位尤度 (rank likelihood)

$p$次元の連続離散混合データ$y_i=(y_{i1},\ldots,y_{ip}) \ (i=1,\ldots,n)$ が観測されているとします．このとき，$y_i$の同時分布をモデル化するため，潜在変数$z_i=(z_{i1},\ldots,z_{ip})\sim N(0, {\bf \Psi})$を導入し，$y_i$の各要素を$y_{ik}=G_k(z_{ik})$の形で表現します．ここで，${\bf \Psi}$は潜在変数の相関行列で，$G_k$は適当な非減少関数です．例えば，$y_{ik}$が3つの値をとる順序データの場合，$G_k(x)=1 \ (x\leq g_1)$, $G_k(x)=2 \ (g_1<x\leq g_2)$および$G_k(x)=3 \ (x\geq g_2)$となり，未知パラメータ$g_1, g_2$が含まれた関数になります．

上記の設定では，$y_i$の従属構造を推定することは$z_i$の相関行列${\bf \Psi}$を推定することに相当し，$G_k(\cdot)$は$y_i$の各変数の周辺分布を決める未知の関数です．したがって，従属構造${\bf \Psi}$のみを推定したい状況では，$G_k(\cdot)$は局外パラメータ（興味のないパラメータ）となります．

順位尤度の考えを用いることで，${\bf \Psi}$の尤度を直接的に与えることができ，局外パラメータを除いた${\bf \Psi}$の事後分布を直接的に与えることができます．重要な事実としては，$G_k(\cdot)$は非減少関数のため，$(y_{1k},\ldots,y_{nk})$の順序と$(z_{1k},\ldots,z_{nk})$の順序が同じになる点です．これは，$(y_{1k},\ldots,y_{nk})$の大小関係によって潜在変数$(z_{1k},\ldots,z_{nk})\in \mathbb{R}^n$の存在領域に制約がかかることを意味しています．

データと潜在変数の$n\times p$行列をそれぞれ${\bf Y}=(y_1,\ldots,y_n)$および${\bf Z}=(z_1,\ldots,z_n)$と定義しておきます．データ${\bf Y}$が得られたもとで，潜在変数${\bf Z}$が存在する集合は

R({\bf Y})=\{{\bf Z} \ |  \ z_{ik}< z_{ij} \ {\rm if}\ y_{ik}< y_{ij}, \ \ i=1,\ldots,n, \ \ \ j,k=1,\ldots,p \}

となります．また，この事象の確率は

P({\bf Z}\in R({\bf Y})\ |\ {\bf \Psi})=\int_{R({\bf Y})} \prod_{i=1}^n\phi(z_i; 0, {\bf \Psi}) dz_i

で与えられます．これが相関行列${\bf \Psi}$に対する順位尤度(rank likelihood)です．

${\bf \Psi}$に対する事前分布を$\pi({\bf \Psi})$とすると，事後分布は

\pi({\bf \Psi}\ |\ {\bf Y})\propto \pi({\bf \Psi})P({\bf Z}\in R({\bf Y})\ | \ {\bf \Psi})

で与えられます．このような順位尤度によるベイズ推定法はHoff (2007)で提案されています．

ギブスサンプラー

${\bf \Psi}$の事後分布は潜在変数${\bf Z}$に関する積分が含まれており直接的に計算するのは困難です．そこで，潜在変数${\bf Z}$からも同時にサンプリングをすることで，簡単にギブスサンプリングを行うことを考えます．

${\bf \Psi}$が所与のもとで，$z_i$は領域$R({\bf Y})$上に制限された多変量正規分布に従うため，比較的簡単にサンプル生成を行うことができます．(詳細は後述します．)

一方で，潜在変数${\bf Z}$が所与のもとでの${\bf \Psi}$の完全条件付き分布は$\pi({\bf \Psi})\prod_{i=1}^n\phi(z_i; 0, {\bf \Psi})$に比例します．これは相関行列${\bf \Psi}$をもつ多変量正規分布のデータから${\bf \Psi}$の事後分布を計算することに相当しますが，一般に相関行列$\Psi$に対する単純な共役事前分布のクラスは存在しないため，簡単そうに見えて意外と難しい問題に直面してしまいます．

もちろん，メトロポリスヘイスティングスアルゴリズムを使うことで，サンプリングのアルゴリズムを行うこと自体は可能ですが，非効率なサンプリングになってしまったり，相関行列の制約を満たしながらサンプル生成を行う必要があるため，テクニカルな作業が必要になってしまいます．

今回の問題に関しては，パラメータ拡大と呼ばれるテクニックを用いることで，簡便なギブスサンプラーを構成することができます．

パラメータ拡大

相関行列${\bf \Psi}$を直接対象とせずに，共分散行列${\bf \Sigma}$を経由した事前分布の特定化・サンプリングアルゴリズムを考えます．まず，$z_i$に対するモデルを$z_i\sim N(0, {\bf \Sigma})$と書き換えます．共分散行列${\bf \Sigma}$に対しては自然な共役事前分布として逆ウィシャート分布を用いることができます．

順位尤度$P({\bf Z}\in R({\bf Y})\ |\ {\bf \Psi})$は$z_{ij}$の大小関係のみに依存しているため，${\bf \Sigma}$の対角成分（各変数のスケール）を識別することができません．その意味で，$z_i\sim N(0, \Sigma)$は不良な定式化です．しかし，${\bf \Sigma}$を${\bf \Psi}=h({\bf \Sigma})\equiv \{\sigma_{ij}/\sqrt{\sigma_i^2\sigma_j^2}\}$と変換することで${\bf \Psi}$は一意に定まり，${\bf \Psi}$は順位尤度から識別可能です．したがって，アルゴリズム上ではサンプル生成がしやすい${\bf \Sigma}$の生成を行い，最終的な結果は${\bf \Psi}=h({\bf \Sigma})$のみについて考えることで，効率的に${\bf \Psi}$の事後推論を行うことができます．このように，パラメータの取り方を拡張することでMCMCの計算を簡単にする方法はパラメータ拡大として知られています．

アルゴリズム

${\bf \Sigma}$の事前分布として，逆ウィシャート分布$IW(\nu_0, {\bf S}_0^{-1})$を用います．これは密度関数が$|{\bf \Sigma}|^{-(\nu_0+p+1)/2} \exp\{-{\rm tr}({\bf S}_0^{-1}{\bf \Sigma}^{-1})/2 \}$に比例する分布です．

ギブスサンプリングの具体的な手順は以下の通りです．

${\bf Z}$が所与のもとで，${\bf \Sigma}$を$IW(\nu_0+n, ({\bf S}_0+{\bf Z}^\top {\bf Z})^{-1})$から生成し，${\bf \Psi}=h({\bf \Sigma})$とする．
${\bf \Sigma}$が所与のもとで，各$z_{ij}$を$(a_{ij}, b_{ij})$上で切断された正規分布

TN_{(a_{ij}, b_{ij})}\big({\bf \Sigma}_{j, -j} ({\bf \Sigma}_{-j,-j})^{-1}z_{i,-j}, \Sigma_{j,j}-{\bf \Sigma}_{j, -j} ({\bf \Sigma}_{-j,-j})^{-1}{\bf \Sigma}_{-j, j}\big),\\
a_{ij}={\rm max}(z_{kj}: y_{kj}<y_{ij}), \ \ \ b_{ij}={\rm min}(z_{kj}: y_{kj}>y_{ij})

から生成する．ただし，${\bf \Sigma}_{j,-j}$は${\bf \Sigma}$の$j$番目の行から$j$番目の列を取り除いたものを表す．

上のアルゴリズムでは，逆ウイシャート分布と一次元の切断正規分布からのサンプリングを繰り返すだけのため，既存の乱数生成パッケージなどを使って簡単に実行することができます．潜在変数$z_{ij}$に対する下限$a_{ij}$と上限$b_{ij}$は，それぞれ「$y_{ij}$より小さい値をもつデータに対応する潜在変数より$z_{ij}$は大きい」および「$y_{ij}$より大きい値をもつデータに対応する潜在変数より$z_{ij}$は小さい」という条件が反映されたものになっています．

ギブスサンプラーによって得られた${\bf \Psi}$の事後サンプルを用いることで，${\bf \Psi}$の事後推論を行うことができます．例えば，${\bf \Psi}$の逆行列${\bf \Psi}^{-1}$の事後平均や信用区間を見ることで，グラフィカルモデリングのように変数間の従属関係について推論を行うことが可能です．

おわりに

今回は従属構造をベイズ推定する方法について紹介しました．上記の説明では$z_i$の分布に多変量正規分布を用いているため，従属構造のモデルとして正規コピュラを推定する方法になっています．他のコピュラモデルを推定したい場合は，$z_i$に対する分布の仮定を適切に変える必要がありますが，MCMCのアルゴリズムは複雑になる可能性があります．

後半の記事では，上記のアルゴリズムをRで実装し，簡単な適用例を紹介したいと思います．

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．