不均一関数データに対する主成分分析と手書き文字データへの応用

不均一関数データに対する主成分分析と手書き文字データへの

応用

Principal Component Analysis for sparse functional data with

application to handwriting data

茅野光範

1∗

_堂園剛司

1

_小西貞則

2∗

Mitsunori Kayano

1

_{, Koji Dozono}

1

_{and Sadanori Konishi}

2

1

_{九州大学大学院数理学府}

1

_{Graduate school of Mathematics, Kyushu University}

2

_{九州大学大学院数理学研究院}

2

_{Faculty of Mathematics, Kyushu University}

Abstract: We consider principal component analysis for multi-dimensional sparse functional data. The mixed effect model and reduced rank model have been used for analyzing the sparse func-tional data. In this paper, we introduce a principal component method for the multi-dimensional sparse functional data based on the reduced rank model, and model selections will be performed by using Akaike information criterion (AIC) and Bayesian information criterion (BIC). Further-more, the use of the proposed method is illustrated through the analysis of human gait data and handwriting data.

1

はじめに

図 1 は，各年齢毎に観測された背骨のミネラル濃度 をプロットした図であり，図中の線分で結んだ点が各 個体の観測点，太線が全体の平均関数である．この図 のように観測・測定される多変量データは一般にアン バランスである．このようなデータは，スパースデー タ (sparse data) とも呼ばれ，スパースデータから各個 体に対する全体の曲線と主成分（重み）関数を推定す る方法が提案されている（Shi et al. (1996)，Rice and Wu (2001)，James et al. (2000) ）． 平滑化して関数主成分分析を行う通常の方法につい ては，いくつかの問題点がある．まず，小数の観測しか 持たない個体があるとき，各曲線に対する一意な表現を 作り出すことは出来ない．また，複数の特性に関する関 数データ（多次元の関数データ）を解析すると，主成分 の個数は次元数に比例して増加する（茅野他 (2006)）． 通常の方法の詳細や応用，拡張については Ramsay and Silverman (2005)を参照されたい． 曲線データの解析において，混合効果モデル（Mixed effect model, Laird and Ware (1982)）は広く用いられ ており，Shi et al. (1996) や Rice and Wu (2001) は

∗_{連絡先：九州大学大学院数理学府/数理学研究院}       〒 812-8581 福岡県福岡市東区箱崎 6-10-1        E-mail: [email protected] (茅野)        [email protected] (小西) 関数主成分問題を解くために混合効果モデルを用いる ことを提案した．混合効果モデルでは，各個体の関数 xi(t)を推定するために i 番目の個体だけでなく全ての 個体に対する観測値を用いることや，固定効果におけ る係数ベクトルとランダム効果ベクトルにおける分散 共分散行列の推定に最尤法を用いることの利点がある． しかし，このモデルでは，多くの局所的な最大値が存 在すること，データ点よりもパラメータ数が多くなる ために一意なパラメータ推定が出来ないことがある． 図 1: 不均一な成長曲線 (James et al. (2000))． 線分：各個体の観測点，太線：全体の平均関数．

人工知能学会研究会資料

SIG-DMSM-A603-05 (2/27)

80

120

160

−140 −80 −40

X Coordinate

Y Coordinate

0

2

4

6

0 50 150

Time (sec)

X Coordinate

0

2

4

6

−140 −80 −20

Time (sec)

Y Coordinate 図 2: 手書き文字データ．中図，右図のような XY 座標データを取り扱う． これに対して，James et al. (2000) は，特別な混 合効果モデルである縮小ランクモデル (reduced rank model)を提案し，これらの不適切問題を回避した．本 稿では，この縮小ランクモデルとその推定方法につい て述べ，多次元関数データへの拡張と情報量規準やベイ ズ型モデル評価基準によるモデル選択手法を提案する． 2節では縮小ランクモデルの定義と基底展開法につ いて述べ，このモデルが混合効果モデルの一部である ことを示し，多次元関数データへの拡張を試みる．EM アルゴリズムを用いた縮小ランクモデルのパラメータ 推定法を 3 節で述べ，4 節でモデル評価基準について 述べる．最後に，本稿で提案する手法を，特に図 2 の ような手書き文字データに適用する．

2

縮小ランクモデル

2.1

1

次元関数データに対するモデル

xi(t) (i = 1,· · · , N, t ∈ T ⊂ R) を時間 t における第 i曲線に対する値，µ(t) を全体の平均関数，ξj(t) (j = 1,· · · , k, t ∈ T ) を第 j 主成分（重み）関数とし，ξ(t) = (ξ1(t),· · · , ξk(t))′とする．このとき，各 xi(t)に対して 以下のモデルを仮定する． xi(t) = µ(t) + k ∑ j=1 ξj(t)αij+ εi(t) = µ(t) + ξ(t)′αi+ εi(t) . (1) ここで，ξj(t)は正規直交性の条件 ∫ T ξ(t)ξ(t)′dt = Ik， つまり， ⟨ξj, ξj′⟩ = ∫ T ξj(t)ξj′(t) dt = δjj′ (2) を満たすとする．ただし，Ikは k 次単位行列で，δjj′ = 1 (j = j), 0 (j ̸= j) である．ランダムベクトル αi = (αi1,· · · , αik)′ は i 番目の個体に対する主成分関数の 重みを与え，εi(t)はランダムな観測誤差である．また， αi, εi(t)は全て平均 0 で，αiは共通の共分散行列 D を 持ち，観測誤差は互いに独立で分散 σ2を持つとする． このモデルの適用において，データが限られた観測 時点でしか観測されていない場合には平均関数と主成 分関数に対して制約をおく必用がある．つまり，平均 関数と主成分関数に対して以下のように基底展開の仮 定をおく．ϕν(t) = (ϕν 1(t),· · · , ϕνM(t))′を時間 t におけ る基底関数のベクトルで，Θ と θµをそれぞれ基底展 開における係数で M × k 行列，M 次元ベクトルとす る．ここで，ν (> 0) は基底関数に含まれるハイパーパ ラメータである（2.3 節参照）．このとき， µ(t) = M ∑ m=1 θµmϕνm(t) = θ′µϕ ν_{(t) ,} ξj(t) = M ∑ m=1 θjmϕνm(t) = θ′jϕ ν (t) より，(1) 式は， xi(t) = ϕν(t)′θµ+ ϕν(t)′Θαi+ εi(t) , εi(t)∼ (0, σ2) , αi∼ (0, D) (3) となる．ここで，Θ = (θ1,· · · , θk)，ξ(t) = Θ′ϕν(t)で ある．また，(2) 式に示した主成分関数の正規直交性の 条件は， Θ′Θ = Ik , W := ∫ T ϕν(t)ϕν(t)′dt = IM (4) となる．上式のように基底関数の交差積行列（cross-product matrix）∫_T ϕν(t)ϕν(t)′dtを W で表す．

しかし，実際に観測されるのは離散データである．そ こで，各 i に対して，tijを観測値が得られている時点と し，xi = (xi(ti1),· · · , xi(tini))′, Φ ν i = (ϕ ν (ti1),· · · , ϕν(tini))′, εi= (εi(ti1),· · · , εi(tini))′とする．ここで， Φν i は i 番目の個体に対する基底行列であり，個体番号 iに依存していることに注意する．このようにして，縮 小ランクモデルは以下のように表される． xi= Φνiθµ+ ΦνiΘαi+ εi, Θ′Θ = Ik, W = ∫ T ϕν(t)ϕν(t)′dt = IM εi∼ (0, σ2I) , αi∼ (0, D) . (5) このモデルの適用は θµ, Θ, D, σ2の推定に帰着され， その推定方法は 3 節で述べる．また，実際には，基底 関数の個数 M と主成分の個数 k などを適切に選択す る必要がある（4 節参照）．

2.2

多次元関数データに対するモデル

p次元の関数データ{xi1(t),· · · , xip(t); i = 1,· · · , N, t ∈ T ⊂ R} の各 xil(t)に対して以下のモデルを仮定 する． xil(t) = µl(t) + k ∑ j=1 ξlj(t)αij+ εil(t) = µl(t) + ξl(t)′αi+ εil(t) . (6) ただし，主成分関数 (ξ1j(t),· · · , ξpj(t))は正規直交性 の条件∑pl=1 ∫ T ξl(t)ξl(t)′dt = Ik，つまり， p ∑ l=1 ∫ T ξlj(t)ξlj′(t)dt = δjj′ を満たすとする．また，ランダムな係数ベクトル αiは 変数の番号 l には依存しないものとする．(6) 式を l に ついてまとめると多次元関数データに対する縮小ラン クモデルは以下のようになる． xi(t) = µ(t) + Ξ(t)′αi+ εi(t) . (7) ここで，xi(t) = (xi1(t),· · · , xip(t))′, µ(t) = (µ1(t), · · · , µp(t))′, εi(t) = (εi1(t),· · · , εip(t))′, Ξ(t) = (ξ1(t), · · · , ξp(t))である．また，Ξ(t) は主成分関数の正規直 交性の条件から ∫ T Ξ(t)Ξ(t)′dt = p ∑ l=1 ∫ T ξl(t)ξl(t)′dt = Ik を満たす．従来の基底展開に基づく関数主成分分析で は主成分の個数は関数データの次元数 p と基底関数の 個数 M に依存するが，このモデルでは他の要因に依存 しないことに注意する． つぎに，第 l 変数に対応する平均関数と主成分関数 に対して基底展開 µl(t) = M ∑ m=1 θµlmϕνm(t) = θ′µlϕ ν (t) ξlj(t) = M ∑ m=1 θljmϕm(t) = θ′ljϕ ν (t) の仮定をおく．これを用いて (7) 式を変形する．まず， µ(t) = (µ1(t),· · · , µp(t))′, θµ= (θ′µ1,· · · , θ′µp)′, ξl(t) = (ξl1(t),· · · , ξlp(t))′, Θl = (θl1,· · · , θlk), Θ = (Θ′1, · · · , Θ′ p)′，Φν(t) = diag{ϕ ν (t)· · · ϕν(t)} とおくと， µ(t) = Φν(t)′θµ, ξl(t) = Θ′lϕ ν (t), Ξ(t) = Θ′Φν(t) となる．ただし，Φν_{(t)は pM× p ブロック対角行列で} ある．したがって，(7) 式は次のようになる． xi(t) = Φν(t)′θµ+ Φν(t)′Θαi+ εi(t) . 正規直交性の条件は， ∫ T Ξ(t)Ξ(t)′dt = Θ′ ∫ T Φν(t)Φν(t)′dt Θ = Ik , つまり，Θ′Θ = Ik， ∫ TΦν(t)Φν(t)′dt = IpM である． ここで，2 つ目の条件は，∫_Tϕν(t)ϕν(t)′dt = IMと同 等であることに注意する． また，各 i に対して tijを観測値が得られている時点と して，xi = (xi(ti1)′,· · · , xi(tini)′)′, Φ ν i = (Φν(ti1),· · · , Φν_(t ini))′, εi = (εi(ti1)′,· · · , εi(tini)′)′とすれば，多 次元関数データに対する縮小ランクモデルは以下のよ うに 1 次元と同様な形式で与えられ，その推定方法な どは 1 次元の場合と同様である． xi= Φνiθµ+ ΦνiΘαi+ εi Θ′Θ = Ik, ∫ T Φν(t)Φν(t)′dt = IpM εi∼ (0, σ2Ipni) , αi∼ (0, D) . (8)

2.3

正規直交な基底関数の構成

基底関数として，安道・井元・小西 (2001) によって 提案されたハイパーパラメータ ν (> 0) が付いたガウ ス型動径基底関数 ϕν,m∗(t) = exp { −(t− cm)2 2νs2 m } (m = 1,· · · , M) に基づいて正規直交な基底関数を構成する．ただし， cmは基底関数の位置，smは基底関数の広がりを表し，

ハイパーパラメータ ν は基底関数の位置をさらに調整 する． まず，上式で与えられるガウス型動径基底関数に対 して交差積行列を計算し，それをコレスキー分解した 結果得られた行列を Uν とする．つまり，ガウス型動 径基底関数のベクトル ϕν,∗(t) = (ϕν,₁∗(t),· · · , ϕν,_M∗(t))′ に対して， W := ∫ T ϕν,∗(t)ϕν,∗(t)′dt = U_ν′Uν となる行列 Uνを求める．ただし，ガウス型動径基底関 数の交差積行列 W は次式によって近似的に計算する． Wmn= ∫ T ϕνm(t)ϕνn(t) dt≈ ∫ ∞ −∞ ϕνm(t)ϕνn(t) dt = √ 2πνs2 ms2n √ s2 m+ s2n exp { − (cm− cn)2 2ν(s2 m+ s2n) } . ここで，Wmnは交差積行列 W の (m, n) 成分とする． このとき， ϕν(t) = (U_ν′)−1ϕν,∗(t) とすれば，この基底関数 ϕν(t)は正規直交の条件 (4 式) を満たす．このようにして，正規直交な基底関数 ϕν(t) を構成して用いる．行列 Uνの代わりに W−1/2を用い ても正規直交な基底関数が構成できるが，行列 Uνを用 いた正規直交化は，特にグラム・シュミットの正規直交 化に対応していることが知られている (堂園他 (2006))．

3

推定方法

(5)や (8) 式で与えられる縮小ランクモデルに含まれ る未知パラメータ θµ，Θ，σ2，D の推定方法を述べる． ここでは，1 次元縮小ランクモデルに対する最尤法を 述べる．多次元縮小ランクモデルに対する最尤法は 1 次元の場合と同様である． ランダムベクトル αiと誤差ベクトル εiはそれぞれ 正規分布に従っているとすると， xi ∼ N (Φνiθµ, σ2I + ΦνiΘDΘ′(Φνi)′) となる．したがって，{xi; i = 1,· · · , N} の同時分布か ら尤度が与えられるが，この尤度を θµ, Θ, σ2, Dにつ いて最大化することは，困難な最適化問題となる．そ こで，αi を欠測値と見なし EM アルゴリズムを用い て最尤推定量を求める．James et al.(2000) による EM アルゴリズムを用いた最尤推定の手順が http://www-rcf.usc.edu/ gareth/で公開されている．

4

モデル選択

観測データ{xi; i = 1,· · · , N} の同時分布は，基底 関数の個数 M と主成分の個数 k，基底関数に含まれる ハイパーパラメータ ν に依存する．そこで，モデル評 価基準として情報量規準 AIC とベイズ型モデル評価基 準 BIC を与え，それらに基づいて M と k，ν を選択 する．xiの密度を f (xi|ti; θµ, Θ, σ2, D)，P をモデル に含まれるパラメータ数としたとき，情報量規準 AIC とベイズ型モデル評価基準 BIC は，それぞれ次式で与 えられる． AIC(M, k, ν) =−2 N ∑ i=1 log f (xi|ti; ˆθµ, ˆΘ, ˆσ2, ˆD) + 2P , BIC(M, k, ν) =−2 N ∑ i=1 log f (xi|ti; ˆθµ, ˆΘ, ˆσ2, ˆD) + P log N . 様々な M と k に対して上式で与えられる AIC と BIC を計算し，最小値をとるときの M と k，ν を最適な値 として選択する．

5

適用例

5.1

歩行データ

人間の歩行時における腰角度や膝角度などを同時に 計測した離散データ (Olshen et al. (1989)) のうち，腰 角度データ (図 4) に対して 1 次元の縮小ランクモデル を適用する．このデータは個体毎に観測時点が共通で あることに注意する． 0.0 0.2 0.4 0.6 0.8 1.0 0 20 40 60

Proportion of a Gait Cycle

Hip angle (degree)

図 4: 腰角度データと推定曲線． 個体数 N = 39，観測数 n = ni= 20．

0.0

0.4

0.8

0.0

0.4

0.8

Time

X Coordinate

0.0

0.4

0.8

0.0

0.4

0.8

Time

Y Coordinate

図 3: 標準化後の手書き文字データ．個体数 N = 12，観測数 ni≃ 300∼600． 0.0 0.2 0.4 0.6 0.8 1.0 −2 −1 0 1 2

Proportion of a Gait Cycle

Hip Angle (degree)

1st PC curve 2nd PC curve PC Curves 図 5: 腰角度に対応する主成分関数． 実線が第 1，点線が第 2 主成分関数． 1次元の縮小ランクモデルを腰角度データに適用し， モデル選択を行った．基底関数の個数 M を M = 3, 4,· · · , 9，ガウス型動径基底関数のハイパーパラメータ ν を ν = 5, 10,· · · , 50，主成分の個数 k を k = 2, 3, · · · , 9 と したとき，AIC と BIC によってモデル選択を行った． AICを最小にしたのは M = 9, ν = 35, k = 7，BIC を最小にしたのは M = 7, ν = 20, k = 6 であった．そ こで，より簡単なモデルとなる M = 7, ν = 20, k = 6 を最適な値として選択した．このときの推定曲線と主 成分関数を図 4，5 に示す．各個体の腰角度の変動がう まく捉えられていることや，主成分関数がうまく推定 できていることがわかる．実際，ここで求めた主成分 関数は従来のように離散データを平滑化してから関数 主成分分析を行った結果と一致する．また，歩行デー タのような個体毎に観測時点が共通であるデータに対 しても縮小ランクモデルによって有益な情報を抽出で きると考えられる．

5.2

手書き文字データ

本稿で提案する多次元の縮小ランクモデルを図 2 に 示した手書き文字データに適用する．手書き文字デー タを関数データ解析の枠組みで取り扱った研究として は Ramsay (2000) などがある． データ取得のために，脇 (2005) を参考に手書き文字 データ取得のソフトウェアを作成した．その際，WA-COM社製ペンタブレット intuos 3 を用い，VisualBa-sic.NETによって作成した．図 2 では，ペンが空中に あるときの XY 座標を便宜上，共に 0 としたが，この 部分の XY 座標を観測されていないと見なし，2 次元 のスパース関数データに基づく主成分分析を適用する． 所得したデータは，九大数理統計グループの学生 12 人（男性・右利き・21∼26 歳）の「あ」のデータである． 個体数は N = 12，観測数は個体ごとに異なり ni≃ 300 ∼600 である．ここでは，X 座標と Y 座標と時間（お よそ 3∼6 秒）を全て [0, 1] として標準化を行った．図 3に標準化した後のデータを示す． つぎに，標準化を行なったデータに対して，まず，AIC と BIC の最小化によってモデル選択を行った．基底関 数の個数を M = 3, 4,· · · , 12，ハイパーパラメータを ν = 5, 10,· · · , 50，主成分の個数を k = 2, 3, · · · , 10 (≤ M )としたとき，M = 12，ν = 20，k = 10 のとき AICと BIC が共に最小となった．つぎに，M = 12， ν = 20，k = 10 として，手書き文字データの XY 座 標に対して 2 次元の縮小ランクモデルを適用した．こ こで，James らによる EM アルゴリズムを元にモデル の推定を行った．彼らの手順では，EM アルゴリズム のステップが終了してから行列 Θ の正規直交化を行っ

ているが，ここでは，ステップ毎に正規直交化をする ことによって，推定の安定，個体差を捉えられる，分 散の大きい順に主成分がとれるという利点があったの で，この方法を用いた．その結果，個体毎の推定曲線 は個体の変動を捉えていることが確認でき，推定され た平均関数も各個体の平均的な変動を捉えていた． いま，モデル評価基準 AIC と BIC によって最適な主 成分の個数は k = 10 個が選択された．この 10 個の主 成分の中での第 2 主成分までの累積寄与率を求めてみ ると，70.3 ％であった．そこで，第 1，第 2 主成分関数 を求め，それぞれの主成分に対する X 座標と Y 座標 の寄与率は，第 1 主成分に対して X の寄与が 34.7 ％， Y の寄与が 65.3 ％，第 2 主成分に対しては X の寄与 が 48.7 ％，Y の寄与が 51.3 ％となった．このことか ら，第 1 主成分に対しては Y 座標の寄与が大きく，第 2主成分に対しては両方の寄与を考える必要がある．こ れらを考慮して結果の解釈を行う．まず，第 1 主成分 に対しては Y 座標の特に時間 0.3∼0.7 の部分で 3ヶ所 大きな重みがかかっていた．また，第 2 主成分に対し ては，時間 0.2∼0.8，0.9 の間の所々に大きな重みがか かっていた．したがって，全体として，時間 0.2∼0.8， 0.9の部分に大きな重みがかかっている，つまり，個人 の癖が現れていると考えられる． ここでは，手書き文字データの解析結果に対して全 体的な解釈を行ったが，主成分関数によって重み付け られた所を部分的に解釈して，主成分の意味を考える ことも可能である．また，時間を [0, 1] に標準化せずに 時間方向（横方向）のシフト・レジストレーションを してから解析することや，手書き文字データの XY 座 標だけでなく筆圧などの情報も考慮した解析を行うこ とも考えられる．

6

おわりに

本稿では，縮小ランクモデルに基づいて，多次元の スパース関数データに対する関数主成分分析を提案し， 情報量規準 AIC やベイズ型モデル評価基準 BIC によっ てモデル選択を行った．さらに，提案手法を歩行データ と手書き文字データに適用し，分析を行った．歩行デー タへの適用ではスパースでない関数データに対しても 縮小ランクモデルは有効に機能し，手書き文字データ への適用では提案手法によって個人の癖を捉えること を試みた． 今後は，縮小ランクモデルに基づく正則化関数主成 分分析や正則化最尤法による推定方法を提案し，一般 化情報量規準 GIC や一般化ベイズ型モデル評価基準 GBICによるモデル選択を行いたい．

参考文献

[1] 安道知寛，井元清哉，小西貞則: 動径基底関数ネット ワークに基づく非線形回帰モデルとその推定．応用統計 学，Vol. 30, No. 1，pp. 19–35 (2001) [2] 堂園剛司，松井秀俊，川野秀一，小西貞則: 正規直交基 底に基づく関数データ化と自己組織マップによるクラス タリング．2006年度統計関連学会連合大会講演報告集， pp.272 (2006)

[3] James, G., Hastie, T., and Sugar, C.: Princi-pal component models for sparse functional data.

Biometrika，Vol. 87, pp.587–602 (2000)

[4] 茅野光範，小西貞則，平川英樹，久原哲: 正則化基底展

開法に基づく関数主成分分析とその応用．応用統計学，

Vol. 35, No. 1, pp. 1–16 (2006)

[5] Laird, N. and Ware, J.: Random-effects models for longitudinal data．Biometrics，Vol.38，pp.963– 74 (1982)

[6] Olshen, R. A., Biden, E. N., Wyatt, M. P. and Sutherland, D. H.: Gait analysis and the bootstrap．

Annals of Statistics，Vol. 17，pp. 1419–1440 (1989)．

[7] Shi, M., Weiss, R. and Taylor, J.: An analysis of paediatric cd4 counts for acquired immune defi-ciency syndrome using flexible random curves． Ap-plied Statistics，Vol. 45，pp. 151—64 (1996) [8] Ramsay, J. O.: Functional components of variation

in handwriting．Journal of American Statistical As-sociation，Vol. 95，pp. 9–15 (2000)

[9] Ramsay, J. O. and Silverman, B. W.: Functional

Data Analysis (Second Edition)．Springer (2005) [10] Rice, J. and Wu, C. . Nonparametric mixed effects

models for unequally sampled noisy curves． Biomet-rics，Vol. 57，pp. 253–259 (2001)

[11] 脇一真:自己組織化マップを用いたペン入力による個人

認証に関する研究．平成17年度大阪電気通信大学卒業