はじめに
東京大学/株式会社Nospare リサーチャーの栗栖です．
この記事から数回に分けて，多様体上に分布するデータの統計分析に用いられる手法について紹介します．今回の記事では特にBhattacharya and Patrangenaru (2005)を参考に，多様体データの「平均」あるいは「期待値」に対応するフレシェ平均に関して知られている結果について紹介します．
多様体上に分布するデータのフレシェ平均
$(M,d)$を距離関数$d$をもつ($m$次元)リーマン多様体とし，$X_1,\dots, X_n$を$M$上に値をとり，共通の分布$Q$に従う独立なデータとします．このとき，フレシェ関数$\tilde{F}(p) = \mathbb{E}[d^2(p,X_1)]$を最小化するような$M$の要素の集合
\tilde{E}= \left\{\tilde{p} \in M: \tilde{F}(\tilde{p}) = \inf_{p \in M}\tilde{F}(p)\right\}
を考えましょう．この集合はフレシェ平均集合とも呼ばれます．$\tilde{E}$の要素が一点のみとなるとき($\tilde{E} = \{\mu\}$)，$\mu$を$Q$のフレシェ平均と呼びます．$M$がユークリッド空間の場合($M = \mathbb{R}^p$)，距離関数をユークリッド距離とするとフレシェ平均は確率ベクトルの期待値$\mu = \mathbb{E}[X_1]$に一致することから，フレシェ平均は通常の意味での確率ベクトルの平均値(期待値)の概念の拡張となっています．
$M$上のデータ$X_1,\dots, X_n$からフレシェ平均を推定する場合は$\tilde{E}$の標本対応を考えます．
\tilde{E}_n = \left\{\tilde{p} \in M: \tilde{F}_n(\tilde{p}) = \inf_{p \in M}\tilde{F}_n(p)\right\},\ \tilde{F}_n(p) = {1 \over n}\sum_{i=1}^n d^2(p,X_i)\ (\text{標本フレシェ関数}).
以下では$\tilde{E}_n$要素が一点のみとなる場合($\tilde{E}_n = \{\mu_n\}$)を考え，$\mu_n$の性質についてみていきます．
標本フレシェ平均の中心極限定理
$\mu$を$M$上の分布$Q$のフレシェ平均とします．また$\mu$と$0 \in \mathbb{R}^m$の適当な近傍$\tilde{U}$，$U$に対して，$\mu$における$M$の接空間から$M$への指数写像$\rm{Exp}_\mu: U \to \tilde{U}$を考えます．このとき，指数写像は$\rm{Exp}_\mu(0) = \mu$を満たし，さらに適当な回数連続微分可能な微分同相写像とします．
指数写像を用いて(標本)フレシェ関数を改めて定義しておきます．
F(x) = \mathbb{E}[d^2({\rm Exp}_\mu(x),X_1)],\ F_n(x) = {1 \over n}\sum_{i=1}^n d^2({\rm Exp}_\mu(x),X_i)
指数写像を用いたフレシェ関数を考えることで，$M$上でのフレシェ関数$\tilde{F}$の最小化問題($M$の幾何学的な構造に依存して複雑になる)がユークリッド空間上での最小化問題に変わり，計算が容易になることに加え，理論的にも扱いが簡単になります．
$F(x)$に対して以下を仮定します：
- $F(x)$は$U$上で２回連続微分可能．
- $F$のヘッセ行列$H=\mathbb{E}[H(0,X_1)]$逆行列をもつ．ここで，
$H(x,y) = ({\partial^2 \over \partial x_j \partial x_k}d^2({\rm Exp}_\mu(x),y))_{1\leq j,k \leq m}$です．
以上の設定の下で以下の結果が成り立つことが知られています(Bhattacharya and Patrangenaru (2005), Corollary 2.1)
\sqrt{n}x_n \stackrel{d}{\to} N_m(0, H^{-1}GH^{-1}).
ここで，$x_n = {\rm Exp}_\mu^{-1}(\mu_n)\in \mathbb{R}^m$, $G = \mathbb{E}\left[G(0, X_1)G(0,X_1)'\right]$，$G(x,y) = ({\partial \over \partial x_1}d^2({\rm Exp}_\mu(x),y),\dots,{\partial \over \partial x_m}d^2({\rm Exp}_\mu(x),y))'$です.
多様体の例
多様体の代表例として(m次元)単位球面が挙げられます．$m=1$の場合は円周上のデータとなり，風向のデータや犯罪が発生した時刻などデータが分析の対象となります．$m=2$の場合は球面上のデータとなり，地球の磁気の方向データや宇宙線が地球に飛来する方向などデータが分析の対象となります．
M = \mathbb{S}^m = \{p \in \mathbb{R}^{m+1}: \|p\| = 1\}.
ここで$\|\cdot\|$はユークリッド距離です．単位球面上の距離としては以下の測地線距離(geodesic distance)を用います.
d_g(p,q) = \text{arccos}(\langle p,q \rangle),\ p,q \in \mathbb{S}^m.
ここで$\langle \cdot, \cdot \rangle $はユークリッド空間における通常の内積です．また$p \in \mathbb{S}^m$における指数写像は以下の形で与えられることが知られています．
{\rm Exp}_p(x) = \cos(\|x\|)p + \sin(\|x\|){x \over \|x\|},\ x \neq 0.
その他の多様体データの例については Bhattacharya and Patrangenaru (2014)で紹介されているので興味にある方は参照してください．
統計分析へ応用
上記の結果は，多様体上のデータであってもフレシェ関数のヘッセ行列が正則であれば，確率ベクトルに対する通常の中心極限定理と同じ結果が成り立つことを意味しています．この結果により多様体上のデータに対してもフレシェ平均に関する信頼領域(confidence region)の構成や2標本検定を行うことができます．
信頼領域
標本フレシェ平均の中心極限定理が成り立つ場合に漸近分散に出てくる行列$G$, $H$はその標本対応で推定することができます．
\hat{G} = {1 \over n}\sum_{i=1}^nG(\mu_n,X_i)G(\mu_n,X_i)',\ \hat{H} = {1 \over n}\sum_{i=1}^n H(\mu_n, X_i)
適当な条件の下で$\hat{G}$, $\hat{H}$がそれぞれ$G$, $H$の一致推定量であることがわかります．このとき，$\mu_n$の中心極限定理と合わせて以下が成り立ちます．
nx_n'\hat{\Gamma}^{-1}x_n \stackrel{d}{\to} \chi_m^2.
ここで$\hat{\Gamma}^{-1} = \hat{H}\hat{G}^{-1}\hat{H}$, $\chi_m^2$は自由度$m$のカイ二乗分布です．この結果を用いてフレシェ平均$\mu$の$100(1-\alpha)$%信頼領域$U_{n,1-\alpha}$が以下のように構成できます．
\begin{align*}
U_{n,1-\alpha}&={\rm Exp}_\mu(C_{n,1-\alpha}),\\
C_{n,1-\alpha} &= \{z \in {\rm Exp}_{\mu_n}^{-1}(\tilde{U}): n(x_n - z)'\hat{\Gamma}^{-1}(x_n - z) \leq \chi_{m,1-\alpha}^2\}.
\end{align*}
ここで$\chi_{m,1-\alpha}^2$は自由度$m$のカイ二乗分布の$1-\alpha$分位点です．また$C_{1-\alpha}$を用いることでフレシェ平均に関する帰無仮説$H_0: \mu = \nu$の検定を有意水準$1-\alpha$で実行することができます．実際，$\nu \notin U_{n,1-\alpha}$であれば帰無仮説$H_0$を棄却すればよいことになります．Bhattacharya and Patrangenaru (2014)では古地磁気学への応用として，地球上における仮想地磁気極(virtual geomagnetic pole)の平均方向の推定と信頼領域の計算結果が紹介されています．
2標本検定
$\{X_i\}_{i=1}^n$を$M$上の分布$Q_1$に従う独立同分布のデータ，$\{Y_i\}_{i=1}^n$を$M$上の分布$Q_2$に従う独立同分布のデータとし，$\{X_i\}_{i=1}^n$と$\{Y_i\}_{i=1}^n$は独立であるとします．
さらに$\mu_n, \nu_n$をそれぞれ$Q_1$, $Q_2$の標本フレシェ平均とし，それぞれに対して中心極限定理が成り立つとします．このとき，$\Gamma_1, \Gamma_2$をそれぞれ$\mu_n, \nu_n$の漸近分散として以下が成り立ちます．
\sqrt{n}\{(x_n-y_n) -(x-y)\} \stackrel{d}{\to} N(0, \Gamma_1 + \Gamma_2).
ここで$x = {\rm Exp}_{\mu_0}^{-1}(\mu), y = {\rm Exp}_{\mu_0}^{-1}(\nu), x_n = {\rm Exp}_{\mu_0}^{-1}(\mu_n), y_n = {\rm Exp}_{\mu_0}^{-1}(\nu_n) \in \mathbb{R}^m$，$\mu_0 \in \tilde{U}$です．この結果を用いて$Q_1, Q_2$のフレシェ平均$\mu,\nu$に関する2標本検定を実行することができます．実際，帰無仮説を$H_0: \mu = \nu$の検定を有意水準$1-\alpha$で実行するには$\gamma_n=x_n - y_n$として集合
\begin{align*}
U_{n,1-\alpha}&={\rm Exp}_{\mu_0}(C_{n,1-\alpha}),\\
C_{1-\alpha} &= \{p \in \mathbb{R}^m: n(\gamma_n - p)'(\hat{\Gamma}_1 + \hat{\Gamma}_2)^{-1}(\gamma_n - p) \leq \chi_{m,1-\alpha}^2\}.
\end{align*}
について$0 \notin U_{n,1-\alpha}$であれば帰無仮説$H_0$を棄却すればよいことになります．Bhattacharya and Patrangenaru (2014)では生物学への応用として，オスとメスのゴリラで頭蓋骨の形状が異なるかどうかの検定結果が紹介されています．
まとめ
この記事ではBhattacharya and Patrangenaru (2005)をもとに多様体上に分布するデータの統計分析について紹介しました．株式会社Nospareには多様体データ分析に限らず，統計学の様々な分野を専門とする研究者が所属しています．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．
参考文献
[1] Bhattacharya, R. and Patrangenaru, V. (2005) Large sample theory of intrinsic and extrinsic sample means on manifolds. Annals of Statistics 33, 1225-1259.
[2]Bhattacharya, R. and Patrangenaru, V. (2014) Statistics on manifolds and landmarks based image analysis: A nonparametric theory with applications. Journal of Statistics Planning and Inference 145, 1-22.