数理解析研究所講究録
第 1834 巻 2013 年 26-44

26

双対平坦空間の情報幾何を利用した統計的推定
東京大学・情報理工学系研究科

廣瀬善大

東京大学・情報理工学系研究科，理化学研究所・脳科学総合研究センター 駒木文保
Yoshihiro Hirose
Graduate School of Information Science and Technology, University of Tokyo
Fumiyasu Komaki

Graduate School of Information Science and Technology, University of Tokyo
Brain Science Institute, RIKEN
概要
文献 [5], [6], [7] において提案された統計的推定手法 bisector regression について解説する．これらの論文で考えられている手法は，Efron らによる least angle
regression の拡張を目指したものである．拡張にあたっては双対平坦空間の情報幾何を利用している．拡張により得られた bisector regression アルゴリズムは，least
angle regression の厳密な拡張ではないが，least angle regression 同様に角の二等分線を利用したパラメータ推定手法となっている．一般化線形回帰，ガウシアングラ
フィカルモデル，分割表の 3 つの問題を考え，特にガウシアングラフイカルモデルについてのアルゴリズムを紹介する．また，bisector regression の適用例も紹介する．

1 導入
本稿では，文献 [5], [6], [7] において提案された統計的推定手法 bisector regressionについて解説する．これらの論文で考えられている手法は，Efron らによる least angle

regression (LARS, [3]) の拡張を目指したものである．拡張にあたっては双対平坦空間の情報幾何 ([1], [2], [8]) を利用している．拡張により得ちれた bisector regression アルゴ
リズムは，LARS の厳密な拡張ではないが，LARS 同様に角の二等分線を利用したパラメータ推定手法となっている．

LARS は線形回帰問題におけるパラメータ推定 (回帰係数の推定) のために提案されたアルゴリズムであり，同時にモデル推定も行う．アルゴリズムの出力としてパラメータの推定値とモデルの系列が得られる．説明変数間の相関を内積 (角度) とみなし，推定量の構成において反応変数との相関が大きな (角度の小さな) 説明変数を順次追加することで推定値の系列を構成する．LARS アルゴリズムはユークリツド幾何により記述され，幾何学的に単純であり統計学的にも自然に解釈される．また，LARS は正則化の代表的な手法である lasso ([13], [14]) とも関連がある．正則化とは，最尤推定量の過学習を避けるために，尤度に正則化項 (罰則化項) を加えたものを最大化する手法である ([4]). LARS 自27
体は正則化ではないが，LARS アルゴリズムを修正することで lasso にょる推定値を出力できることが知られている．

双対平坦空間の情報幾何はユークリッド幾何を一般化したものである．確率分布族のなす空間を考え，計量や接続といった構造を入れる．統計学でょく現れる Fisher 情報を計量とし，-接続と呼ばれる接続の族を考える． -接続には互いに双対な構造をもつ組がある．本稿では，その中でも特に，1-接続と $-1$ -接続を利用する．1-接続と $-1$ -接続はそれ$\alpha$ぞれ$e$$\alpha$-接続， -接続とも呼ばれる．リーマン接続に対応するのは -接続であり，自己双
$m$

$0$

対な接続である．接続が与えられると測地線や平坦さを扱うことができ，-接続と -接続のそれぞれに関する測地線等があらわれる．また，-接続に関して平坦な空間は -接続の意味でも平坦であることが知られており，その逆も成り立つ．このことから互いに双対な接続に関して平坦な空間は双対平坦空間と呼ばれる．与えられた接続に関して平坦な空間には便利な座標系 (アファイン座標系) が存在し，双対平坦空間には -接続と -接
続というふたつの接続に関するアファイン座標系がそれぞれ存在する．指数型分布族は1-接続に関して双対平坦空間になることが知られている．自然パラメータが アファイン座標系になっており，期待値パラメータが m- アファイン座標系になっている．情報幾
何は，統計学においては高次の漸近論などで有用であり，さらに制御論や最適化などの分
$m$

$e$

$m$

$e$

$m$

$e$

$\pm$

$e-$

野でも有用である．
文献 [5], [6], [7] において bisector regression による推定が扱われた問題は 3 つある．

一般化線形回帰，ガウシアングラフィカルモデル，分割表である．一般化線形回帰は線形回帰を含むより広い問題のクラスであり，LARS の直接的な拡張を目指した問題である．
3 つの問題の共通点は，いずれも指数型分布族の自然パラメータの推定問題であることである．本稿では特にガウシアングラフィカルモデルについての bisector regression を説明する．

本稿の構成は以下の通りである．第 2 節において bisector regression を適用する 3 つの問題を紹介する．それぞれの問題で推定するパラメータを定義し，後の節のための設定を確認する．第 3 節では bisector regression アルゴリズムの解説に必要な情報幾何の説明を行う．指数型分布族は双対平坦空間をなし，bisector regression アルゴリズムはその双対平坦空間内の入れ子構造を利用して推定を行う．アルゴリズムの記述とその解説は第4 節において行う．ガウシアングラフィカルモデルをとりあげ，まずアルゴリズムの幾何学的な説明を行い，その後 bisector regression アルゴリズムを紹介する．第 5 節で本稿のまとめをする．

28

2 問題
提案手法 bisector regression が適用される問題について説明する．文献 [5] にある通り，3 つの問題を考える．扱う問題は，一般化線形回帰，ガウシアングラフイカルモデル，

分割表である．いずれの場合も，bisector regression はパラメータ推定とモデル推定を同時に行う．

2.1 一般化線形回帰
一般化線形回帰におけるパラメータ推定は，一般化線形モデル (generalized linear
model, GLM) のパラメータ推定であり ([10]), 対応する指数型分布族の中から分布をひとつに決定する問題である．
のとする．計画行列 $X$ は行列である．ただし，$\{y_{a}, x^{a}=(x_{1}^{a}, x_{2}^{a}, \ldots, x_{d}^{a})\}_{a=1,2,\ldots,n}$が観測されているも$X=(x_{i}^{a})_{1}\leq a\leq n,$ $1\leq i\leq d=(x_{1}, x_{2}, \ldots, x_{d})$$x_{i}=(x_{i}^{1}, x_{i}^{2}, \ldots, x_{i}^{n})^{T}(i=1,2, \ldots, d)$で与えられ，$n\cross d$

とおいた．1 を$n$

個の 1 を$n\cross(d+1)$ 行列を $\tilde{X}=(1|X)$ と定もつベクトル，すなわち $1=(1,1, \ldots, 1)^{T}$ とし，義する．
指数型分布族は$f(y| \xi)=\exp(\sum_{a=1}^{n}y_{a}\xi^{a}+\sum_{b=1}^{r}u_{b}(y)\xi^{b+n}-\psi(\xi))$と表される．ただし，$\xi’:=(\xi^{1}, \xi^{2}, \ldots, \xi^{n})^{T},$$(\xi^{1}, \xi^{2}, \ldots, \xi^{n+r})^{T}=(\xi^{\prime T}, \xi^{\prime/T})^{T},$$\xi":=(\xi^{n+1}, \xi^{n+2}, \ldots, \xi^{r+n})^{T},$$\xi:=$$u(y);=$
とおくことにする．
に関する十分等計量であり， は正規化定数$y=(y_{1}, y_{2}, \ldots, y_{n})^{T}$$\{y_{1}, \ldots , y_{n}, u_{1}(y), u_{2}(y), \ldots , u_{r}(y)\}$はに対応する凸関数である．パラメータ

$\xi$

$\xi$

$\psi(\cdot)$

は自然パラメータと呼ばれ，$\psi(\cdot)$

は

$\xi$

のポテン

シャル関数と呼ばれる．一般化線形回帰問題においては，さらに$\xi’=\tilde{X}\theta’, \xi"=\theta", \theta\in\Theta$

として$\theta’$を推定する．ただし，$\theta’:=(\theta^{0}, \theta^{1}, \ldots, \theta^{d})^{T},$$\theta:=(\theta^{0}, \theta^{1}, \ldots, \theta^{d+r})^{T}=(\theta^{\prime T}, \theta^{J/T})^{T}$

る．十分等計量 $u(y)$ の期待値$\mu_{b+n}=E[u_{b}(y)](b=1, \ldots, r)$$\theta":=(\theta^{d+1}, \theta^{d+2}, \ldots, \theta^{r+d})^{T},$とおく．

$\Theta\subset R^{d+r}$$\mu=(\mu_{1}, \ldots, \mu_{n+r})^{T}$

はパラメータ空間であは $\mu_{a}=E[y_{a}](a=1, \ldots, n)_{\backslash },$で与えられ，期待値パラメータと呼ばれる．

$r>0$ の場合，推定したいのは$\theta’$

であり$\theta"$は局外パラメータである．パラメータのうち bisector regression アルゴリズムにより推定されるのは$\theta’$だけである．局外パラメー

29
タの扱いについては，第 4 節においてガウシアングラフィカルモデルを例にとって説明する．一般化線形回帰においても同様な手続きが必要になる．

一般化線形回帰の例としてロジスティック回帰を考える．考える指数型分布族は$f(y| \xi)=\prod_{a=1}^{n}\frac{\exp(y_{a}\xi^{a})}{1+\exp\xi^{a}}=\exp(\sum_{a=1}^{n}y_{a}\xi^{a}-\psi(\xi))$$y\in\{0,1\}^{n},$ $r=0$ でである．ただし，は自然パラメータであり，ある． に関するポテンシャル関数は $\psi(\xi)=\sum_{a=1}^{n}\log(1+\exp\xi^{a})$ である．一般化線形$\xi=(\xi^{1}, \xi^{2}, \ldots, \xi^{n})^{T}$$\xi$

回帰問題として $\xi=X\theta,$ $\Theta=R^{d},$$\theta=(\theta^{0}, \theta^{1}, \ldots, \theta^{d})^{T}$

は $\mu_{a}=E[y_{a}]=\exp\xi^{a}/(1+\exp\xi^{a})(a=1,2, \ldots, n)$を考える．期待値パラメータで与えられる．パラメータ$\theta$

$\mu$

を推

定することが問題である．

2.2

ガウシアングラフィカルモデル

ガウシアングラフィカルモデル (Gaussian graphical model, GGM) を考える．問題

は独立グラフと呼ばれるグラフにおける辺選択である．独立グラフとは，二変数の間の条件付き独立性に対応して辺が与えられたグラフである．多変量正規分布の精度行列の非対角成分の推定により独立グラフにおける辺の判定を行う．本稿では無向グラフのみを考える．グラフィカルモデルについては，文献 [9] などを参照のこと．
$X$ は多変量正規分布にしたがうものと変量確率変数 $X=(X_{1}, \ldots, X_{p})^{T}$ を考え，$\Sigma^{-1}=(\sigma^{ab})$

$p$

が他の変
する．二変数の条件付き独立性は次のように定義される．すなわち， と$\sigma^{ab}=0$数が与えられた下で条件付き独立であるとは，が成り立つことである．ただし，$X_{a}$$\Sigma=(\sigma_{ab})$は共分散行列であり，$\Sigma^{-1}=(\sigma^{ab})$は共分散行列$\Sigma$$X_{b}$の逆行列で精度行列と呼$X$ の各成分に対してグラフの頂点をひとつずつ対応させたもばれる．独立グラフとは，ので，二頂点$X_{a}$と$X_{b}$

り決定される (図 1).

$p$の間の辺の有無は条件付き独立の場合には辺を与えないことによ
変量正規分布にしたがう確率変数 $X$ に対応する独立グラフは，$P$個の頂点をもち最大 $p(p-1)/2$ 本の辺をもつグラフである．確率分布に多変量正規分布(ガウス分布) を仮定しているので，この独立グラフを利用して得られるモデルはガウシ

アングラフィカルモデルと呼ばれる．つまり，GGM とは，独立グラフにおいて辺のない
二変数が条件付き独立であるような多変量正規分布の集合である．グラフにおいて辺がある二変数については，対応する GGM に含まれる分布に条件付き独立性に関する制約は課されない．
問題は GGM における辺選択である．条件付き独立性の定義から，精度行列

対角成分 $\sigma^{ab}(a<b)$ の推定を行うことになる．観測として

$p$

変量変数が

$n$

$\Sigma^{-1}$の非サンプル分与30$\Sigma^{-1}=$

図1$(00**$$****$$0***$$0***)$$\Leftrightarrow$

確率変数が 4 変量の場合 $(p=4)$ における精度行列と独立グラフの例． は任意
$*$

において
の値を表す．精度行列と独立グラフにおいては，2 つの辺となっている．

$\sigma^{13}=\sigma^{14}=0$$\Sigma^{-1}$$e_{13}$$e_{14}$$\Sigma^{-1}$に対応するが除去されている．

えられているものとする．あるいは観測の標本平均と標本共分散行列が与えられているものとする．平均と$\Sigma^{-1}$の対角成分は局外パラメータである．興味のあるパラメータと局外パラメータの扱いの違いについては，GGM を例にとり第 4 節で説明する．

2.3 分割表
分割表の例を表 1 に示す．ここでは二元分割表を例にして説明する．因子を 2 つ以上もつ場合についても同様に扱うことができる．
2 つの因子$X_{1}$と$X_{2}$からなる分割表を考える．

$X_{1}$$X_{1}=0,1,$個の水準をもつものとする．すなわち，

は $m+1$ 個の水準を，$\ldots,$$X_{2}$$m,$$X_{2}=0,1,$ $\ldots,$ $n$

は $n+1$. 考えるのは多項分布であり，その確率分布は$f(y|p)= \frac{N!}{y_{00}!y_{01}!\ldots y_{mn}!}p_{00}^{y00}p_{01^{1}}^{y0}\ldots p_{mn}^{y_{mn}}.$ただし， は総観測数，陶はセル $(i,j)$ の観測数でありを満たす．
また， はセル $(i,j)$ の確率で$N$$\sum_{i=0}^{m}\sum_{j=0}^{n}y_{ij}=N$を満たす．

$\sum_{i=0}^{m}\sum_{j=0Pij}^{n}=1$$Pij$

確率分布の対数をとると，$\log f(y|p)=\sum_{i=0}^{m}\sum_{j=0}^{n}y_{ij}\log p_{ij}+\log\frac{N!}{y_{00}!y_{01}!\ldots y_{mn}!}$$= \sum_{i=1}^{m}\sum_{j=1}^{n}y_{ij}\log\frac{p_{ij}p_{00}}{p_{i0}p_{0j}}+\sum_{i=1}^{m}(\sum_{j=0}^{n}y_{ij})\log\frac{p_{i0}}{p_{00}}$$+ \sum_{j=1}^{n}(\sum_{i=0}^{m}y_{ij})\log\frac{p_{0j}}{p_{00}}+N\log p_{00}+\log\frac{N!}{y_{00}!y_{01}!\ldots y_{mn}!}$31
表1$y_{ij}$$X_{2}:n+1$ 個の水準をもつ因子，$X_{1}:m+1$ 個の水準をもつ因子，
分割表の例．
$N$: セル $(i,j)$ の観測数， : 総観測数．

と変形できる。 $i=1,2,$$\ldots,$$m,$$j=1,2,$$\ldots,$

$n$

に対して

$\theta_{X_{1}}^{i}=\log\frac{p_{i0}}{p_{00}}, \theta_{X_{2}}^{j}=\log\frac{p_{0j}}{p_{00}}, \theta_{X_{1}X_{2}}^{ij}=\log\frac{p_{ij}p_{00}}{p_{i0}p_{0j}}$とおき，自然パラメータ$\theta=(\theta_{X_{1}};\theta_{X_{2}};\theta_{X_{1}X_{2}})$$=(\theta_{X_{1}}^{1}, \theta_{X_{1}}^{2}, \ldots, \theta_{X_{1}}^{m};\theta_{X_{2}}^{1}, \theta_{X_{2}}^{2}, \ldots, \theta_{X_{1}}^{n};\theta_{X_{1}X_{2}}^{11}, \theta_{X_{1}X_{2}}^{12}, \ldots, \theta_{X_{1}X_{2}}^{mn})$を推定する．ここで， は 次元ベクトル， は 次元ベクトル，
ベクトルである．bisector regression により推定したいパラメータは$\theta_{X_{1}}$$m$$\theta_{X_{2}}$$n$$\theta_{X_{1}X_{2}}$$\theta_{X_{1}X_{2}}$と$\theta_{X_{2}}$は $mn$ 次元

であり，$\theta_{X_{1}}$は局外パラメータである．

3 情報幾何
本節では bisector regression アルゴリズムに必要な情報幾何の説明を行う．双対平坦空間の情報幾何はユークリッド幾何の一般化のひとつである．双対平坦空間においては，
ユークリッド空間と同様に直線や距離に対応する量を扱うことができる．そのためユーク

リッド幾何により記述されていた LARS アルゴリズムの拡張を目指すことができる．本稿では指数型分布族のみを扱う．情報幾何の詳細については，文献 [1], [2], [8] などを参照のこと．

指数型分布族のなす空間 $S=\{f(x|\theta)|\theta\in\Theta\}$ を考える．指数型分布族は$f(x| \theta)=\exp(\sum_{a=1}^{n}u_{a}(x)\theta^{a}-\psi(\theta))$32と表される． は自然パラメータ， は$\theta$$\psi$$\theta$$u(x)=(u_{a}(x))$ はに関するポテンシャル関数，

十分統計量と呼ばれる．指数型分布族のパラメータとして，十分統計量の期待値をとった
期待値パラメータ $\eta=E[u(x)]$ も有用である．ここで，$\partial_{i}=\partial/\partial\theta^{i},$$\partial^{i}=\partial/\partial\eta_{i}$と略記することにする．計量として Fisher 情報$g_{ij}=E[(\partial_{i}\log f(X|\theta))(\partial_{j}\log f(X|\theta))]$を用いる．-接続と呼ばれる接続の族を考え，特に $\alpha=\pm 1$ に注目する．1-接続と $-1$ -接$\alpha$$m$ -接続とも呼ばれ，互いに双対な構造をもつ．自然パラメータ
続はそれぞれ -接続，
$e$

に関する

$\alpha$$\theta$-接続の係数螺は$\Gamma_{ijk}^{(\alpha)}=E[(\partial_{i}\partial_{j}\log f(X|\theta)+\frac{1-\alpha}{2}\partial_{i}\log f(X|\theta)\partial_{j}\log f(X|\theta))\partial_{k}\log f(X|\theta)]$

で与えられる．

一般に，あるパラメータに関して接続の係数がすべてあると言われ，そのパラメータは$-\alpha$$\alpha$$0$であるとき空間 $S$ は$\alpha$-平坦で-アファイン座標系と呼ばれる． が -平坦であれば$S$$\alpha$-平坦であり，その逆も成り立つ．このことから， -平坦な空間を双対平坦空間と呼ぶ．いま$\pm\alpha$$S$は指数型分布族であり， -接続に関して双対平坦空間になることが知られて$\pm 1$いる．自然パラメータが $e-$ アファイン座標系であり，期待値パラメータが m- アファイン
座標系である．

期待値パラメータ$\eta$

に関するポテンシャル関数$\phi$

を $\phi(\eta)=E[\log f(x|\eta)]$ で定義する．

このとき，関係式$\partial_{i}\psi=\eta_{i}, \partial^{i}\phi=\theta^{i}, \partial_{i}\partial_{j}\psi=g_{ij}, \partial^{i}\partial^{j}\psi=g^{ij},$$\phi(\eta)+\psi(\theta)-\eta\cdot\theta=0$が成り立つ．ただし， は Fisher 情報行列の逆行列，-接続に関する測地線，-測地線を定義する． の曲線 :$(g^{ij})$$\alpha$$(g_{ij})$$S$$\alpha$るとは，任意の $t\in[0,1],$$k$$\eta\cdot\theta=\sum\eta_{i}\theta^{i}$$t\mapsto\gamma(t)$$\gamma$が$\alpha$である．

-測地線であに対して$\ddot{\gamma}^{k}(t)+\sum_{i,j}\dot{\gamma}^{i}(t)\dot{\gamma}^{j}(t)(\Gamma_{ij}^{(\alpha)k})_{\gamma(t)}=0$が成り立つことである．ただし，は$\Gamma_{ij}^{(\alpha)k}$$(\gamma^{1}(t), \gamma^{2}(t), \ldots, \gamma^{d}(t))$

である．

$\alpha=1$

は曲線

とすると，

$\theta_{(1)}$と$\gamma(t)$$\theta_{(2)}$$\Gamma_{ij}^{(\alpha)k}$$=$$\sum_{\iota}\Gamma_{ij\downarrow}^{(\alpha)}g^{\iota k}$

で与えられ，の成分である．測地線は直線に対応するものを結ぶ -測地線は，e- アファイン座標系である自然$e$

パラメータが$\theta_{\gamma}(t)=t\theta_{(1)}+(1-t)\theta_{(2)}, t\in[0,1]$33と表される曲線である．同様に，$\eta_{(1)}$と$\eta_{(2)}$を結ぶ $m$ -測地線は，期待値パラメータが$\eta_{\gamma}(t)=t\eta_{(1)}+(1-t)\eta_{(2)}, t\in[0,1]$と表される曲線である．
点$\eta_{(1)}$から点$\theta_{(2)}$への Kullback-Leibler ダイバージェンス ( $KL$ ダイバージェンス) は$D( \eta_{(1)}|\theta_{(2)})=E_{\eta_{(1)}}[\frac{\log f(X|\eta_{(1)})}{\log f(X|\theta_{(2)})}]$$KL$ ダイバージェンスはユークリッド距離の二乗に対応する量である．
で定義される．
$KL$

ダイバージェンスはポテンシャル関数$\psi,$$\phi$

を用いて$D(\eta_{(1)}|\theta_{(2)})=\phi(\eta_{(1)})+\psi(\theta_{(2)})-\eta_{(1)}\cdot\theta_{(2)}$と表すことができる．

多変量正規分布$(\mu_{1}, \Sigma_{1})$と$(\mu_{2}, \Sigma_{2})$との $KL$ ダイバージェンス $D(\mu_{1}, \Sigma_{1}|\mu_{2}, \Sigma_{2})$ は$D( \mu_{1}, \Sigma_{1}|\mu_{2}, \Sigma_{2})=-\frac{1}{2}\log|\Sigma_{1}|+\frac{1}{2}\log|\Sigma_{2}|+\frac{1}{2}tr(\Sigma_{1}\Sigma_{2}^{-1})-\frac{p}{2}$$+ \frac{1}{2}(\mu_{2}-\mu_{1})^{T}\Sigma_{2}^{-1}(\mu_{2}-\mu_{1})$で与えられる．第 4 節で考える bisector regression では期待値が同じ多変量正規分布については，$\Sigma_{1}$から$\Sigma_{2}$$\mu$は固定される．期待値

への $KL$ ダイバージェンスは$D( \Sigma_{1}|\Sigma_{2})=-\frac{1}{2}\log|\Sigma_{1}|+\frac{1}{2}\log|\Sigma_{2}|+\frac{1}{2}tr(\Sigma_{1}\Sigma_{2}^{-1})-\frac{p}{2}$で与えられる．
$S’$ を $S$
最後に $m$ -射影を定義する．
$p\in S$ とする．点の部分空間とし，$p$の $S’$ への m-射影とは， からのダイバージェンスが最小であるような $S’$ の点である．
$P$

4 推定手法
第 2 節において bisector regression を適用する問題として 3 つの問題を紹介した．ア

ルゴリズムの基本的なアイデアはどの問題においても同じであるが，問題ごとに異なる点がある．本節で bisector regression のアルゴリズムを解説するに当たっては，3 つの問
題のうち特にガウシアングラフィカルモデル (GGM) に注目することにする．線形回帰，

一般化線形回帰では確率変数の期待値を扱っていたが，第 2.2 節で確認した通り GGM の辺選択においては確率変数の共分散行列を扱う．したがって，線形回帰と同じく正規分布

34を扱うが，考える問題の構造は異なっている．多変量正規分布は指数型分布族であり双対平坦空間になっていた．GGM の辺選択問題は，情報幾何を用いることで，一般化線形回帰とは異なった形で自然に扱うことのできる問題である．

推定の手続きとしては，まず局外パラメータを固定し，残されたパラメータを bisector

regression アルゴリズムで推定するという手続きをとる．局外パラメータの固定は，線形回帰における原点の設定に対応するものであり，計画行列の正規化作業を一般化したものと考えてよい．ただし，bisector regression においてはパラメータのスケーリングについては不変性があり，bisector regression を適用するモデルを設定するということが局外パ
ラメータの固定の幾何学的意味である．

$u=\Sigma^{-1}\mu$と$\Sigma^{-1}$と共分散行列をパラメータとしてもつ．あるいはをパラメータとすることもできる．前者が期待値パラメータであり，多変量正規分布は平均ベクトル$\mu$$\Sigma$

後者が自然パラメータである．2 種類のパラメータの双対性を利用して局外パラメータを固定する．つまり，期待値パラメータである

$\mu$

と

$\Sigma$の対角成分を最尤推定値に固定し，の非対角成分を推定する．パラメータの数は $d=p(p-1)/2$
である．これにより局外パラメータと推定するパラメータが直交し，bisector regression
自然パラメータである$\Sigma^{-1}$による推定を自然に行うことができるようになる．

以下では，まず局外パラメータの固定について説明を行い，その後 bisector regression

アルゴリズムを紹介する．GGM に特化した bisector regression アルゴリズムなので，便宜上 BRGGM (Bisector Regression for Gaussian Graphical Models) と呼ぶことにする．

4.1 局外パラメータの固定
BRGGM アルゴリズムは，多変量正規分布モデルのうち，期待値と分散が固定されたサブモデルにおいて推定値を更新するアルゴリズムである．期待値ベクトル と共分散
$\mu$

行列

$\Sigma$の対角成分 (分散$\sigma_{aa}$) をそれぞれ $\mu=\hat{\mu}MLE$ と $\sigma_{aa}=(\hat{\Sigma}_{MLE})_{aa}(a=1, \ldots,p)$として固定する．ただし，$\hat{\mu}_{MLE}$と$(\hat{\Sigma}_{MLE})_{aa}$はそれぞれ最尤推定量であり，実際には標本平均と標本分散である．標本共分散行列を$M=\{\Sigma|\sigma_{aa}=(\hat{\Sigma}_{MLE})_{aa}(a=1, .

\ell.

,p)\}$

とおく

とおく．サブモデル $M$ を$\hat{\Sigma}_{MLE}$

なお，平均は固定されている(図 2).

ので表記を省略する．サブモデル $M$ は双対平坦空間であり，BRGGM アルゴリズムが動く空間である．表記上， と$\Sigma$$\Sigma^{-1}$はともに共分散行列$\Sigma$の多変量正規分布を表すものとする．BRGGM アルゴリズムは $M$ における推定量の更新 (移動) として記述される．

最初の推定量$\Sigma_{(0)}^{-1}^{\wedge}$はフルモデルの最尤推定量$\Sigma_{MLE}^{-1}^{\wedge}$

であり，独立グラフ

完全グラフ GMLE に対応する．最後に出力される推定量$\Sigma^{-1}^{\wedge}$$(d)$$\hat{G}_{(0)}$としてはは各変数がそれぞれ独立で

35

図2

$M$
BRGGM アルゴリズムが適用されるサブモデル．

$(\hat{\Sigma}_{MLE})_{aa}(a=1, \ldots,p)\}.\hat{\Sigma}_{MLE}$ ;ある． ;$\hat{\Sigma}_{0}$$=$$\{\Sigma|\sigma_{aa}$標本共分散，BRGGM の最初の推定量$\hat{\Sigma}_{(0)}$

$=$

でも

各変数が独立であるというモデル $\{\Sigma|\sigma_{ab}=0(a\neq b)\}$ での最尤推定量であり，BRGGM が最後に出力する推定量$\hat{\Sigma}_{(d)}$はこれと一致する．BRGGM アルゴ

リズムが適用されるサブモデル $M$ と独立モデル $\{\Sigma|\sigma_{ab}=0(a\neq b)\}$ は直交する．
$\hat{\Sigma}_{MLE}$と$\hat{\Sigma}_{0}$はともに $M$ に含まれる．

あるというモデル $\{\Sigma|\sigma_{ab}=0(a\neq b)\}$ での最尤推定量
は辺をひとつももたないグラフ$\hat{G}_{0}$$\Sigma_{0}^{-1}^{\wedge}$

であり，独立グラフとなる．なお， の対角成分と$\hat{\Sigma}_{0}$$\hat{\Sigma}_{MLE}$$\hat{G}_{(d)}$の対角成分の

値は一致する．また，BRGGM アルゴリズムが適用されるサブモデル $M$ と独立モデル$\{\Sigma|\sigma_{ab}=0(a\neq b)\}$4.2が直交することが知られている．

Bisector regression アルゴリズム

GGM の辺選択に対して bisector regression アルゴリズムを説明する．本節の説明はサブモデル $M$ におけるものであり，平均と分散は固定されているので省略する．精度行
列$\Sigma^{-1}$の非対角成分 $\sigma^{ab}(a\neq b)$ を推定し，対応する独立グラフ $G$ を推定することが目$(a, b)\in I\subseteq\{(a’, b’)|1\leq a’<b’\leq p\},$ $s\in R$ に対して，ふたつのサブモデル
的である．

$M(I)$ と $M((a, b), s, I)$ を$M(I)=\{\Sigma|\sigma^{a’b’}=0((a’, b’)\not\in I)\},$$M((a, b), s, I)=\{\Sigma|\sigma^{ab}=s, \sigma^{a’b’}=0((a’, b’)\not\in I)\}$との間に辺をもたなにより定義する．サブモデル $M(I)$ は $(a’, b’)\not\in I$ に対して
い独立グラフに対応するモデルである．サブモデル $M((a, b), 0, I)$ は $(a’, b’)\not\in I$ に対し$X_{a’}$て $X_{a’}$ と$X_{b’}$

が辺をもたず，さらに$X_{a}$と$X_{b}$$X_{b’}$

も辺をもたないような独立グラフに対応す36
るモデルである．
まず BRGGM アルゴリズムで何をしようとしているのかについて幾何学的な説明をし，その後に BRGGM アルゴリズムを紹介する．BRGGM アルゴリズム，あるいはガ
ウシアングラフィカルモデル以外の問題における bisector regression アルゴリズムについても，線形回帰問題におけるパラメータ推定モデル推定の手法である least angle
regression (LARS) の拡張を目指したものである．アルゴリズムの反復ごとにパラメータの推定値を出力し，推定値の入るモデルが入れ子状に順に小さくなっていくのが特徴である．これによりアルゴリズムの反復はパラメータの次元と同じ回数だけであることが分かり，考えるモデルの候補数が大幅に削減される．

LARS では，説明変数と反応変数の相関を内積 (角度) として扱い，反応変数となす相関の大きい (角度の小さい) 説明変数を推定量の構成に順次取り込んでいた．すでに推定量を構成している説明変数同士は反応変数となす相関 (角度) が互いに一致しており，幾何学的には推定量が角の二等分線上を動いているものとみなせる．双対平坦空間の情報幾何を利用して LARS の拡張を行うにあたって，角の二等分線を距離を用いた形で一般化

し，bisector regression はその曲線上を推定量が動くようなアルゴリズムになっている．
ただし，bisector regression は LARS の厳密な拡張ではない．LARS では原点を出発点にしてアルゴリズムが始まり，フルモデルの最尤推定量を推定量の最終的な到達点にしていた．しかし，bisector regression では，推定量はフルモデルの最尤推定量を出発し，最終的に原点にたどりつくようになっている．これはユークリッド空間を双対平坦空間に一般化する際に生じる困難を避けるためであった．推定量の移動 (更新) によりパラメータの推定値とモデルの組の系列を出力する点は同じである．

4.1 節において局外パラメータの固定について説明した．BRGGM アルゴリズムはサ$M$ 内を動く推定量の点を考え，ブモデル $M$ に対して働くアルゴリズムである．
$M(I)$の形で表される $M$ のサブモデルに推定量が順次動いていく手続きを記述する．図 3

は推定量の点が角の二等分線に対応する曲線に沿って移動する様子を表したものであり，bisector regression アルゴリズムではこの手続きを反復することになる．推定量の点が$k$

回目の反復で得られた推定値$\hat{\Sigma}_{(0)}=\hat{\Sigma}_{MLE}$である．精度行列に対して，サブモデル$M(I\backslash (a, b))$にいるものとする．アルゴリズム開始時には個の非零成分をもち，$|I|=d-k$ を満たに対してとなっている．各 $(a, b)\in I$$\Sigma_{(k)}^{-1}^{\wedge}$

すある $I\subseteq\{(a, b)|1\leq a<b\leq p\}$$\hat{\Sigma}_{(k)}$は $d-k$$\Sigma_{(k)}^{-1}^{\wedge}\in M(I)$を考える．サブモデル$M(I\backslash (a, b))$

$M(I)$ よりもは，一次元だけ次元の小さなサブモデルで $M(I)$ に含まれている．現在の推定値
サブモデル $M(I\backslash (a, b))$ への $KL$ ダイバージェンスを測り，その最小値をにする．

$\hat{\Sigma}_{(k)}$からの $KL$ ダイバージェンスが$t^{*}$$t^{*}$$\hat{\Sigma}_{(k)}$

から各

とおくこと

になるように各サブモデル $M(I\backslash (a, b))$

37
$M(I)$ における各を平行移動したものは，それぞれ $M((a, b), s_{ab}^{*}, I)$ の形で表される．
$M((a, b), s_{ab}^{*}, I)$ の交点をに対と決める．このとき， を達成した
してであるから，推定値を含むサブモデルの次元がひとつ小さ$\hat{\Sigma}_{(k+1)}$$(a^{*}, b^{*})\in I$$t^{*}$$\hat{\Sigma}_{(k+1)}\in M(I\backslash (a^{*}, b^{*}))$

くなったことが分かる．なお，上述の説明では角の二等分線に対応する曲線を陽に用いていないが，曲線の両端がそれぞれ推定値$\hat{\Sigma}_{(k)}$と$\hat{\Sigma}_{(k+1)}\iota_{\llcorner}^{-}$

対応している．

角の二等分線に対応する曲線について図 3 にしたがって簡単に説明する．現在の推定値から各サブモデル $M(I\backslash (a, b))$ への $m$-射影を$\hat{\Sigma}_{(k)}$への $KL$$\overline{\Sigma_{(k)}}^{ab}$$(a^{*}, b^{*})\in I$

ダイバージェンスのうち最小値がである．

$\hat{\Sigma}_{(k)}$から$\tilde{\Sigma}_{(k)}^{-a’b’}$$t^{*}$

とおく．

$\overline{\Sigma_{(k)}}^{ab}$$\hat{\Sigma}_{(k)}$から各であり，最小値を達成するのはへの $KL$ ダイバージェンスが$\hat{\Sigma}_{(k)}$から$\overline{\Sigma_{(k)}}^{a^{*}b^{*}}$

へを平行移動しダイバージェンスと一致するように，サブモデルて $M((a’, b’), s_{ab’}^{*}, I)$
はの $M((a’, b’), s_{ab’}^{*}, I)$ へのとする．ただし，の $KL$$M(I\backslash \{(a’, b’)\})$$\tilde{\Sigma}_{(k)}^{-a’b’}$

影であり，$s_{ab’}^{*}$$\hat{\Sigma}_{(k)}$を満たす値として決まる．

は $D$ $(\Sigma_{(k)}^{\wedge}|\tilde{\Sigma})=t^{*}$$(a^{*}, b^{*})$

同様のことを考えると，$s_{a^{*}b^{*}}^{*}=0,\tilde{\Sigma}_{(k)}^{-a^{*}b^{*}}=\overline{\Sigma_{(k)}}^{a^{*}b^{*}}$

定値$\hat{\Sigma}_{(k+1)}$$m-fl1\backslash$に対してであることが分かる．新しい推は $M((a’, b’), s_{a’ b’}^{*}, I)$ と $M((a^{*}, b^{*}), 0, I)=M(I\backslash \{(a^{*}, b^{*})\})$ との交点で$\ovalbox{\tt\small REJECT}$あった．ここで，$0\leq t\leq t^{*}$

なる

$t$

に対して， の代わりに を使って
$t^{*}$$t$$\tilde{\Sigma}_{(k)}^{-a’b’}$と$\tilde{\Sigma}_{(k)}^{-a^{*}b^{*}}$に対応する点を作ることができる．そして，それらをそれぞれ含むようなサブモデル$M((a’, b’), s_{a’b’}(t), I)$と $M((a^{*}, b^{*}), s_{a^{*}b}*(t), I)$ の交点

3 における点線上の点になる).
$\hat{\Sigma}(t)=\hat{\Sigma}_{(k+1)}$$t=0$ のとき
特に，$\hat{\Sigma}(t)$を考えることができる (図$\hat{\Sigma}(t)=\hat{\Sigma}_{(k)}$であり，$t=t^{*}$のときである．拡張ピタゴラスの定理 ([1], [2]) より，$D(\Sigma_{(k)}^{\wedge}|\tilde{\Sigma}^{-a’b’}(t))=D(\Sigma_{(k)}^{\wedge}|\tilde{\Sigma}^{-a^{*}b^{*}}(t))$,$D(\Sigma^{-a’b’}(t)\sim|\hat{\Sigma}(t))=D(\Sigma^{-a^{*}b^{*}}(t)\sim|\hat{\Sigma}(t))$

が成り立つ．このことから$\{\hat{\Sigma}(t)|0\leq t\leq t^{*}\}$

が角の二等分線に対応する曲線になっていることが分かる．

BRGGM アルゴリズムは以下の通りである．ステップ 2 から 6 が反復され，反復ごとにパラメータの推定値とモデルが出力される．独立グラフを $G=(V, E)$ とおき， は頂$V$

点集合を， は辺集合を表すものとする． は固定されており，辺集合
$E$

$V$

される．ステップ 4 の$l_{(k)}^{-ab}$は，$\hat{\Sigma}_{(k)}$と$\overline{\Sigma_{(k)}}^{ab}$$E$の要素が推定を結ぶ $m$ -測地線である．

BRGGM
入力 :

標本平均 $\mu_{sample}=\hat{\mu}_{MLE}$ , 標本共分散

出力 :推定値の系列1.

$\hat{\Sigma}_{(0)}$$\Sigma_{(1)}^{-1}^{\wedge},$$\Sigma_{(0)}^{-1},$$\Sigma_{(d)}^{-1}^{\wedge}\wedge$$\ldots,$$:=\hat{\Sigma}_{MLE},\hat{G}_{(0)}$$:=\hat{G}_{MLE},$$I$,$\Sigma_{sample}=\hat{\Sigma}_{MLE}$

独立グラフの系列$\hat{G}_{(0)},\hat{G}_{(1)},$$\ldots,\hat{G}_{(d)}$$:=\{(a, b)|1\leq a<b\leq p\},$ $k$ $:=0$ とおく．

38
$M(I\backslash$$\tilde{\Sigma}_{(k)}^{-a^{*}b^{*}}$$\hat{\Sigma}_{0}$$M(I\backslash \{(a^{*}, b^{*})\})=M((a^{*}, b^{*}), 0, I)$

図3

推定量の更新．
$I$$\{(a, b) 1 \leq a < b \leq p\}$ :$\subseteq$$\Sigma^{-1}$の非対角$M((a, b), s, I)$
成分のうち推定されるべきものを指定するインデツクス．
$\{\Sigma|\sigma^{ab}=s,$$\sigma^{a"b"}=0((a", b")\not\in I)\}.\hat{\Sigma}_{0}$ ;番目の推定量．

$\hat{\Sigma}_{(k+1)}:k+1$$M((a, b), 0, I)$

番目の推定量．

影でも同じ．

独立モデルの最尤推定量．

$\overline{\Sigma_{(k)}}^{ab}$における最尤推定量，あるいは;$\hat{\Sigma}_{(k)}$サブモデル $M(I\backslash \{(a, b)\})=$の $M(I\backslash \{(a, b)\})$ への $m$ -射$D(\Sigma_{(k)}^{\wedge}|\overline{\Sigma_{(k)}}^{a’b’})>D(\Sigma_{(k)}^{\wedge}|\overline{\Sigma_{(k)}}^{a^{*}b^{*}})$$M((a’, b’), s_{a’b’}^{*}, I)$ への $m$より決まる値．

$\hat{\Sigma}_{(k)}$$\tilde{\Sigma}_{(k)}^{-a^{*}b^{*}}=\overline{\Sigma}_{(k)}^{-a^{*}b^{*}}$

は-射影． : 条件$s_{ab}^{*}$$\hat{\Sigma}_{(k+1)}$である．

$\tilde{\Sigma}_{(k)}^{-a’b’}$:$\hat{\Sigma}_{(k)}$の$D(\Sigma_{(k)}^{\wedge}|\tilde{\Sigma}_{(k)}^{-a’b’})=D(\Sigma_{(k)}^{\wedge}|\overline{\Sigma_{(k)}}^{a^{*}b^{*}})$$M(I\backslash \{(a’, b’)\})$

である．

$=$

$\hat{\Sigma}_{(k)};k$

よりも$M(I\backslash \{(a^{*}, b^{*})\})$に近いのでは $M((a^{*}, b^{*}), 0, I)$ と $M((a’, b’), s_{a’b’}^{*}, I)$ の交点として定義される．点線は角の二等分線に対応する曲線である．BRGGM アルゴリ
ズムでは，推定量がこの曲線に沿って動いていると見ることができる．

2. 各 $(a, b)\in I$ に対してモデル $M(I\backslash \{(a, b)\})$ の最尤推定量$(a^{*}, b^{*}):=$ arg min3.

$\overline{\Sigma_{(k)}}^{ab}$を計算する．

$(a,b)\in ID(\Sigma_{(k)}^{\wedge}|\overline{\Sigma_{(k)}}^{ab})$$t^{*}:= \min_{(a,b)\in I}D(\Sigma_{(k)}^{\wedge}|\overline{\Sigma_{(k)}}^{ab}),$

とおく．

4.

各 $(a, b)\in I$

に対して，二条件を満たすような

5.

$(\mathfrak{a}, b)\in I$$s_{ab}^{*}$

に対して

と

$\tilde{\Sigma}_{(k)}^{-ab}\in l_{(k)}^{-ab}$$\hat{\sigma}_{(k+1)}^{ab}:=s_{ab}^{*},$$\hat{E}_{(k+1)}:=\hat{E}_{(k)}\backslash \{e_{a^{*}b^{*}}\}$

6.

7.

$D(\hat{\Sigma}_{(k)}|\tilde{\Sigma}_{(k)}^{-ab})=t^{*},\tilde{\Sigma}_{(k)}^{-ab}\in M((a, b), s_{ab}^{*}, I)$を求める．
$(a, b)\not\in I$

に対して$I:=I\backslash \{(a^{*}, b^{*})\}$$:=\emptyset$

とする．

とする．

$k:=k+1,$$k+1<d-1$ の場合，
る．
$k+1=d-1$ の場合，ステツプ 7 へ進む．
$\hat{\Sigma}_{(d)}:=\hat{\Sigma}_{0},\hat{E}_{(d)}$$\hat{\sigma}_{(k+1)}^{ab}:=0$としてステップ 2 へ戻とし，アルゴリズム終了．

$(a^{*}, b^{*})\in I$ であるがステップ 5 では，のことにより，1 回の反復において$\Sigma^{-1}$$\hat{\sigma}_{(k+1)}^{ab^{*}}=0$であることに注意が必要である．この非対角成分がひとつ$0$になり，モデルが入れ

39
子状に小さくなっていくことになる．

4.3 例
一般化線形回帰とガウシアングラフィカルモデルについての例を紹介する．計算には統計解析のフリーソフトウェア $R$ を利用した ([12]).

4.3.1 一般化線形回帰の例
一般化線形回帰のひとつとしてロジスティック回帰を考える．ロジスティック回帰については 2.1 節で簡単に説明した．一般化線形回帰におけるリンク関数は $g(\mu)=$
で与えられ，考える確率分布は二項分布のサブモデルである．ロジス
ティック回帰は，線形回帰問題に対する正則化の代表的手法である lasso を形式的に適$\log\{\mu/(1-\mu)\}$

用することのできる問題である ([11]). しかし，lasso の適用は最適化問題として定式化することはできるものの，ロジスティック回帰問題の構造を取り入れたものではないことに注意が必要である．一方，bisector regression は確率分布のなす双対平坦空間の情
報幾何に基づいており，問題の構造を利用して推定を行っている．一般化線形回帰に対する bisector regression アルゴリズムを BRGLM (Bisector Regression for Generalized
Linear Models) と呼ぶことにする．
データとして South African Heart Disease (SAHD) データを使用した ([4]). 正則化

手法 lasso をロジステイック回帰に適用した文献 [11] においても SAHD データを解析している([11], Fig. 2). 使用したデータはる．SAHD データは 9 個の説明変数$x_{1},$ $x_{2},$分のデータが集められている．反応変数$x_{2}$: tobacco,

alcohol,$x_{9}$: ldl,: adiposity,

$R$

$y$

の ElemStatLearn パッケ$\ldots$, xg と反応変数$y$$-\grave{\backslash }\grave{\grave{}}$に含まれてい$n=462$ 人からなり，は chd である．説明変数はそれぞれ: family history,$x_{1}$: sbp,: obesity,:: age である．BRGLM により出力される推定値は 10 個であるが，説明変数$x_{3}$$x_{4}$$x_{5}$$x_{6}$: type $A,$$x_{7}$$x_{8}$の組合せは全部で $2^{9}=512$ 通りある．BRGLM により考えるモデル数が劇的に減少していることが分かる．

SAHD データに対する BRGLM の結果を図 4 に示した．横軸は BRGLM が出力した各推定値から原点までの $KL$ ダイバージェンスの平方根を示している．原点については 4.1

節において説明した．縦軸は説明変数 $x_{i}(i=1,2, \ldots, 9)$ の回帰係数$\theta^{i}(i=1,2, \ldots, 9)$の値を表している． GLM の出発点であるフルモデルの最尤推定値が図 4 の右端に対$B\grave{R}$応する．BRGLM アルゴリズムは図の右端から左端に向かって進む．図に入っている縦線はそれぞれ BRGLM アルゴリズムが出力した推定値に対応している．左に進むにしたがって非零成分の個数がひとつずつ減っており，縦線の上の数字は非零成分の個数を示し

-

40

2

$0$

6

$\sqrt{div}4$

SAHD データのロジステイック回帰に対する BRGLM の結果．横軸は各推定値

図4から原点までの $KL$ ダイバージエンスの平方根を示しており，ユークリツド距離に対応するものである．縦軸は説明変数 $x_{i}(i=1,2, \ldots, 9)$ の回帰係数 $\theta^{i}(i=1,2, \ldots, 9)$ の値を表す．図の右端はフルモデルの最尤推定値であり，BRGLM アルゴリズムの出発点を表している．左端は BRGLM アルゴリズムの出力した最後の推定値であり，推定
したいパラメータがすべて になっている．図に入っている縦線はそれぞれ BRGLM$O$

アルゴリズムが出力した推定値に対応している．左に進むにしたがって非零成分の個数がひとつずつ減っていく．縦線の上の数字は非零成分の個数を示している．パラメー
タが$0$になる順番は$\theta^{8},$$\theta^{4},$$\theta^{1},$$\theta^{7},$$\theta^{3},$$\theta^{2},$$\theta^{6},$$\theta^{5},$$\theta^{9}$

であった．

である．
ている．図 4 より， の成分が になった順番はPark and Hastie の方法 ([11]) を SAHD データに適用した結果を図 5 に示す．Park$\hat{\theta}$$0$$\theta^{8},$$\theta^{4},$$\theta^{1},$$\theta^{7},$$\theta^{3},$$\theta^{2},$$\theta^{6},$$\theta^{5},$$\theta^{9}$and Hastie の方法は lasso を一般化線形回帰に拡張した正則化の方法である． の$R$

glmpath パッケージを利用して計算を行った．図 5 の横軸は推定値の$l_{1^{-\nearrow}}$

ルムを，縦軸はパラメータ $\theta^{i}(i=1,2, \ldots, 9)$ の値を示している．縦線の上の数字は非零成分の個数を示している．図 5 の右端から左端に見ていくと，パラメータが$\theta^{8},$$\theta^{4},$$\theta^{7},$$\theta^{1},$$\theta^{6},$$\theta^{3},$$\theta^{2},$$\theta^{5},$$\theta^{9}$$O$になる順番はであったことが分かる．

BRGLM と Park and Hastie による出力は似ているが，以下のような違いがある．
BRGLM により出力されるパラメータの推定値は単調性をもつことが分かる．Park andの推定値の変化から明らかである．まHastie の方法は単調性をもたない．このことは$\theta^{9}$

た，BRGLM による推定値は符号が変化しない．Park and Hastie の方法においては推定

41

00

05

1.0

1.5

2.

$0$

$|\theta|$

SAHD データに対する Park and Hastie の方法の結果．横軸は推定値のノルムを，縦軸はパラメータの値を示している．縦線の上の数字は非零成分の個数を示している．図の右端から左端に見ていくと，パラメータが 0 になる順番は
図5$l_{1^{-}}$$\theta^{i}$$\theta^{8},$$\theta^{4},$$\theta^{7},$$\theta^{1},$$\theta^{6},$$\theta^{3},$$\theta^{2},$$\theta^{5},$$\theta^{9}$であった．

値の符号は変化する．

上述の点はアルゴリズムの定義から分かる点であるが，結果には次の違いもある．ふたつの手法で出力されるモデル (説明変数の組合せ) は異なっている．このことは，例えば，6 個の説明変数をもつモデルを比較すると分かる．BRGLM にょるモデルには$x_{7}$が含まれているが，Park and Hastie の方法によるモデルには が含まれていない．実際，図 4と図 5 において 6 と書かれた縦線に対応する推定値を確認すると，BRGLM においては$x_{7}$$\hat{\theta}_{(3)}^{7}\neq 0$であるが Park and Hastie の方法では$\hat{\theta}_{PH}^{7}=0$となっていることが分かる．

4.3.2 ガウシアングラフィカルモデルの例
GGM の辺選択の例を紹介する．精度行列の非対角成分を推定し，対応する独立グラ

フを確認する．アルゴリズムの出力として，独立グラフの系列を紹介する．比較のため
graphical lasso ([15]) と呼ばれる手法の出力も載せている．

データは数学成績データを利用した． の SIN パッケー に含まれる mat 化 arks デー
タである．このデータは 88 人の学生の 5 科目の試験成績を集めたものである．5 科目と$R$$\grave{}\sqrt{}\grave{}\grave{}\grave{}$は，1: 力学，2: 線形代数，3: 代数，4: 解析，5: 統計である．辺の総数は $d=5\cdot 4/2=10$

42

数学成績データに対する BRGGM の結果．5 つの科目は，1: 力学，2: 線形代
数，3: 代数，4: 解析，5: 統計．BRGGM により出力される系列では，左上から右下ま
図6で辺がひとつずつ少なくなっていく．

である．

数学成績データに対する BRGGM アルゴリズムの出力を図 6 に示す．出力としては 11個のグラフが得られるが，最後のグラフは辺をもたないグラフであり図 6 では省略した．
BRGGM アルゴリズムは左上から右下に進む．同じデータに対する graphical LASSOの出力 (の一部) を図 7 に示す．ふたつの手法による出力は，4 本の辺をもつグラフが異なっておりまったく同じではないが，似た系列になっている．線形回帰の LARS は lassoと関連があることは第 1 節で述べた．この結果は，LARS の拡張を目指した bisector

regression の lasso との関連を示唆するものであると考えられる．

5

まとめ

本稿では，least angle regression (LARS) を双対平坦空間の情報幾何にもとついて拡張した結果 ([5], [6], [7]) についての解説を行った．扱った問題は，一般化線形回帰，ガウ
シアングラフィカルモデル，分割表の 3 つであった．推定手法である bisector regression

アルゴリズムの説明は，特にガウシアングラフイカルモデルの場合について行った．推定
手続きは，局外パラメータの固定と bisector regression アノ 1/ ゴリズムによる推定の二段階であった．LARS では角の二等分線が重要な働きをしたが，bisector regression においても双対平坦空間において対応する曲線を考え，その曲線を利用して推定を行った．アルゴ
リズムの紹介とその幾何学的説明を行った．また，bisector regression アルゴリズムの適

43

34
図7

数学成績データに対する graphical lasso の結果．1: 力学，2: 線形代数，3: 代

数，4: 解析，5: 統計．Lagrange 乗数が大きくなるにしたがって，graphical LASSO
による出力は左上から右下に進む．

用例を紹介した．

参考文献[1] S. Amari (1985). Differential-Geometrical Methods in Statistics. Springer Lecture

Notes in Statistics, vol. 28, Springer.
[2] S. Amari and H. Nagaoka (2000). Methods of Information Geometry. Translationsof Mathematical Monographs, Vol. 191, Oxford University Press.
[3] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani (2004). Least Angle Regression (with discussion), The Annals of Statistics, vol. 32, 407-499.
[4] T. Hastie, R. Tibshirani, and J. Friedman (2009). The Elements of Statistical
Learning; Data Mining, Inference, and Prediction, 2nd Edition. Springer, NewYork.
[5] Y. Hirose (2012). Dimension Reduction Based on the Geometry of Dually Flat
Spaces. Doctoral thesis, University of Tokyo.
[6] Y. Hirose and F. Komaki (2010). An Extension of Least Angle Regression Basedon the Information Geometry of Dually Flat Spaces, Journal of Computational andGraphical Statistics, vol. 19, 1007-1023.
[7] Y. Hirose and F. Komaki (2012). Edge Selection Based on the Geometry of Dually44
Flat Spaces for Gaussian Graphical Models, Statistics and Computing,[8] R. Kass and P. Vos (1997). Geometrical Foundations of Asymptotic Inference.
John Wiley, New York.
[9] S. L. Lauritzen (1996). Graphical Models. Clarendon Press, Oxford.
[10] P. McCullagh and J. A. Nelder (1989). Genemlized Linear Models, 2nd Edition.
Chapman and Hall/CRC, Boca Raton.
-Regularization Path Algorithm for Gen[11] M. Y. Park and T. Hastie (2007).
eralized Linear Models, Journal of the Royal Statistical Society, Series , vol. 69,$L_{1}$$B$

659-677.
[12] R Development Core Team (2012).

$R$: A Language and Environment for Statis-

tical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN
3-900051-07-0, http: $//www.R$-project.org/.
[13] R. Tibshirani (1996). Regression Shrinkage and Selection via the Lasso, Journalof the Royal Statistical Society, Series$B$, vol. 58, 267-288.

[14] R. Tibshirani (2011). Regression Shrinkage and Selection via the Lasso: a Retrospective, Journal of the Royal Statistical Society, Series , vol. 73, 273-282.
$B$[15] M. Yuan and Y. Lin (2007). Model Selection and Estimation in the GaussianGraphical Model. Biometrika, vol. 94, 19-35.