多変量データ解析と時系列解析 吉田 亮 (統計数理研究所; モデリング研究系) e-mail: [email protected] スライドのダウンロード： http://daweb.ism.ac.jp/~yoshidar/index_j.htm 1 １日目： 記述統計、確率分布 ２日目： 推定・検定 ３日目： 回帰分析 ４日目： 多変量解析、時系列解析 10:00～11:10 主成分分析 11:20～12:30 判別分析 13:30～14:40 クラスタ分析 14:50～16:00 時系列解析 前日までの講義内容と比較して難易度はやや高めです。 統計学の実践をデモンストレーションすることで、統計学がどのように使われてい るかを実感して頂きたい。 フォローしきれない人は、各々の要素技術をどのような局面で使うべきか、使うこと でどのような情報抽出ができるか、最低限これらのイメージを持って帰って下さい。 解析ツールの単なるユーザーになるのでなく、ゆくゆくは解析的なメカニズムを きっちりと理解し、自由自在に統計学を操れるデータ解析スペシャリストに 2 多変量データ解析の例 サンプル（個体）が複数の変数（特徴変数）によって特徴付けられるデータ (例) 患者（個体）400人の約20,000個の遺伝子の発現量（変数） p = 20,000 個の特徴変数 1 2 p 1 x1,1 x1,2 x1, p 2 x2 ,1 x2 , 2 x2 , p n xn ,1 xn ,1 xn , p n = 400人の個体 3 多変量データの直観的理解は困難 データファイルを開き、数万個の遺伝子の発現値を一つずつ眺めてみ ましょう。並はずれた根気強さを持つ人にしかできません。 ちなみに古いエクセルでは、 65,536行、256列におさまるデータしか開 けません。 データをヒートマップ表示してみよう。 個体: n = 400 変数: p = 20,000 データ量が膨大で何も見えてこない 個体と変数をそれぞれパターンに基づき分類して、データを並べ替える 4 多変量データ解析技術は実践で間違いなく役に立つ 人間の脳は多変量データを直観・統合的に理解することが苦手 多変量解析は、大規模データを処理する上で、脳の情報処理機能の代替として使う そもそも、多くの局面において、目で見て直観的に理解できる「単変数」のデータに対 して、わざわざ統計解析を行う必要がない？ 本講義では「初等的な」多変量解析の手法を概観 本講義で取り扱う多変量解析の手法 • 主成分分析 • 判別分析（教師つき分類） • クラスタ分析（教師なし分類） • 多重回帰分析 データの特徴を捉える データを分類する 多変数データの関係を捉える 5 準備: 用語と記号 特徴変数ベクトル（ i 番目のサンプルを特徴付ける p 個の変数 ） xi x1,i , x2,i , , x p 1,i , x p ,i T 転置記号 p 個の特徴変数が並んでいる 添え字 i は個体番号を表す ※ xi = ( 個体 i の体脂肪率, 個体 i の血糖値, 個体 i の所得 ) 変数の数は p で表す。 総サンプル数は n で表す。 xi i 1, , n （サンプル）平均ベクトル（各変数の平均値を p 個並べたもの） x x1 , x2 , , x p 1 , x p T ベクトルの要素は各変数の平均値 1 n x j x j ,i n i 1 j 1, , p 6 準備: 用語と記号 （サンプル）分散共分散行列 • p × p 個の変数ペアの共分散（共変動性を表す尺度）と p 個の変数の分散（ばらつき 尺度）を要素に持つ p × p の対称行列 • 単変数データの分散の概念を多変数に拡張したもの • 多変数データの変動特性を表す最も基本的な統計量 対角要素には分散 2 1 n x x 1,i 1 n i 1 1 n x x2 x1,i x1 S n i 1 2 ,i n 1 n x p ,i x p x1,i x1 i 1 非対角要素 ( i, j ) には変数 i と j の共分散 1 n x x2 x1,i x1 n i 1 2 ,i 2 1 n x x 2 n i 1 2 ,i 1 n x xp n i 1 p ,i x p 2 ,i x2 x1 n 1 x x x x p 2 ,i 2 n i 1 p ,i 2 1 n x xp n i 1 p ,i x 1 n x xp n i 1 p ,i 1, i 7 分散共分散行列の直観的理解（2変数の場合） 変数1の分散 変数1と 変数2の共分散 S 変数 1 と 変数2の共分散 変数 2 の分散 変数1と2の間に正の相関がある場合 変数1と2の間に負の相関がある場合 8 準備: 用語と記号 （サンプル）相関行列 • p × p 個の変数ペアの相関係数を要素に持つ p × p の対称行列 • 対角要素は1 1 1 n x2 ,i x2 x1,i x1 n s s R 1 2 i 1 n 1 x p ,i x p x1,i x1 n s s i 1 1 p sj 1 n x xj n i 1 j ,i 1 n s1 s2 x n i 1 2,i x2 x1,i x1 1 n s p s2 1 1 n s p s2 x n i 1 p ,i xp n 1 x xp n s p s1 i 1 p ,i x 2 ,i x2 x n i 1 p ,i xp 1 x 1, i 1 x x 2,i 2 x 標準偏差 2 9 主成分分析 （ PCA: Principal Component Analysis ） 10 主成分分析で何ができるか？ 相関のある多変数データが保持する「情報」をできるだけ失うことなく、少数の変数 （主成分）に要約することが目的 (データの次元圧縮、特徴抽出) 主成分分析が対象とする情報 = 分散や変数間の相関関係（共分散） 4変数、150個体のデータ分布。 150個体は3種類 のグループに分類されている。 PCAを適用して2次元の合成変数（主成分）に圧縮。元 のデータの分布特性（グループ構造）が保存されている。 11 (例) 手書き文字画像データの特徴抽出 100枚の手書き文字画像 ：256変数(ピクセル) 復元された100枚の手書きデータ 次元圧縮 復元 PCAを適用して20パターンの特徴を抽出（固有ベクトル） 12 PCAの基本的な考え方 数学 (x1,i) と英語 (x2,i) の50人 (n=50) の得点からなる2変数デ-タを重み係 数 w = (w1, w2) を使って1変数に要約（射影）する。 2変数データを1変数（主成分）に変換 yi w1 x1,i w2 x2,i wT xi i 1, , n 重み係数 w = (w1, w2) をどのように決めるか？ データの変動特性をできるだけ保存するように！ 13 PCAの設計原理 p 個の変数の場合も同様に yi w1 x1,i w2 x2,i w p x p,i wT xi i 1, , n PCAは yi の分散が最大になるような w を求める。 1 n 2 max yi y w n i 1 14 分散を最大にする射影の方向、直観的理解 n 個のデータ点を方向 w に射影する w w この例では、w’ に射影する方が得られる合成 変数の分散が大きい。 • 元の２変数データのばらつきをより適切に 反映できる変換。 • 相関の高い変数群を一つの合成変数にま とめる。 15 主成分分析の定式化 ① 計算手続きは、分散共分散行列の固有値問題を解くだけです。 射影したデータ点、すなわち主成分の分散に対して行列表現を行います。 1 n 2 [ 主成分の分散 ] S y w yi n i 1 1 n p p wk wh xk ,i xh,i n i 1 k 1 h 1 ( y の式を代入 ) wT S x w ( 2次形式で表現 ) ※ データの分散共分散行列を方向 w に射影したもの ※相関行列を使っても構いません ここで、Sx は元のデータの分散共分散行列 n 2 x1,i 1 i 1 Sx n n x2,i x1,i i 1 x x 2 , i 1, i i 1 n 2 x2 , i i 1 n 16 主成分分析の定式化 ② 分散の最大化 max S y w w ベクトルの長さに対す る拘束条件が必要 w T Sx w s.t. w 1 ※ ここでは長さを1とするが、任意の値に設定しても良い。 ※ 拘束条件の必要性は Syの定義式に立ち返ればすぐに分かる。 c2 S y c w Sx c w T 係数ベクトルに任意の定数 c を掛けることで、いくらでも分散が大きくすることが できる。したがって、係数ベクトルの長さに対して適当な制約を課す必要がある。 17 主成分分析の定式化 ③ 最適化問題を解いて、係数ベクトルを求めてみよう。 max S y w w s.t. w 1 ラグランジェ未定乗数法を使います。 max w T S x w w T w w ラグランジェ乗数 ラグランジェアンの微分を計算すると次の固有方程式が得られる。 Sx w w 係数ベクトルはデータの分散共分散行列の固有ベクトルのいずれかになる。 ラグランジェ乗数は固有値 18 主成分分析の定式化 ④ 係数ベクトルの「候補」として、m 個の固有ベクトルが得られました。 それらの中でどれを選べば良いか？（分散を最大にするものは？） (解の候補) (固有方程式を満たす) Sx w*j j w*j 固有ベクトル 固有値 j 1, *j j 1, w*j , m , m 最大固有値に相当する固有ベクトルが、分散を最大にする係数ベクトル 固有値の大きさ = 固有ベクトルで射影したデータの分散 左から固有ベクトルを掛ける Sx w*j j w*j w*jT Sx w*j j w*jT w*j j 射影したデータ点の分散 固有値 19 第1主成分、第2主成分、第3主成分 最大固有値の固有ベクトルは、分散を最大にする係数ベクトル これに基づく合成変数を第1主成分と言います y1,i w1,1 x1,i w1,2 x2,i w1, p x p,i w1T xi i 1, , n 第１主成分だけでデータの変動特性を十分説明できないとき、第２主成分、第３主成 分と順に合成変数を追加する。 大きさの順序が2番目以降の固有ベクトルを使って、第2主成分から第m主成分を求 める。 第2主成分 y2,i w2T xi 2 第3主成分 y3,i w3T xi 3 ym,i wmT xi m 第m主成分 固有値(=分散)を降順に並べる 20 第1主成分と第2主成分 w1 , , wm は p 次元空間の正規直交基底 つまり、元データの座標変換に相当する。 第1主成分ベクトル w1 第2主成分ベクトル w2 21 各主成分の解釈の仕方 一般に、p 個の主成分スコアの大きいものから数個の変数を取り出し、各主成分にど のような変数が寄与しているか調べ、各々の科学的知識に基づき主成分を解釈する。 係数の大きさ（主成分スコア）は各主成分に対する寄与率を表す y1,i wmath,1 xmath,i wlang,1 xlang,i 第1主成分は数学の得点を反映 第1主成分は国語の得点を反映 wmath,1 wlang,1 wmath,1 wlang,1 第1主成分ベクトル 第1主成分ベクトル 22 解析例：米国50州、 居住者100,000人当たりの犯罪統計 • 殺人 (Muder) • 暴行 (Assault) • 性犯罪 (Rape) • 都市部居住者人口の割合 変数間で数値データの単位（ス ケール）が異なることに注意 このような場合、分散共分散に基づ くPCAは不適切 23 分散共分散行列の固有ベクトル Murder Assault UrbanPop Rape PC1 0.041704 0.995221 0.046336 0.075156 PC2 -0.04482 -0.05876 0.976857 0.200718 PC3 0.079891 -0.06757 -0.20055 0.974081 PC4 -0.99492 0.038938 -0.05817 0.072325 各主成分は個々の変数の情報だけを反映 相関行列の固有ベクトル Murder Assault UrbanPop Rape PC1 -0.5359 -0.58318 -0.27819 -0.54343 PC2 0.418181 0.187986 -0.87281 -0.16732 PC3 -0.34123 -0.26815 -0.37802 0.817778 PC4 0.649228 -0.74341 0.133878 0.089024 相関の高い変数ペアを合成することができた 第１主成分： Murder + Assault + Rape 第2主成分： UrbanPop + Murder 24 主成分寄与率: 各主成分はデータの全変動どのくらいの割合を占めるか 第1主成分の寄与率 第2主成分の寄与率 第3主成分の寄与率 第4主成分の寄与率 1 0.62 2 0.25 3 0.09 4 0.04 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 第１主成分と第2主成分でおよそ87％の変動を説明できる 残りの成分は、観測誤差と解釈することとする 25 低 都市部への人口集中度 多 犯罪件数 少 高 26 各主成分の解釈の仕方 ② 理化学研究所の鎌谷直之グループ ディレクターは、遺伝子のわずかな 個人差によって、日本人が二つの 集団に大別できることを、理化学研 究所が明らかにした。沖縄の人の 大部分が含まれる「琉球クラスター （集団）」と、本土の人の大部分が属 する「本土クラスター」があるという。 ２５日付の米国人類遺伝学会誌に 発表した。 (毎日新聞 2008年9月26日 東京朝刊) 理化学研究所 ゲノム医科学研究センター 鎌谷直之グループディレクター 27 主成分分析のまとめ 多変数データの次元削減や特徴抽出、視覚的理解に活用できる統計技法 合成変数の分散（ばらつき）が大きくなるように変数に重みを与える（データ の分布特性に関連する変数を自動的に同定することに相当） 数万次元のデータでも、統計解析ソフトウェア（例えば R）を利用すれば簡 単に計算できます。 参考図書 ※ 初学者向け 中村 永友 (著) 『多次元データ解析法 (Rで学ぶデータサイエンス 2) 』 共立出版 ※ 中級者以上 小西 貞則 (著) 『多変量解析入門―線形から非線形へ』 岩波書店 C. M. ビショップ (著) 『パターン認識と機械学習 上・下』 シュプリンガー・ジャパン 29 分類の統計学 判別分析 + クラスタ分析 30 判別分析 ( Discriminant Analysis ) 31 判別分析とは？ ある個体の特徴量（複数の観測値）から、その個体が、あらかじめ与えられたいくつ かの群のどれに属するかを判断したい。例えば、 ・ 血圧やバイオマーカーの検査値から、疾患・非疾患の診断を行う。 ・ 財務データから、企業がデフォルト起こすかどうかを判断する。 Ａさん 甲薬 Ａさんのデータ 効果有？無？ B さん 甲薬 B さんのデータ 効果有？無？ 32 判別ルールを作る 属性クラス G1 p 個の変数で特徴付けられた個体 x x1 , x2 , , xp ・ ・ ・ ・ ・ 判別ルール 属性クラス GK 変数 x1 データ x がこの領域 に入ればG1に分類 判別境界(ルール) G1 G3 G2 変数 x2 33 「訓練データ」を使って判別ルールを作り、「テストデータ」にもとづきルールの 良さを評価する。 現在手元にあるデータを分割 訓練データ ① 「訓練データ」を使って判別ルールを設計 ① 属性ラベルの分かっているデータから判 別境界を学習 テストデータ 未来のデータ ② 性能評価用 ③ 実運用 ② 属性ラベルの分からないデータのクラス を予測し、判別ルールの性能を評価 G1 G3 G2 グループ１ グループ２ グループ３ 34 確率モデルにもとづく判別分析 35 確率モデルを導入する データはある確率分布から生成されていると仮定する。 属性ラベル yi があたえられたもとでのデータ xi の条件付き分布 yi , xi ~ p yi , xi p xi yi p yi i 1, 属性ラベルと特徴変数の同時分布 yi 1, xi ,n 属性ラベル yi の実現確率 , K 個体 i の属性ラベル ラベル yi =1 のデータ生成分布 p xi yi 1 p 個の特徴変数を要素にもつベクトル ラベルは確率的に決まる p yi 1 0.6 p yi 2 0.4 p xi yi 2 ラベル yi =2 のデータ生成分布 36 事後確率にもとづく判別 ベイズの定理 にもとづき，データが各クラスに属する（事後）確率を計算する。 「特徴 xi をもつ個体がラベル yi =k に属する確からしさ」 (ラベルの事後確率) p yi k x i p xi yi k p yi k p xi (例) P ( ラベル=企業の倒産 | 特徴=企業の財務データ ) p yi 1 x i ? p yi 2 x i ? k 1, ,K 観測された財務データ xi の 事後確率を計算する 倒産リスクは？ 37 ベイズルール 事後確率が最も大きい属性に分類する判別方式を「ベイズルール」という。 k* arg max k p yi k xi 直観的で自然な判別方式 • 判別の平均的な誤り率（ベイズリスク）を最小にする判別方式 問題は p yi , x i に対してどうようなモデルをおくのか、また、仮定 されたモデルを訓練データからどうやって推定するのか？ このあと正規分布モデルを例に、実際の判別分析の手続きを導出します。 38 多変量正規分布にもとづく判別 39 多変量正規分布 1次元の正規分布 : 2 平均 分散 １次元正規分布の確率密度関数 x 2 p , 2 x exp 2 2 2 2 1 多変量正規分布 : 平均ベクトル p ,V x 2 p/ 2 V 1/2 分散共分散行列 V T 1 exp x V 1 x 2 2変量正規分布の確率密度関数 高 低 相関 40 最尤法によるパラメータの推定 最尤法（さいゆうほう: maximum likelihood estimation) 確率分布 p x のパラメータ を n 個の標本 x1 , , xn から推定する際、最尤法は尤度 と呼ばれる「標本の確率分布に対する適合度」を最大にするパラメータを推定値とする。 n (対数) 尤度 l log p xi i 1 尤度を最大にするパラメータを求めるのが最尤法 1次元正規分布の平均パラメータの推定（イメージ） 0 の正規分布 高 低 尤度 3 の正規分布 100個の標本 41 多変量正規分布の最尤推定 多変量正規分布の最尤推定量は、 平均ベクトル 分散共分散行列 1 n ˆ xi n i 1 （標本平均） n 1 T Vˆ xi ˆ xi ˆ n i 1 （標本分散共分散行列） 最尤推定量の導出に関する参考資料： • 自然科学の統計学 (基礎統計学), 東京大学出版会 42 線形判別と非線形判別 以下では、多変量正規分布を利用して、線形判別と2次判別を導く。 線形判別 －判別の識別面が１次関数－ 非線形判別 －判別の識別面が非線形関数－ クラス1と判定 クラス2と判定 識別面が2次曲線になる場合を「2次判別」という 43 多変量正規分布を使った線形判別 (モデリング) 属性ラベルと特徴変数の生成過程として、多変量正規分布を仮定する。 p x i , yi k p x i | yi k p y i k ① ② K 面体のサイコロを振って属性 k を決定し、クラス k の正規分布からデータを生成する。 ① 各ラベルに属するデータの生成モデルとして正規分布を仮定する。 p xi yi k 1 2 p/2 V 1/ 2 T 1 exp x k V 1 x k 2 ※ ここで，分散共分散行列はクラス間で共通と仮定 ② 各属性の(事前)確率 p yi k k k 1, ,K 44 線形判別の導出 ① ベイズの定理から、各属性ラベルの事後分布は次のような形であたえられる。 log p yi k xi log p xi yi k p yi k p xi 1 T const . log k xi k V 1 xi k 2 マハラノビス距離 あとは、事後分布が最も大きくなる属性 k に分類すればよい。 「近い」クラスに分類 クラス1の分布 クラス２の分布 xi 1 2 45 線形判別の導出 ② データ xi が二つの属性 k と h のどちらに属するかを判別するには，事後分布の対数比 が０よりも大きいか，あるいは小さいかを評価すればよい． f( k ,h ) xi log p yi k log p x i yi k p yi h p xi yi h 1 T log k k h V 1 k h xi T V 1 k h h 2 w0 xi T w1 識別面は1次関数で表される 0 x1 • 分散共分散行列はクラス間で共通と仮定する ことで、線形の識別面が得られる。 • この仮定を外すと識別面は二次曲線になる。 x２ 46 最尤推定量の計算 n 個の「訓練データ」を使って最尤法でパラメータを推定する。 n l log p yi , xi | 1 , i 1 , K 1 , , K , V そのあとで、推定したパラメータを事後分布の式にプラグインして判別を行う。 混合比率 平均 分散共分散行列 ˆ k nk n ˆ k 1 nk （クラス k に属する個体の割合） xi （クラス k に属する個体の標本平均） i in class k 1 K T ˆ V xi ˆ k xi ˆ k n k 1 i in class k （クラスごとに計算した標本分散共分散行列を さらに平均化したもの） nk : クラス k に属する個体数 47 多変量正規分布にもとづく2次判別 48 2次判別と線形判別 各属性の特徴変数の分布として「共通の分散共分散行列」をもつ多変 量正規分布を仮定することで線形判別を導出した。 次に、「非共通の分散共分散行列」をもつ多変量正規分布を仮定するこ とで2次判別関数を導出する。 識別面は2次関数 識別面は線形 クラス1と2の正規分布が構成 する中心からのマハラノビス 距離は同一と仮定 異なるマハラノビス距離 49 多変量正規分布を使った2次判別 (モデリング) 属性ラベルと特徴変数の生成過程として、多変量正規分布を仮定する。 p x i , yi k p x i | yi k p y i k ① ② ① 各ラベルに属するデータの生成モデルとして正規分布を仮定する。 p x i yi k 1 2 p/ 2 1 1/ 2 k V T 1 exp x k Vk1 x k 2 クラス間で異なる分散共分散行列を仮定する。線形判別との 相違点はこの点のみ。 ② 各属性の(事前)確率 p yi k k k 1, ,K 50 2次判別の導出 ① ベイズの定理から、各属性ラベルの事後分布は次のような形であたえられる。 log p yi k xi log p x i yi k p y i k 変更点はここ p xi 1 1 T const . log k log Vk xi k Vk1 xi k 2 2 マハラノビス距離 あとは、事後分布が最も大きくなる属性 k に分類すればよい。 「近い」クラスに分類 クラス1の分布 クラス２の分布 xi 1 2 51 2次判別の導出 ② データ xi が二つの属性 k と h のどちらに属するかを判別するには，事後分布の対数比 が０よりも大きいか，あるいは小さいかを評価すればよい． log p yi k p yi h log p x i yi k p x i yi h Vh k 1 T T 1 1 log x k Vk xi k xi h Vh xi h log h 2 i Vk w0 xi T w1 xi T Wxi 識別面は2次曲面で表される 0 • 分散共分散行列はクラス間で共通と仮定することで、線形の識別面が得られる。 • この仮定を外すと識別面は二次曲線になる。 52 Fisher (E. Anderson) のアヤメデータ Fisher (1936) が線形判別分析の例題として使用 ３種類のアヤメ －setosa, versicolor, virginica setosa と virginica を併せて１つの群とする 計測した特徴量 ヒオウギアヤメ（setosa） －がく片の長さと幅、花弁の長さと幅 （４変数） 採取されたデータ －各アヤメをそれぞれ５０サンプル収集（計150サンプル） 花片 問題 －これら４つの特徴変数から、アヤメの種類を 判別するための規則（線形判別関数）を構成する。 がく片 53 ① setosa, ② versicolor, ③ virginica 特徴量の分布（箱ひげ図） がく片の長さ ①＆③ vs ② を判別したい 花弁の長さ がく片の幅 Sepal.Length Sepal.Width 花弁の幅 Petal.Width 3 5.5 1.0 3.0 6.0 4 1.5 6.5 3.5 5 7.0 2.0 4.0 6 7.5 2.5 7 8.0 Petal.Length setosa ① 0.5 1 2.0 4.5 2 5.0 2.5 変数をクラス①&③と②の違い？ versicolor virginica setosa versicolor virginica setosa versicolor virginica setosa versicolor virginica ② ③ ① ② ③ ① ② ③ ① ② ③ 54 線形判別の結果 －setosa & versicolor vs virginica－ ２種のアヤメの全データを訓練用とテスト用に分ける。 （ランダムに訓練データを選んで、7:3の割合になるようにデータを分割） 訓練データから得られた 線形判別関数 0.1x1 2.14x2 1.06x3 2.37x4 5.75 （ x1 , x2 : がく片の長さと幅、 x3 , x4 : 花弁の長さと幅 ） 誤判別率が高い。特にversicolorの判別性能が悪い。 誤判別の個数 PREDICTION 訓練データ の判別結果 TRUE (訓練エラー 0.3) setosa & virginica versicolor setosa & virginica 58 12 versicolor 19 16 PREDICTION テストデータ の判別結果 (予測エラー 0.2) TRUE setosa & virginica versicolor setosa & virginica 28 2 versicolor 7 8 55 2次判別の結果 －setosa & versicolor vs virginica－ ２種のアヤメの全データを訓練用とテスト用に分ける。 （ランダムに訓練データを選んで、7:3の割合になるようにデータを分割） 訓練データから得られた 線形判別関数 0.80x12 3.43x22 5.27x32 32.19x42 （ x1 , x2 : がく片の長さと幅、 x3 , x4 : 花弁の長さと幅 ） 2次判別を使うことで誤判別率が大幅に減少 PREDICTION 訓練データ の判別結果 TRUE (訓練エラー 0.02) setosa & virginica versicolor setosa & virginica 69 1 versicolor 1 34 PREDICTION テストデータ の判別結果 (予測エラー 0.07) TRUE setosa & virginica versicolor setosa & virginica 30 0 versicolor 3 12 56 線形判別と2次判別の結果に関する考察 データを見れば線形判別がうまくいかないことはすぐに理解できます。 「データを見る」ことが重要！ 線形判別では、setosa & virginica （赤と 青）とversicolor（緑）のグループに対して 「共通の分散共分散行列」を仮定する。 アヤメデータの特徴量を対散布図行列で図示したもの しかしながら、2群の分布形状（ばらつき の方向）は明らかに異なる。また、線形の 識別面で分離することは無理がありそう である。 この場合は、データを見ることで、線形判 別ではパフォーマンスが出ないことはす ぐに分かる。 57 判別分析のまとめと補足 判別分析 ― ある個体を、その特徴量から、予め与えられたいくつかの群の どれに属するかを判定したいときに、既に属性の分かっているいくつかの個 体の特徴量（訓練データ）を用いて判別ルールを構成し、そのルールに基 づいて判定を行うこと 正規分布を利用した線形判別と2次判別分析 上記の方法以外にも実に様々な判別分析の手法が整備されている。 • ロジスティック判別 • サポートベクタマシン • アダブースト データが非連続量（例えば、カテゴリ型の特徴変数や化合物の構造など）の場合や、 クラス間の識別面が複雑な形状をしている場合は、本講義で紹介した正規分布によ る判別分析は直接適用できません。 また、分布の推定手法として最尤法を紹介したが、サンプル数が小さかったり、特徴 データの次元が極端に高い場合は、最尤法がうまく機能しないことがあります。 58 より詳しく勉強したい人のための参考図書 小西 貞則 (著) 『多変量解析入門―線形から非線形へ』 岩波書店 C. M. ビショップ (著) 『パターン認識と機械学習 上・下』 シュプリンガー・ジャパン 金森 敬文, 竹之内 高志, 村田 昇 (著) 『パターン認識 (Rで学ぶデータサイエンス 5)』 共立出版 59 クラスタ分析（ Cluster Analysis, Unsupervised Learning ） 60 データのかたまり（クラスタ）を見つける統計学 データのパターンに応じてサンプルを分類するための統計技術 同一クラスターに分類されるサンプルは類似性が高い集団 クラスタ分析 与えられるのは、各個体の特徴データ（多変量データ）のみで、それらを用いて各 個体をいくつかの群（クラスタ）に分類する。 教師なし学習(unsupervised learning) 判別分析 あらかじめ与えられた複数の群に対して、どの群に属するか既に分かっている各個 体の特徴データを用いて、判別法を構成する。 教師あり学習(supervised learning) 61 教師付き分類（判別分析） 今手元にあるデータの属性ラベル（教師）と特徴量のパターンか ら判別ルールを作り、未来のデータの属性ラベルを予測する 変数1 正常細胞 A癌細胞 B癌細胞 訓練データ 将来のデータ（テストデータ） 変数2 62 教師なし分類（クラスタ分析） 属性ラベルの情報を利用せずに分類のルールを作る サンプルの属性に関する情報がない 一部の属性ラベルが欠損 敢えて利用しない（仮説の検証） 変数1 変数2 63 大量データの視覚的理解や潜在クラスの発見 階層的クラスリングを適用して、イースト菌の6000遺伝子の発現パターンを分類 関連する遺伝子群の発見？ 64 多変量データの直観的理解は困難 エクセルでデータファイルを開いて、数万個の遺伝子の発現値を一つ ずつ眺めてみましょう。超人的な根気強さを持つ人にしかできません。 ちなみに古いエクセルでは、 65,536行、256列におさまるデータしか開 けません。 まずはデータをヒートマップ表示してみましょう。 ? 個体: n = 400 変数: p = 20,000 データ量が多すぎて何も分からない。 このような状況で多変量解析は役立ちます。 個体と変数をそれぞれパターンに応じて分類し、データを並べ替える 65 データの類似度 66 特徴変数ベクトルと個体間の類似度 各個体に対して p 個の特徴変数を観測 個体 i 個体 j xi xj 変数B Cluster 2 xj 距離は？ 類似度 Cluster 1 xi 変数A 類似度を評価する尺度を どのように設計するかが鍵 67 標準的な類似度： (擬)距離 (a) ユークリッド距離 ※ 距離の公理を満たす必要はない ※ その他、ピアソンの相関係数など d E u, v p u i i 1 vi 2 p (b) City-Block 距離計量 d CB u, v ui vi i 1 (c) 最大距離計量 d MD u, v max ui vi i u 変数 B 最大距離計量 v City-Block 距離計量 変数 A 68 K平均法 ( K-means algorithm ) 69 K平均法の概要 特徴空間上に n 個のサンプル(個体)が配置されている。 K 個のクラスタ中心（平均）を求める。 各サンプルを最も近いクラスタ中心に分類する。 Cluster 1 中心からの距離を測る Du,v d u, v Centroid Cluster 3 Centroid 最も近いクラスタ中心 Centroid Cluster 2 70 K平均法のアルゴリズム Step 0. （初期化）n 個の個体を適当にグルーピング Step 1. (中心の計算) K 個のクラスタ平均を計算 Step 2. (再グルーピング) データ点から最も近いクラスタ平均に割り付ける Step 3. Step 1に戻る (1) 初期化+クラスタ平均 (2) 再グルーピング (3) クラスタ平均の更新 変数B 変数B 変数A 変数B 変数A 変数A クラスター1 クラスター2 71 K平均法は何をしているか？ サンプルを表すインデックス: i 1, 分類関数を求める: C i 1, ,n , K (サンプル i はどのクラスタに属するか？) グループ内距離の総和 W(C) を最小化する学習方法 1 K W C d xi , x j 2 k 1 i:C i k j:C j k ※ グループ内総距離は集団内の同一性を測る尺度 変数B 変数B グループ内距離の総和を より小さくする分類を探索する 変数A 変数A 72 グループ数決定の指針 素朴な方法: GAP 統計量に基づく選択 GAPK W CK 1 W CK • クラスタ数を K 個から K+1 個に増やしたとき、グループ内総距離はどれく らい減少するか？ • クラスタ数を増やしてもGAP統計量が「ほとんど減少しない」最小のクラスタ 数が有効クラスタの目安 グループ内距離 変数B クラスタ数を増やしてもグループ内総 距離はあまり減少しない 変数A 2 3 4 5 クラスタ数 73 階層型クラスタリング ( hierarchical clustering ) 74 階層型クラスタリング 階層状のクラスタを探索・表現 樹形図（tree diagram）によるクラスタ構造の視覚化 米国50州の犯罪統計のクラスタ解析 階層的クラスリングを適用して、イースト菌の6000 遺伝子の発現パターンを分類 2 groups 3 groups 75 樹形図によるクラスタ表現 最も下の階層に n 個の個体が並ぶ（これらをn個のクラスタと考える） 段階的に最も「近い二つのクラスタの組」を選び出し、それらを一つの頂点で結んでいく 頂点の高さはグループ内総距離 米国50州の犯罪統計（殺人、暴行、都市人口、レイプ）のクラスタ解析 大 3 clusters に分割 4 clusters に分割 小 グループ内総距離 サブクラスタに分割 n 個の個体（米国50州） 76 クラスタの類似度 階層型クラスタリングでは、個体間の類似度とは別に、「近隣クラスタ」を定義するため のクラスタ間の類似度が必要 以下、代表的なクラスタ間類似度を挙げる Single linkage Centroid linkage 最も近い二つの個体間の距離 クラスタ中心間の距離 Cluster 1 Cluster 2 Cluster 1 Cluster 2 Complete linkage Average linkage 最も遠い二つの個体間の距離 クラスタに属する個体の全組合わせ から計算される平均類似度 Cluster 1 Cluster 2 Cluster 1 Cluster 2 77 アルゴリズム： 凝集型 (agglomerative) i. 単一の個体からなる n個のクラスタから開始 ii. 最も近い二つのクラスタの組を選び、それらを一つのクラスタに併合する。 iii. 手続き ii をクラスタ数が1つ（全個体を含む）になるまで繰り返す。 他にも分岐型 (divisive) の方法もある。その場合、全個体を含むクラスタから開始して、 逐次的にクラスタを分割しながら階層構造を構築していく。 凝集型アルゴリズムの例 Step 5 Step 4 Step 6 Step 3 Step 1 B A G Step 2 D Step 5 F Step 4 Step 3 Step 2 E Step 1 C A B C D E F G 78 クラスタ分析のまとめと補足 クラスタ分析 ― 特徴量のパターンが類似している集団を見つけ出し、一つのグルー プに分類するための統計手法 K平均法と階層型クラスタリングについて紹介 大量のデータがあたえられたとき、まずはじめにクラスタ解析を行い、クラスタ毎に データを図示すると、データの全体像を理解しやすくなる。 個体間・クラスタ間の類似尺度の設計が鍵 データの形式が離散値（ゲノムの配列）であったり、グラフのような構造化された特徴 量の場合にも、何らかの方法で類似尺度を定義できれば、K平均法や階層型クラスタ リングを行うことができる。 ※ 初学者向け 神嶌 敏弘, "データマイニング分野のクラスタリング手法(1) − クラスタリングを使ってみよう！ −", 人工知能学会誌, vol.18, no.1, pp.59-65 (2003) 中村 永友 (著) 『多次元データ解析法 (Rで学ぶデータサイエンス 2) 』 共立出版 ※ 中級者以上 小西 貞則 (著) 『多変量解析入門―線形から非線形へ』 岩波書店 C. M. ビショップ (著) 『パターン認識と機械学習 上・下』 シュプリンガー・ジャパン 79 時系列データ解析入門 80 本講義で概観する時系列解析の技法 時系列データの記述統計的な処理（データの視覚化、自己相関や相互共分散を 使った動的変動特性の要約） 自己回帰モデル (AR model: autoregressive model) : 時系列データに対する最も 標準的なモデリング技法 多次元自己回帰モデル (VAR model: vector autoregressive model): ARモデルを 多変数モデルに拡張したもの 時系列データ（Time Series Data）: 時間とともに不規則に変動する現象の記録 81 時系列データの例 ① (a) 月毎に測定した太陽黒点数の変動 (1749-1983, 2988観測値) (b) 1974-1979年英国の呼吸器系疾患 (肺癌、肺気腫、喘息) による月別 死亡者数 82 時系列データの例 ② 薬剤を投与 薬剤を投与 薬剤耐性を獲得した細胞 0-24 hours 薬剤に対して感受性のある細胞 がん細胞の17,654遺伝子の働きを24時間モニタしたもの 83 様々な時系列データ 定常時系列と非定常時系列 定常: 確率的特性が一定で時間的に変化しないもの 非定常: 確率的特性が時間とともに変化するもの 連続量と離散測定値 1変量時系列と多変量時系列データ • 株価・為替レート • 地震波の東西成分 • 小売食品業に従事する労働者人口 • 細胞内のタンパク質発現量の時間変化 • 船舶の横揺れ、縦揺れ、エンジンの回転数 • 消費者の購買履歴 84 時系列データ解析の目的 記述（Description） 時系列データを図示したり、基本記述統計量を計算してデータの特性を表現する。 モデリング (Modeling) 時系列モデルを設計・推定し、データの時系列特性を理解したり、将来のデータに対する予測式 を構築する。 予測 (Prediction) 現在までにえられたデータをもとに将来のデータを予測する。 制御 (Control) 操作可能な変数を適当に変化させ、制御変数の望ましい変動を実現する。 分析 データの測定 統計処理 モデリング 予測 制御 85 時系列データの1次処理 実際の時系列解析では多くの場合、モデリングや予測などを行う前に、計測値に対し て適切な変換を施す必要がある。 (a) 変数変換、 (b) 階差、 (c) 前年比・前期比、(d) 移動平均平滑化、など (例) あるサーバー経由でインターネットに接続したユーザの数（毎分）（観測数100) 原データには、強いトレンド成分が観測される 原データの１階差分をとることで、ある程度トレンド成分を除去できる 86 変数変換が有効な例 1960-1980年の Johnson-Johnson の一株当たり四半期毎利益(単位ドル) ベースラインの上昇に伴い分散が増大する傾向がある。 分散不均一性のあるデータの時系列解析は取り扱い が難しいのであらかじめこれを除去しておきたい。 対数変換を施すことで、分散不均一性を緩和 対数変換を施しても、トレンドは除去できない トレンド: 時間によって連続的に変化するベースラインの変化 87 差分によるトレンドの除去 データのトレンドはほぼ1次関数で表現可能 yt t 1, , n 元のデータは時間 t の1次式で表される。 yt a bt t : 時間; a: ベースライン; b: 勾配 元のデータは時間 t の2次式で表されるなら yt a bt ct 2 xt yt yt 1 １階差分をとると、トレンドを除去できる。 xt yt yt 1 b t に依存しない 2階差分をとると、トレンドを除去できる。 xt xt 1 yt 2 yt 1 yt 2 2c 88 移動平均フィルタ 現時刻を中心に前後（2K+1個）のデータの平均を計算し（移動平均）、時系列データのノイ ズを除去する（フィルタ）。 xt 1 y 2K 1 tK yt yt K 原データ(赤)に対して、ラグ数の K=1からK=12の移動平均フィルタ を適用した結果 ラグ数Kを大きくとるにつれ滑らかな補 正値がえられる。 89 平均、分散、自己共分散/自己相関関数 時系列データはある確率過程にしたがう確率変数の n 個の実現値と見なす 確率過程 Yt t 1, 2, , 時系列データ yt 確率過程を特徴付ける量 平均 t E Yt 分散 2 St Var Yt E Yt t 自己共分散関数 自己相関関数 C t ,k Cov Yt , Yt k E Yt t Yt k t k Rt ,k Corr Yt , Yt k E Yt t Yt k t k Var Yt Var Yt k 一般には、時間 t やラグ数 k の関数であることに注意 90 (弱) 定常性の定義 確率変数の平均と分散の値が時点 t に依存せず, 自己相関係数（自己共分散 関数）の値が 2時点の差 k のみに依存するとき、確率過程は定常であるという。 平均 E Yt 分散 2 S Var Yt E Yt t 時刻 t に対して不変 自己共分散関数 自己相関関数 C k Cov Yt , Yt k E Yt Yt k Rk Corr Yt , Yt k E Yt Yt k Var Yt Var Yt k ラグ数 k の関数であることに注意 91 自己共分散/自己相関関数の推定 1変量時系列データに関する最も基本的な記述統計量は 自己共分散(auto-covariance)と 自己相関係数(auto-correlation) である。定常時系列があたえられたとき、平均、自己共分 散および自己相関関数の推定値は次式により計算される。 標本平均 標本自己共分散関数 標本自己相関関数 1 n ˆ yt n t 1 1 n ˆ C k yt ˆ yt k ˆ n t k 1 Cˆ k ˆ Rk Cˆ 0 確率過程が定常の場合、サンプル数が増加するにつれ、真の平均、 自己共分散/自己相関関数に収束することが保証される。 92 （例） 自己相関関数の推定 米国四半期毎の大統領支持率 (ギャラップ社調査) 女性の10分間隔の血液中の黄 体形成ホルモン量の時系列 正の自己相関が 次第に減衰 周期的減衰？ 93 相互共分散/相互相関関数（多変量時系列の基本統計量） 多変量時系列データの場合には、平均、自己共分散および自己相関関数以外に、相互共分 散や相互相関関数と呼ばれる変数間の時間依存性を表す基本統計量が重要である。 ( 時刻 t の変数 i ) yi , t ( 時刻 t - k の変数 j ) y j ,t k 変数 i と j の相互共分散関数（ラグ k ） Ck i , j Cov Yi ,t , Y j ,t k E Yi ,t i Y j ,t k j スケーリング 変数 i と j の相互相関関数（ラグ k ） E Yi ,t i Y j ,t k j Rk i , j Corr Yi ,t , Y j ,t k Var Yi ,t Var Y j ,t k 94 相互相関関数の推定 標本相互相関関数 標本相互共分散関数 1 n ˆ Ck i , j yi ,t ˆ i y j ,t k ˆ j n t k 1 Rˆ k i , j Cˆ k i , j Cˆ 0 i , i Cˆ 0 j , j 1974-1979年英国の呼吸器系疾患(肺癌、肺気腫、喘息)による月別 死亡者数 男性 女性 標本相互共分散 95 時系列データの予測 これまでは、時系列データの統計処理や基本統計量にもとづく動的特性の抽出方法を 概説してきた。次に、時系列解析の実運用において最も重要な「予測」の概念を説明す る。 例えば、天気予報における気温予測、主要先進国の次四半期GDP ..... M 次線形予測子 時刻 t -1 にえられた観測値集合 y s s 1, 2, , t 1 時刻 t の値 yt を予測することを考える 簡単化のため、時系列の平均は時刻によらず０としておく 直近のM個のデータの線形結合を使った予測方式（線形予測子） M yˆ t am yt m （M 個の係数） m 1 a m m 1, , M 予測誤差 t yt yˆ t 96 最良線形予測 最良の予測方式をどのように構成するか（係数の決定） 直近のM個のデータ y t m m 1, 2, , M t 1, 2, これらが無相関かつ平均0で分散一定になるように結 合係数をあたえる。これを「最良線形予測」という。 予測誤差の列 t yt yˆ t t 1, 2, 直観的理解―過去M 時点のデータが持つ情報を余すことなく利用して yt を予測する ということは、「誤差系列」の予測をこれ以上改善することができないことを意味する。 「過去のデータと誤差系列の間」に無相関性が成り立つように予測子を構成すると、「予 測誤差の系列は互いに相関しない平均０、分散一定の確率変数（白色雑音）になる。 97 自己回帰モデル（Autoregressive Models） 観測値と線形予測子を結ぶ関係を表現したものを考える。 M yt a m y t m t m 1 E t 0, Var t 2 例えば、誤差項に正規分布を仮定 パラメータ： 2 , a1 , , aM (分散は時刻によらず一定) t ~ N 0, 2 定常性の条件： 上の式から生成される時系列が定常確率過程であるためには、M 個 の係数よってあたえられる z に関する特性方程式の根が全て１より大きくなることが必要。 1 a1z a2 z 2 aM z M 0 誤差項を0においたとき、どんな初期値から出発しても 時刻を無限にすれば yt は0に収束する。 特性根に関する条件が満たされるとき、M 次の自己回帰モデルと呼ぶ。 98 自己回帰モデルの推定 最尤推定は以下のARモデルの対数尤度を最大化すればよい。 log p y1 , , yn log p y1 , , yM n t M 1 log p yt yt 1 , , yt M しかしながら、厳密な最尤推定値を求めるためには、数値的に最適化問題を解か なくてはならないので、近似解を求めるために以下のような簡易計算手法が提案さ れている。 • • • • 最小二乗法 Yule-Walkerアルゴリズム Levinson-Durbinの逐次公式 バーグアルゴリズム 詳しく勉強したい人向けの参考図書 • 北川 『時系列解析入門』 岩波書店 • 赤池(監修) 尾崎（編） 北川（編） 『時系列解析の方法』 朝倉書店 99 自己回帰モデルの適用例 AR次数は M = 1～10 yt M a m 1 m yt m t t ~ N 0, 2 月毎に測定した太陽黒点数の変動 (1749-1983, 2988観測値) 実測値 フィッティング 予測（25時点） 次数が小さいと、予測値系 列に周期が出ない 高い次数のモデルを使うことで、 時系列データの周期特性を捉 えることに成功 100 AICによる次数選択 自己回帰モデルにもとづく予測では次数選択がとても重要 赤池情報量規準 ( AIC: Akaike Information Criterion ) を最小にする次数を選択 AIC (M) = － 2 ( M 次ARモデルの最大対数尤度 ) + 2 ( パラメータ数 ) = － モデルの適合度 + 高次モデルを使うことに対する罰則 次数を1から20に変化させたときのAICの推移 M = 9 を選択 101 多変量ARモデル（ VAR: Vector Autoregressive Models ） 互いに時間的依存性のある多変数の時系列データが得られる場合、予測する変数 以外の情報を用いて予測を行うことが望ましい。 yt y1,t , , y p ,t T VAR モデル M yt Am yt m t m 1 E t 0, , 0 T E t sT O p p t s Am は p×p の係数行列で、要素は p 個の変数間のラグ m の時間依存関係を表す VARモデルの確率特性や推定方法については、参考図書を参照 ※１変量ARモデルと同様に、最小二乗法やYule-Walkerアルゴリズムを使って簡単に推定できる。 • 北川 『時系列解析入門』 岩波書店 • 赤池(監修) 尾崎（編） 北川（編） 『時系列解析の方法』 朝倉書店 102 多変量ARモデルの適用例 (1/2) 米国の経済指標（4変数）の推移: (1) マネタリベース (M1) ; (2) 実質GNP; (3) 91日物TB金利; (4) 長期債金利 原系列 差分によるデータ処理 103 多変量ARモデルの適用例 (2/2) AICにもとづくモデル選択 M = 1～20 M=2 次数２で最小になるが、次数の変化に 対してAICの値のばらつきが大きく、信 頼性は低い。とくに、次数2と次数5の差 はとても小さい。 M=5 104 時系列解析のまとめと補足 本講義では、(a) 時系列データの統計処理、(b) 基本統計量にもとづく時系列特性の 要約、(c) 自己回帰モデルにもとづく予測について概観してきた。 時系列解析のほんの一部に触れたに過ぎない。時系列解析を本格的に学ぼうという 意思のある方は、本講義では触れることの出来なかった次の内容についてもフォロー して下さい。 • 季節調整法 • トレンド解析 • 周波数領域の分析手法： スペクトル解析、ピリオドグラム • ARMA ( autoregressive moving average ) モデル • 状態空間モデルによる非定常時系列データの解析 • モデル選択 参考図書 • 北川源四郎 『時系列解析入門』 岩波書店 • 赤池 (監修) 尾崎（編） 北川（編） 『時系列解析の方法』 朝倉書店 105 本講義で使ったスライドのダウンロード http://daweb.ism.ac.jp/~yoshidar/index_j.htm 106
© Copyright 2024 ExpyDoc
