関数データに対する主成分分析 ～方法論の紹介とRによる実装～ 滋賀大学データサイエンス学部 松井 秀俊

2 • 関数データ解析手法の１つである 関数主成分分析について Rによる実装法と方法論を紹介 • 主成分の重みが関数として与えられ 各主成分の特徴を表す （その解釈は分析者が行う） • 主成分得点を求めることで 曲線の特徴を低次元のベクトルに 変換できる • 得られた主成分から，少ない情報で 元の関数データを復元できる 本資料の内容 ――：第１関数主成分 - - -：第２関数主成分 子どもの身長のデータと 関数主成分分析による重み関数 主成分得点plot （横：第１，縦：第２） 全体の中で全期間通じて 身長が高い人ほど 横軸の値が大きく 成長期が遅い人ほど 縦軸の値が大きい

3 • 関数データ解析の概要 • データの関数化 • 多変量データに対する主成分分析 • 関数主成分分析 • Rによる関数主成分分析の実装 目次

4 • 右図のように，観測個体（各子供）それぞれが 時間等の経過に伴い繰り返し計測された 形式のデータを考える • このようなデータに対して各個体を時間の関数 として表し，観測データ（ベクトルデータ）の 代わりに関数をデータとして扱う方法は 関数データ解析とよばれている • 関数データ解析を適用することで，観測時点が 不均一であったり，観測時点・時点数が 個体ごとに異なっていたとしても 比較的容易に分析ができる 関数データ解析 (Functional Data Analysis; FDA) 子供の身長の推移データ Ramsay & Silverman (2005)

5 • 関数データ解析では，古典的な統計手法を関数データの枠組みへ拡張した ものが多く研究されている • 関数の特性を利用した，関数データならではの分析手法もある 関数データ解析手法 ✓ 主成分分析 ✓ 回帰分析 ✓ 仮説検定 ✓ クラスター分析 ✓ 判別分析 ✓ 空間データ解析 ✓ 時系列解析 ✓ 曲線アライメント ✓ 微分方程式モデル

6 • 本資料で紹介する関数主成分分析の手順は次の通り 1. 観測されたデータを関数化し，関数データ集合を得る 2. 関数データ集合に対して関数主成分分析を適用 • 観測時点が個体ごとにまばらにしか観測されていない場合は 本資料の方法ではうまく分析できない場合が多い • 次のスライドから，数理的な手法の説明をします（厳密さに欠けるかも） 「数理は置いといてとりあえず実装してみたい」という方は こちらのページへジャンプしてください 関数主成分分析の手順

7 • 関数データ解析の概要 • データの関数化 • 多変量データに対する主成分分析 • 関数主成分分析 • Rによる関数主成分分析の実装 目次

8 • 観測されるデータは（当然ながら）関数として得られる 訳ではなく，離散時点で数値として得られる • そのため，はじめに各個体に対してデータを関数化する ための(前)処理が必要になる • 関数データ解析に関する手法について説明する前に データの関数化法について説明 離散時点観測から関数データへ

9 • 経時測定データを，次のように表現する 𝑡𝑖𝛼 , 𝑥𝑖𝛼 ; 𝑖 = 1, … , 𝑛, 𝛼 = 1, … , 𝑁𝑖 • 𝑖：観測個体番号， 𝛼：時点番号 • 𝑡𝑖𝛼 ：時点, 𝑥𝑖𝛼 ：観測値 • 𝑛：サンプルサイズ，𝑁𝑖 ：時点数 • 時点 𝑡𝑖𝛼 やその数は，観測 𝑖 ごとに 異なっていてもよい • 右の例に当てはめると次の通り： 𝑖：都市番号， 𝛼：月番号， 𝑛：都市の数， 𝑁𝑖 ：第𝑖番目の都市の観測時点数 𝑡𝑖𝛼 ：月， 𝑥𝑖𝛼 ：気温 離散時点観測データの表記 𝑡14 , 𝑥14 𝑡24 , 𝑥24 𝑡34 , 𝑥34 例：3都市の1年間における 月別平均気温

10 • 観測されているデータは，本来存在する連続的な変動にノイズが伴って 得られていると仮定 • そこで，観測値 𝑥𝑖𝛼 は，次のように時間の未知の関数𝑢𝑖 𝑡 に 誤差𝜀が加わった形で与えられると仮定 𝑥𝑖𝛼 = 𝑢𝑖 𝑡𝑖𝛼 + 𝜀𝑖𝛼 , 𝜀𝑖𝛼 ~ 𝑖.𝑖.𝑑. 𝑁 0, 𝜎2 • これは説明変数を𝑡𝑖𝛼 ，目的変数を𝑥𝑖𝛼 とみなした回帰モデルに対応 • 関数𝑢𝑖 𝑡 を推定することで，関数データ𝑥𝑖 𝑡 ≡ ො 𝑢𝑖 𝑡 を得る • 各観測 𝑖 で，𝑢𝑖 𝑡 をどのように推定すればよいだろうか？ ※本節では以降，添え字 𝑖 を省略して説明します データの関数化

11 • 関数𝑢 𝑡 は，基底関数とよばれる既知の関数系の線形結合で 表されると仮定 𝑢 𝑡 = 𝑤1 𝜙1 𝑡 + ⋯ + 𝑤𝑚 𝜙𝑚 𝑡 = 𝒘𝑇𝝓 𝑡 𝒘 = 𝑤1 , … , 𝑤𝑚 𝑇, 𝝓 𝑡 = 𝜙1 𝑡 , … , 𝜙𝑚 𝑡 𝑇 • 𝒘は未知のパラメータ， 𝝓 𝑡 は基底関数からなるベクトル • 𝝓 𝑡 の種類としてはさまざまなものが考えられている • 多項式 • Fourier級数 • B-スプライン • 動径基底関数 • ウェーブレット • etc. 基底関数展開 このあたりの詳細を 勉強したい方向けの本 Green & Silverman (1993) 小西（2010）

12 • 特定の区間における多項式がその 端点（節点という）で，隣り合う 区間の多項式と接続するように 構築された関数はスプライン関数とよばれる • スプラインの1種であるB-スプライン(𝑟次)は 次の式で逐次的に構築される (de Boorのアルゴリズム) 基底関数の一種：B-スプライン 𝜙𝑗 𝑡; 0 = ൝ 1 𝑠𝑗 ≤ 𝑡 < 𝑠𝑗+1 0 その他 (𝑠𝑗 ：節点) 𝜙𝑗 𝑡; 𝑟 = 𝑡 − 𝑠𝑗 𝑠𝑗+𝑟 − 𝑠𝑗 𝜙𝑗 𝑡; 𝑟 − 1 + 𝑠𝑗+𝑟+1 − 𝑡 𝑠𝑗+𝑟+1 − 𝑠𝑗+1 𝜙𝑗+1 𝑡; 𝑟 − 1 B-スプライン基底関数 (𝑚 = 9) (0次) (1次) (2次) (赤の縦破線が節点の位置)

13 • 基底関数展開に基づく回帰モデル 𝑥𝛼 = 𝒘𝑇𝝓 𝑡𝛼 + 𝜀𝛼 , 𝜀𝛼 ~ 𝑖.𝑖.𝑑. 𝑁 0, 𝜎2 の回帰係数𝒘を推定することで得られる関数 𝑥 𝑡 = ො 𝑢 𝑡 = ෝ 𝒘𝑇𝝓 𝑡 を，関数データとして扱う • 右図のように，基底関数𝜙1 𝑡 , … , . 𝜙𝑀 𝑡 それぞれ に重み 𝑤𝑚 をつけて足し合わせることで 関数を表現している • 回帰係数𝒘は， 𝑥𝛼 を目的変数， 𝝓 𝑡𝛼 を説明変数 とした一般的な回帰モデルの枠組みで，最小二乗 法や正則化法などを用いて推定される 曲線（関数データ）の構築 0 0.2 0.4 0.6 0.8 1 0 20 40 60 -140 -90 -40 10 0 20 40 60 -140 -90 -40 10 0 20 40 60 基底 関数 重みを つける 足し 合わせ

14 • オートバイ衝突実験データ(R package “MASS”)に対して，B-スプラインに 基づく曲線当てはめを適用 • 基底関数𝜙𝑗 𝑡 の個数𝑚によって推定結果が大きく変わる • 基底関数の個数𝑚が多いと観測されたデータに過度に当てはまり（過適合） 逆に少ないとデータの構造を捉えきれない （基底関数の個数𝑚の決定方法については，小西(2010)等を参照） 例：B-スプラインによる曲線当てはめ 基底関数の数：多 基底関数の数：少 基底関数の数：中

15 • 関数データ解析の概要 • データの関数化 • 多変量データに対する主成分分析 • 関数主成分分析 • Rによる関数主成分分析の実装 目次

16 • 多変量（多次元）データの持つ情報を可能な限り失わずに 低次元に圧縮するための方法 • 多変量データを２次元・３次元に圧縮できればデータ全体の特徴を 視覚化できるため，そのままでは全体像を掴みにくい多変量データの 解釈がしやすくなったり，特徴を捉えやすくなる • 圧縮されたデータを利用して元のデータを復元することで 元のデータのノイズを除去できる場合がある • データを圧縮するために，多変量データを1次元の軸に射影 • データの持つ情報＝データの分散と考え，射影されたデータの 分散ができるだけ大きくなる軸を探す そもそも主成分分析って？

17 • 𝑥𝑖1 , … , 𝑥𝑖𝑝 , 𝑖 = 1, … , 𝑛 ：𝑝変量の観測値 • 例：𝑛 = 30名, 𝑝 = 5科目の試験の得点 • 右図の例では𝑝 = 2を想定（横軸： 𝑥𝑖1 ，縦軸： 𝑥𝑖2 ） • 各変量で中心化されていると仮定 （前処理として行っておく） • 多変量データを1次元の軸に射影するということは，𝑥𝑖1 , … , 𝑥𝑖𝑝 の線形結合で 与えられる新たな変数 𝑧𝑖 𝑖 = 1, … , 𝑛 を求めることに対応する 𝑧𝑖 = 𝑢1 𝑥𝑖1 + ⋯ + 𝑢𝑝 𝑥𝑖𝑝 = 𝒖𝑇𝒙𝑖 , 𝑖 = 1, … , 𝑛 (𝒖 = 𝑢1 , … , 𝑢𝑝 𝑇：重みベクトル，𝑥 = 𝑥𝑖1 , … , 𝑥𝑖𝑝 𝑇：観測値ベクトル) • 𝑧𝑖 は𝒖と𝒙𝑖 との内積とみなせる： 𝑧𝑖 = 𝒖, 𝒙𝑖 • 𝑧1 , … , 𝑧𝑛 の分散が最大となるように，重み 𝒖 を求める（ 𝒖 = 1 の条件下で） （関数データの前に）多変量データに対する主成分分析 例：𝑧𝑖 = 1 ⋅ 𝑥𝑖1 + 0 ⋅ 𝑥𝑖2

18 • 𝒖 = 𝑢1 , … , 𝑢𝑝 𝑇 , 𝑋 = 𝑥𝑖𝑗 𝑖𝑗 , 𝒛 = 𝑧1 , … , 𝑧𝑛 𝑇とおくと， 𝑧1 , … , 𝑧𝑛 の分散は 次で表される 1 𝑛 ෍ 𝑖=1 𝑛 𝑧𝑖 2 = 1 𝑛 𝒖𝑇𝑋𝑇𝑋𝒖 = 𝒖𝑇𝑆𝒖 𝑆 = 1 𝑛 𝑋𝑇𝑋 • これを 𝒖 = 1の条件下で最大化する問題は，分散共分散行列𝑆の最大固有値 に対応する大きさ1の固有ベクトルを求める問題に帰着される（証明は割愛） 𝑆𝒖 = 𝜆𝒖 𝒖𝑇𝒖 = 1 • 最大固有値に対応する固有ベクトル 𝒖1 を 第１主成分重み，そのとき与えられる値 𝑧𝑖1 を第１主成分得点という 主成分の導出（１） 第１主成分

19 • 続いて，第１主成分の方向に直交するものの中で，𝑧1 , … , 𝑧𝑛 の分散を 最大にする方向 𝒖 = 𝑢1 , … , 𝑢𝑝 𝑇（ 𝒖 = 1の条件下で）を求める • この問題は，分散共分散行列𝑆の２番目に大きい固有値に対応する 大きさ１の固有ベクトルを求める問題に帰着される 𝑆𝒖 = 𝜆𝒖 𝒖1 𝑇𝒖 = 0, 𝒖𝑇𝒖 = 1 • ２番目に大きい固有値に対応する固有ベクトル𝒖2 を第２主成分重み， そのとき与えられる値 𝑧𝑖2 を第２主成分得点という • 変数の数が 𝑝 = 2 の場合は第２主成分が 最後の主成分になる 主成分の導出（２） 第１主成分 第２主成分

20 • 変数の数が 𝑝 の場合は第 𝑝 主成分まで求められる 第𝑝主成分は，第 1, … , 𝑝 − 1までのすべての重み 𝒖1 , … , 𝒖𝑝−1 と 直交するという条件の下で，𝑧1 , … , 𝑧𝑛 の分散を最大にする𝒖を求める • この問題は，分散共分散行列𝑆の 𝑝 番目に大きい（＝最小の）固有値に対応す る大きさ１の固有ベクトルを求める問題に帰着される 𝑆𝒖 = 𝜆𝒖 𝒖1 𝑇𝒖 = ⋯ = 𝒖𝑝−1 𝑇 𝒖 = 0, 𝒖𝑇𝒖 = 1 • 𝑝 番目に大きい固有値に対応する固有ベクトル𝒖𝑝 を第 𝒑 主成分重み， そのとき与えられる値𝑧𝑖𝑝 を第 𝒑 主成分得点という • 結局，第 𝑘 主成分を求める問題は，データの分散共分散行列の 𝑘 番目に大きい固有値に対応する固有ベクトルを求めればよい 主成分の導出（３）

21 • 30名の５科目の得点データに対して主成分分析を適用（𝑛 = 30, 𝑝 = 5） • 第1，第2主成分得点 𝒛 𝑘 = 𝒖𝑘 𝑇𝑋 はそれぞれ次の式で表される 𝒛 1 = 0.58英語 + 0.56国語 + 0.17理科 + 0.54社会 + 0.18数学 𝒛 2 = −0.13英語 − 0.16国語 + 0.68理科 − 0.54社会 + 0.69数学 • 主成分重みベクトル 𝒖1 , 𝒖2 の値から 第１主成分は文系寄りの能力を， 第２主成分は理系寄りの能力を表している 主成分分析の適用例 名前 英語 国語 理科 社会 数学 A 87 90 92 93 85 B 52 53 51 58 41 C 69 69 96 62 100 D 99 96 59 96 53 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ５科目の得点のデータ 第１主成分得点と第２主成分得点の散布図 理系科目のみ 得意なグループ 全ての科目が 平均的なグループ 全ての科目で 好成績なグループ 文系科目のみ 得意なグループ 全ての科目で 成績がよくない グループ

22 • 関数データ解析の概要 • データの関数化 • 多変量データに対する主成分分析 • 関数主成分分析 • Rによる関数主成分分析の実装 目次

23 • 関数データに対する主成分分析は，多変量データに対する主成分分析の 自然な拡張として与えられる • 関数データ 𝑥𝑖 𝑡 , 𝑖 = 1, … , 𝑛 と重み関数𝑢 𝑡 との（関数空間上の）内積 𝑧𝑖 を考える（多変量データのときは 𝑥𝑖 も 𝑢 もベクトルだった） 𝑧𝑖 = 𝑢, 𝑥𝑖 𝐿2 = න𝑢 𝑡 𝑥𝑖 𝑡 𝑑𝑡 , 𝑖 = 1, … , 𝑛 • この 𝑧1 , … , 𝑧𝑛 の分散が最大となるような（𝐿2 ノルムが1の）重み関数𝑢 𝑡 を考える 1 𝑛 ෍ 𝑖=1 𝑛 𝑧𝑖 2 = 1 𝑛 ෍ 𝑖=1 𝑛 න𝑢 𝑡 𝑥𝑖 𝑡 𝑑𝑡 2 subject to න 𝑢 𝑡 2𝑑𝑡 = 1 • 関数主成分分析は，𝑥𝑖 𝑡 , … , 𝑥𝑛 𝑡 の情報をできるだけ保持する𝑧1 , … , 𝑧𝑛 を 求めるために用いられる 多変量データから関数データへ（１）

24 1 𝑛 ෍ 𝑖=1 𝑛 𝑧𝑖 2 = 1 𝑛 ෍ 𝑖=1 𝑛 න𝑢 𝑡 𝑥𝑖 𝑡 𝑑𝑡 2 subject to න 𝑢 𝑡 2𝑑𝑡 = 1 • これが最大となるような重み関数𝑢1 𝑡 を第1関数主成分 (1st Functional Principal Component; FPC)の重み関数，そのときの 𝑧𝑖 の値𝑧𝑖1 を 第1関数主成分得点という • 多変量データに対する主成分分析と同様，直交性の条件 𝑢𝑗 , 𝑢𝑘 = න𝑢𝑗 𝑡 𝑢𝑘 𝑡 𝑑𝑡 = 0 , 𝑗 = 1, … , 𝑘 − 1 の下で𝑧1 , … , 𝑧𝑛 の分散を最大にするような𝑢 𝑡 が第𝒌関数主成分の重み関数， そのときの𝑧𝑖 の値𝑧𝑖𝑘 が第𝒌関数主成分得点である 多変量データから関数データへ（２）

25 • 多変量データではデータの分散共分散行列の固有値を求めることで 主成分を求めたが，関数データでも同様に考える • 関数データ𝑥1 𝑡 , … , 𝑥𝑛 𝑡 の標本共分散関数を考える 𝑐𝑋 𝑠, 𝑡 = 1 𝑛 ෍ 𝑖=1 𝑛 𝑥𝑖 𝑠 𝑥𝑖 𝑡 • また，標本共分散関数から構成される共分散作用素𝐶を考える 𝐶𝑢 = න𝑐𝑋 ⋅ , 𝑡 𝑢 𝑡 𝑑𝑡 関数主成分の導出（１）：共分散関数最大化

26 • 合成変数𝑧𝑖 の最大化問題を改めて考える 1 𝑛 ෍ 𝑖=1 𝑛 𝑧𝑖 2 = 1 𝑛 ෍ 𝑖=1 𝑛 න𝑢 𝑡 𝑥𝑖 𝑡 𝑑𝑡 2 subject to න 𝑢 𝑡 2𝑑𝑡 = 1 • 目的関数を変形すると 1 𝑛 ෍ 𝑖=1 𝑛 𝑧𝑖 2 = 1 𝑛 ෍ 𝑖=1 𝑛 න𝑢 𝑠 𝑥𝑖 𝑠 𝑑𝑠 න𝑢 𝑡 𝑥𝑖 𝑡 𝑑𝑡 = න 𝑢 𝑡 න 1 𝑛 ෍ 𝑖=1 𝑛 𝑥𝑖 𝑠 𝑥𝑖 𝑡 𝑢 𝑠 𝑑𝑠 𝑑𝑡 = න𝑢 𝑡 න𝑐𝑋 𝑠, 𝑡 𝑢 𝑠 𝑑𝑠 𝑑𝑡 = න𝑢 𝑡 𝐶𝑢 𝑡 𝑑𝑡 = 𝑢, 𝐶𝑢 関数主成分の導出（２）

27 • 重み関数𝑢 𝑡 に対して，次を満たすスカラー𝜆を考える න𝑐𝑋 𝑠, 𝑡 𝑢 𝑡 𝑑𝑡 = 𝜆𝑢 𝑠 共分散作用素を用いて，上式は次で表される 𝐶𝑢 = 𝜆𝑢 • 共分散作用素 𝐶 に対して，𝐶𝑢 = 𝜆𝑢を満たす𝑢 𝑡 は固有関数とよばれる （𝜆は対応する固有値） • 前頁より，分散1 𝑛 σ 𝑖=1 𝑛 𝑧𝑖 2の最大化問題は，内積 𝑢, 𝐶𝑢 の最大化問題に対応 制約条件׬ 𝑢 𝑡 2𝑑𝑡 = 1と併せると，関数主成分は共分散作用素𝐶の 固有値と（大きさ１の）固有関数を求めることに対応する • 最大固有値𝜆1 に対応する固有関数𝑢1 𝑡 は，第1関数主成分の重み関数に対応 • 同様に，𝑘番目に大きい固有値 𝜆𝑘 に対応する固有関数𝑢𝑘 𝑡 は 第𝑘関数主成分の重み関数に対応する 関数主成分の導出（３）

28 • 共分散作用素を計算機上で直接評価することは困難（積分が含まれるため） なため，例えば次の方法が用いられる ① 基底関数展開の仮定 ② 関数化したものを細かい時点で離散化し，多変量データに対して 主成分分析実行 ③ 共分散関数を推定し，これに対して主成分計算 • 本資料では①の方法について説明 関数主成分の計算方法

29 • 共分散作用素からなる方程式 න𝑐𝑋 𝑠, 𝑡 𝑢 𝑡 𝑑𝑡 = 𝜆𝑢 𝑠 を計算するために，関数データ𝑥𝑖 𝑡 と重み関数𝑢(𝑡)に対して 基底関数展開の仮定をおく 𝑥𝑖 𝑡 = ෍ 𝑚=1 𝑀 𝑤𝑖𝑚 𝜙𝑚 𝑡 = 𝒘𝑖 𝑇𝝓 𝑡 , 𝑢 𝑡 = ෍ 𝑚=1 𝑀 𝑏𝑚 𝜙𝑚 𝑡 = 𝒃𝑇𝝓 𝑡 𝒘𝑖 = 𝑤𝑖1 , … , 𝑤𝑖𝑀 𝑇, 𝒃 = 𝑏1 , … , 𝑏𝑀 𝑇, 𝝓 𝑡 = 𝜙1 𝑡 , … , 𝜙𝑀 𝑡 𝑇 • ベクトル 𝒘𝑖 は，関数化の段階で推定されているもの（既知）とする • ベクトル 𝒃 は未知とする これを求めることで，重み関数𝑢(𝑡)を求める 関数主成分の計算（１）

30 • 基底関数展開の仮定の下で，共分散関数は次で表される 𝑐𝑋 𝑠, 𝑡 = 1 𝑛 ෍ 𝑖=1 𝑛 𝑥𝑖 𝑠 𝑥𝑖 𝑡 = 1 𝑛 ෍ 𝑖=1 𝑛 𝝓 𝑠 𝑇𝒘𝑖 𝒘𝑖 𝑇𝝓 𝑡 = 1 𝑛 𝝓 𝑠 𝑇𝑊𝑇𝑊𝝓 𝑡 𝑊 = 𝒘1 , … , 𝒘𝑛 𝑇 • これより න𝑐𝑋 𝑠, 𝑡 𝑢 𝑡 𝑑𝑡 = 1 𝑛 න𝝓 𝑠 𝑇𝑊𝑇𝑊𝝓 𝑡 ⋅ 𝝓 𝑡 𝑇𝒃 𝑑𝑡 = 1 𝑛 𝝓 𝑠 𝑇𝑊𝑇𝑊Φ𝒃 Φ = න𝝓 𝑡 𝝓 𝑡 𝑇𝑑𝑡 関数主成分の計算（２）

31 • この結果から， න𝑐𝑋 𝑠, 𝑡 𝑢 𝑡 𝑑𝑡 = 𝜆𝑢 𝑠 subject to න𝑢2 𝑠 𝑑𝑠 = 1 は次のように変形される 1 𝑛 𝝓 𝑠 𝑇𝑊𝑇𝑊Φ𝒃 = 𝜆𝝓 𝑠 𝑇𝒃 subject to 𝒃𝑇Φ𝒃 = 1 • これが任意の 𝑠 について成り立つため， 1 𝑛 𝑊𝑇𝑊Φ𝒃 = 𝜆𝒃 subject to 𝒃𝑇Φ𝒃 = 1 • ෩ 𝒃 = Φ 1 2𝒃 Φ 1 2はΦ = Φ 1 2Φ 1 2をみたすもの とおけば， 1 𝑛 Φ 1 2𝑊𝑇𝑊Φ 1 2෩ 𝒃 = 𝜆෩ 𝒃 subject to ෩ 𝒃𝑇෩ 𝒃 = 1 となるため， 1 𝑛 Φ 1 2𝑊𝑇𝑊Φ 1 2に対する固有値問題を解けば，固有ベクトルから 係数ベクトル𝒃，ひいては主成分重み関数 𝑢 𝑡 が得られる 関数主成分の計算（３）

32 • 関数データ𝑥𝑖 𝑡 は関数主成分得点𝑧𝑖𝑘 ，重み関数𝑢𝑘 𝑡 を用いて 𝐾 今回の場合は𝐾 ≤ 𝑀 個までの線形和 𝑥𝑖 𝑡 ≈ ෍ 𝑘=1 𝐾 𝑧𝑖𝑘 𝑢𝑘 𝑡 で近似できる • この近似は基底関数展開と同様の形で表されているが 基底関数に対応する𝑢𝑘 𝑡 をデータドリブンで構成するため 一般的に汎用的な基底関数よりも少ない数で精度よく関数を近似できる （このような表現はKarhunen-Loéve (KL)展開とよばれる） • 関数主成分の個数𝐾は，第1~𝑀主成分の固有値の和に対する第1~𝐾主成分まで の固有値の和の比（累積寄与率）の値に応じて決める （明確な決まりはないが，80%とか90%とかが多い） 関数データの復元

33 • 関数データ解析の概要 • データの関数化 • 多変量データに対する主成分分析 • 関数主成分分析 • Rによる関数主成分分析の実装 目次

34 • 成長データに対して 関数主成分分析を 実行するRコード • 各行での処理を 次頁から紹介 関数主成分分析の実行コード（全体） library(fda) # 成長データplot par(mar=c(4,4,1,1), mfrow=c(1,1)) matplot(growth$age, growth$hgtm, pch=16) # 成長データ関数化 growthbasis = create.bspline.basis(c(1,18), 10) growthSmooth = smooth.basis(growth$age, growth$hgtm, growthbasis) growthfd = growthSmooth$fd plot(growthfd) # 関数主成分計算 data_fd_pca = pca.fd(growthfd, nharm=4) plot(data_fd_pca$harmonics) data_fd_pca$varprop #寄与率 # 第1，第2主成分得plot plot(data_fd_pca$scores[, 1], data_fd_pca$scores[, 2], pch=16, xlab="PC1", ylab="PC2", main="PC scores for growth data") abline(h=0, v=0, lty=2)

35 関数主成分分析の実行コード（詳細１） library(fda) # 成長データplot par(mar=c(4,4,1,1), mfrow=c(1,1)) matplot(growth$age, growth$hgtm, pch=16) 関数データ解析パッケージ”fda”読み込み （ついでに成長データも） 横軸に年齢(growth$age)， 縦軸に39名の子供の身長(growth$hgtm)の散布図を描画 時点数はすべての子供で共通 年齢は等間隔になっていないことに注意 子供 → 年 齢 ↓ growth$hgtm の中身 growth$hgtm の散布図

36 関数主成分分析の実行コード（詳細２） # 成長データ関数化 growthbasis = create.bspline.basis(c(1,18), 10) growthSmooth = smooth.basis(growth$age, growth$hgtm, growthbasis) growthfd = growthSmooth$fd plot(growthfd) 横軸に対応する区間[1, 18]で 基底関数の個数10の B-spline基底関数を構築 時点をgrowth$ageとして 身長growth$hgtmを B-spline基底関数により関数化 関数データオブジェクト(growthSmooth)のうち 関数化データの情報（fd）を抽出して 散布図に描画 growthfdの散布図

37 関数主成分分析の実行コード（詳細３） # 関数主成分計算 data_fd_pca = pca.fd(growthfd, nharm=4) plot(data_fd_pca$harmonics) data_fd_pca$varprop #寄与率 # 第1，第2主成分得点plot plot(data_fd_pca$scores[, 1], data_fd_pca$scores[, 2], pch=16, xlab="PC1", ylab="PC2", main="PC scores for growth data") abline(h=0, v=0, lty=2) 主成分の個数(nharm)を𝐾 = 4として growthfdに対して関数主成分分析実行 第4関数主成分までの固有関数𝑢𝑘 𝑡 (harmonics)を描画 第4関数主成分までの寄与率(varprop)を出力 第1，第2関数主成分得点𝑧𝑖1 , 𝑧𝑖2 (scores)を描画 data_fd_pca$harmonics（左） data_fd_pca$scores（右） の散布図 黒，赤，緑，青の順に 第1, 2, 3, 4主成分

38 • 第1関数主成分重み関数は 全期間（1~18歳）を通じて 正の値をとっている → 全期間を通じて身長が高い人ほど 第1主成分得点の値は大きい （逆に低い人ほど小さい） • 第2関数主成分重み関数は 12歳～16歳で急激に負の方向へ減少 → 第2次性徴が遅れて現れた人ほど 第2主成分得点の値は大きい （早い人ほど小さい） • 関数主成分重み関数から，関数データ の傾向を2次元に集約できる 分析結果の考察 ――：第１関数主成分 - - -：第２関数主成分 子どもの身長のデータと 関数主成分重み関数 主成分得点plot （横：第１，縦：第２）

39 • 関数主成分分析は関数データ解析手法の一つで，各観測個体が経時的に 観測されたとき，それを関数化したものの特徴を低次元のデータに 集約するために用いられる • 関数主成分分析を適用するためにはデータを関数化する必要がある ここでは関数化の方法として基底関数展開を使用 • 関数主成分分析は多変量データに対する主成分分析において ベクトルとして扱っていたものを関数として扱う 参考文献： – Ramsay, J. O. and Silverman, B. W. (2005) Functional Data Analysis 2nd ed. Springer. – R package “fda” https://cran.r-project.org/web/packages/fda/index.html – 小西貞則 (2010). 多変量解析入門．朝倉書店 – Green, P. and Silverman, B. W. (1993). Nonparametric Regression and Generalized Linear Models. Chapman and Hall/CRC. まとめ