① 多変量解析、② 時系列解析

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

Download Report