1-4 データ分析

東京⼤学 数理・情報教育研究センター
2021年5⽉6⽇

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学
数理・情報教育研究センター
北川源四郎,
寒野善博
CC BY-NC-SA
東京⼤学
数理・情報教育研究センター
北川源四郎
20212021
CC BY-NC-SA

概要•

回帰分析，時系列分析，クラスター分析などの典型的なデータ分析
⼿法を学び，予測，分類，パターン発⾒などへの応⽤例を学びます．

•

⾼次元のデータをより扱いやすいものとする次元削減の⼿法を学びます．また，さまざまなデータ分析の⼿法の基礎となっている最適化の基本的な概念や具体例を学びます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学
数理・情報教育研究センター
北川源四郎,
寒野善博
CC BY-NC-SA
東京⼤学
数理・情報教育研究センター
北川源四郎
20212021
CC BY-NC-SA

2

本教材の⽬次
１．回帰分析

４−21

２．ロジスティック分析

22−32

３．時系列分析

33ー46

４．クラスター分析

47ー55

５．パターン発⾒56ー64

６．次元削減65ー91

７．最適化概論

92ー111

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学
数理・情報教育研究センター
北川源四郎,
寒野善博
CC BY-NC-SA
東京⼤学
数理・情報教育研究センター
北川源四郎
20212021
CC BY-NC-SA

3

１．回帰分析
いくつかの量的変数があるときに，特定の⽬的変数を他の変数によって説明するのが回帰分析で，その結果は予測に⽤いることもできます．特に，ひとつの説明変数で⽬的変数を説明する⽅法は単回帰分析，複数の説明変数で⽬的変数を説明するのが重回帰分析と呼ばれます．説明変数
を⼊⼒，⽬的変数を出⼒と呼ぶこともあります．回帰分析では，予測誤差とその分散，回帰係数の推定法や，次数や変数の選択が重要です．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

44

説明変数との関係をみる
下の棒グラフは左から北海道，⻘森，・・・，沖縄の順に４７都道府県の⼈⼝千⼈あたりの新型コロナウイルスの累積感染者数を⽰しています．都道府県によって感染した⼈の割合が⼤きく異なることがわかります．
東京
感染者数︵千⼈あたり︶

沖縄

⼤阪
福岡

北海道

47都道府県

感染者数︵千⼈あたり︶

このような都道府県（以下では単に県と書きます）による感染者の違いを各県の県庁所在地の⼈⼝密度と関連づけて説明することを考えてみましょう．
実際に，各県の県庁所在地の⼈⼝密度を横軸に，千⼈当たりの感染者数を縦軸にとって散布図を描いてみると，右のように⼈⼝密度と感染者数には明確な関係がみえます．

県庁所在地の⼈⼝密度
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA5

単回帰分析

右の散布図において，注⽬する変数を縦軸に表し，⽬的変数あるいは出⼒と呼びます．⼀⽅，横軸には説明変数あるいは⼊⼒を表します．

感染者数

そこで，⼈⼝密度の情報を利⽤することによって，感染者数の値を説明したり，予測したりすることが考えられます．

【注意】前⾴の散布図と違って，右図では⼈⼝密度の対数を取って表⽰しています．このような変換によって⼆つの変数の関係が直線に近くなり，関係をとらえやすくなることがあることに注意すべきです．

県庁所在地の⼈⼝密度の対数

単回帰分析では，⽬的変数（出⼒） y の値を

y  a  bxのようにひとつの説明変数（⼊⼒） x を⽤いて回帰直線で表します．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

6

回帰モデル
ただし，実際のデータではすべての点がこの直線上にのるわけではありません．そこで回帰直線と説明変数から定まる予測値と実際の値との差 n を予測誤差とみなし，回帰モデル

yn  a  bxn   nを考えることにします．a と b は回帰係数と呼ばれますが，特に a を切⽚,
b を直線の傾きと呼ぶこともあります．

y
以下では回帰係数をデータから推定する⽅法を考えます．

実測値

yn
予測誤差

n

回帰直線

a

予測値

xn
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

説明変数

x
7

最⼩⼆乗法による回帰直線の推定
回帰モデルの係数 a と b の推定にあたっては，なるべく予測誤差 n を⼩さくすることを⽬指します．そのためにN 組のデータ ( x1 , y1 ), ,( xN , y N ) があると
き予測誤差 n の⼆乗の総和

S (a, b) NN    y  a  bx 
2
n 

n 1

n

2

n

n 1

を最⼩とするように係数 a と b を求めます．このように誤差の⼆乗和を最⼩化することによって未知数を推定する⽅法を最⼩⼆乗法と呼びます．
最⼩⼆乗法を適⽤すると係数の最⼩⼆乗推定値が次のように得られます．デー
タから推定された係数などには「＾」をつけて表現します．

 (x  x )( y  y ) (x  x )Naˆ  y  bx ,

bˆ 

n

n1

n

N

n1

2

n

なお，以下のように S ( aˆ , bˆ) をデータ数 N で割ると残差分散が得られます．

ˆ 2  S ( aˆ , bˆ) / N
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

8

【参考】最⼩⼆乗法の解を求める
最⼩⼆乗法では評価関数 S(a,b) を最⼩にしますが，そのために S(a,b) を偏微分したものを０とする a と b を求めます．実際にS ( a, b) 

N



2
n 

n 1

N



 yn  a  bxn 

2

n 1

を a と b について偏微分を計算してみると
NS (a, b) 2  yn  a  bxn   0,
a
n 1



N
S (a, b) 2 xn  yn  a  bxn   0
b
n 1



となるので，この連⽴⽅程式を解いて，前⾴の解が得られます．

回帰係数などの未知数の推定にあたって，残差の⼆乗和を最⼩にすることは，残差 nが正規分布に従っていると仮定しているものと考えることができます．この想定する分布の形が異なると，推定結果も異なってきます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

9

例：新型コロナウイルス感染者数データ
東京a  7.7792,

b  3.1890となります．したがって，回帰直線はy  7.7792  3.1890 x

感染者数︵千⼈あたり︶

右図のCovid-19データの場合に a と b を求めてみるととなって，⾚線で⽰した直線が得られます．
県庁所在地の⼈⼝密度の対数

最も右上の点(東京)は誤差が⼤きいですが，それ以外の県についてはよく傾向を表現できているといえます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

10多項式回帰モデル
単回帰モデルはデータに直線をあてはめますが，⼊⼒と出⼒の関係が直線では⼗分に表現できない場合には， 多項式回帰モデル

yn  a0  a1 xn    am xnm   n
を使うとよいことがあります．ここで，m は多項式の次数，aj は回帰係数と呼ばれます．m=1のときは直線，m=2の時には2次式となります．

yn  11.499  9.026 xn  1.891xn2
です．⾚の直線よりも，データに良く適合しているように⾒えます．

感染者数︵千⼈あたり︶

右図の⻘⾊の曲線は m=2 として得られた2次の多項式

県庁所在地の⼈⼝密度の対数
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

11表現⼒が⾼いモデルを使う時の注意：

実際，8次の多項式を使ってみると右上の図の緑⾊の曲線のように，データの両端ではちょうどデータの上を通っていてデータを忠実に表現しています．

感染者数︵千⼈あたり︶直線を⽤いた単回帰の代わりに2次式を⽤いると良い当てはまりが得られたので，もっと⾼い次
数を使ってみたくなるのは⾃然です．

県庁所在地の⼈⼝密度の対数
3

実際，残差分散をプロットしてみると，次数が⾼くなると残差分散は単調に減少しており，⾼次のモデルほどデータへの適合がよいことが分ります．しかし，この8次の回帰曲線は本当に良いモデルでしょうか？

残差分散

2

最⼩
1

0
0

1

2

3

4

5

6

7

8

次数

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

12表現⼒が⾼いモデル利⽤時の注意：

過学習

表現⼒が⾼すぎるモデルを使うと，データにはよく適合するが観測されたデータ以外の予測，特に外挿を⾏ったときの予測精度が著しく悪くなる恐れがあります.
実際，説明変数の範囲を2.0から4.5に広げて回帰曲線を描いてみると図のようになります．外挿区間を広げても直線や2次式は安定していますが，緑⾊の8次式は両端で跳ね上がっていて，外挿能⼒が低いことがわかります．

感染者数︵千⼈あたり︶

データをよく表現しているが，予測能⼒が低くなる現象は過学習あるいは過剰適合(オーバーフィッティング)と呼ばれます．学習に使ったデータだけはよく表現しているが，それ以外のデータにも使える汎⽤性がないのです．

データがある区間

県庁所在地の⼈⼝密度の対数区間を広げた東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA13

次数選択，モデル選択
多項式の次数は低すぎると現象を⼗分に表現できませんが，⾼すぎると過剰適合を起こします．したがって次数を適切に選ぶことが重要で，次数選択の問題と呼ばれます．
予測の観点からモデルの良さを評価し，次数を⾃動的に選択するための評価基準として情報量規準 AICやBICが⽤いられます．
候補となるすべての次数 m に対してAICm (あるいは BICm ）を計算してその中で最⼩にするmを探すと，それが最適な次数となります．
AIC
(BIC)

AIC(または)BIC

AICは⾚池情報量規準，BICはベイズ型
情報量規準と呼ばれます．

次数
最適な次数
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

14情報量規準AIC と BIC
m 次の多項式回帰モデルの場合， N をデータ数， m はm次の多項式回帰モデルの残差分散とするとき，AIC と BIC は（共通の定数を無視すると）次のように定義されます．
2

AIC m  N log 2 m2  2( m  2)
BIC m  N log 2 m2  ( m  2) log N
データからの乖離

モデルの複雑さへのペナルティ

ここで，右辺の第１項はデータからの乖離を，第２項はパラメータ数を増やすことによるモデルの複雑さへのペナルティ項と考えることができます．
AIC
(BIC)

データからの乖離は次数を増やすと減少する⼀⽅，ペナル
ティ項は増⼤するので，通常適当な次数で最⼩となります.

データから
の乖離

【注意】上記の定義は回帰モデルの場合で，情報量規準はもっと⼀般のモデルに適⽤できます．また，AICやBICはモデルの次数選択だけでなく，モデルの選択や変数選択にも利⽤できます.

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

AIC(または)BIC

ペナルティ

次数

15

次数選択の例
表は新型コロナウイルス感染者数（千⼈当たり）データに様々な次数の多項式をあてはめた場合の残差分散, AIC, BICを計算した結果です．
次数が増えると，データへの適合度に対応する残差分散は単調に減少しますが,AICやBICは次数2で最⼩となっており，このデータの場合は 2次曲線が適当であることがわかります．
13ページの外挿結果を⾒ると，AIC や BIC が残差分散はある程度⼤きくても，低い次数の多項式を選択しているのは妥当な判断であることがわかります．

東京⼤学 数理・情報教育研究センター

次数

残差分散

AIC

BIC

0

2.776

185.37185.221

0.650

119.14120.84

2

0.403

98.62

102.173

0.401

100.39105.794

0.391

101.27108.525

0.366

100.08109.18

6

0.349

99.88

110.83

7

0.329

99.10

111.908

0.325

100.49115.144-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

16

例：実効再⽣産数について（簡易）実効再⽣産数
1.8

上の図は東京におけるコロナ感染者の(簡易)実効
再⽣産数の変化を⽰します．⼀⽅，真ん中の図は
都⼼の⼈出の感染拡⼤前⽐です．

1.4
1.0
0.6
9/1

9/29 10/27 11/24 12/22 1/19

【注意】実効再⽣産数や都⼼⼈出は独⽴な観測値ではなく，時間的な相関がある時系列なので，厳密には後に学ぶ時系列
データの分析法を適⽤すべきです．

4/13

3/16

4/13

-20
-30
-40
-50
9/1

9/29 10/27 11/24 12/22 1/19

2/16

1.4

実効再⽣産数

そこで，回帰直線を引いてみると，⼈出が感染拡
⼤前と⽐較して30%以上減少しているときには実効再⽣産数が 1 以下となって，新規感染者数が減少していることが⽰されます．

3/16

都⼼⼈出（感染拡⼤前⽐）-10下の図は都⼼の⼈出率を横軸に，実効再⽣産数を縦軸とした散布図です．⾚の点線で⽰すように⼈出が減少すると実効再⽣産数が⼩さな値をとる傾向が⾒えます．

2/16

1.2

1

0.8

0.6
-40

-35

-30

-25

-20

⼈出率

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

17

17

重回帰分析
単回帰分析の例では，各県の千⼈あたりの感染者数を県庁所在地の⼈⼝密度で単回帰しましたが，感染者数に関連すると思われる変数は県全体の⼈⼝密度
⾼齢者率
交通便利度など，その他にもいろいろ考えられます．したがって，これらの変数を説明変数と考えて，⽬的変数を
1番⽬の
説明変数

m番⽬の
説明変数

残差

yn  a0  a1 xn1    am xnm   n
m 個の説明変数のように，定数項 a0 とm 個の説明変数 x1,…,xmの影響として表現するモデル

が重回帰モデルです．単回帰の時と同様に a1 ,,…,am は回帰係数で，残差 n は正規分布 N(0, 2)に従うものと仮定します．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

18

重回帰分析：都道府県別感染者数の例
以下ではコロナウィルスの都道府県別感染者数を例に考えてみます．説明変数の候補として，表のように xn1 から xn8 までの 8変数を考えることにします．xn1は定数項に対応します．xn8 は県庁所在地の⼈⼝密度 xn4 の2乗です．この変数は多項式回帰で2次のモデルが最もよかったことから，採⽤したものです．

これらの説明変数は，コロナウイルスの感染が⼈⼝密度や年齢あるいは⼈の移動に関連して発⽣しやすいと⾔われていることから，候補として採⽤したものですが，実際にはこのほかにもいろい
ろな変数を考えることができます．

東京⼤学 数理・情報教育研究センター

表：⽬的変数と説明変数
yn

コロナウイルスの累積患者数

xn1

定数項

xn 2

県全体の⼈⼝密度

xn 3

⾼齢化率

xn 4

都道府県庁所在地の⼈⼝密度

xn 5

交通便利度(2時間）

xn 6

（4時間）
交通便利度

xn 7

交通便利度
（6時間）

xn8

xn2（⼈⼝密度の2乗）
44-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

19

重回帰分析：説明変数の選択
m 個の変数について，モデルの変数として利⽤する場合と利⽤しない場合を考えると，全部で 2m 種類のモデルがあります．現在の例では
m=8 なので，256通りのモデルがありますが，部分回帰という⽅法を⽤いると，説明変数の数
1,…,8 のそれぞれについて，すべての候補の中で誤差分散が最も⼩さなモデルを⾃動的に求めることができます．
AICで判断すると，表から x1 (定数項）, x3 (⾼齢化率), x4 (県庁所在地の⼈⼝密度), x8 (⼈⼝密度の2乗) の 4変数を⽤いたモデル

変数
の数

残差分散

AIC

1

1.241

147.53

x8

2

0.460

102.89

x3,x8

3

0.403

98.62

x1,x4,x8

4

0.312

88.70

x1,x3,x4,x8

5

0.309

90.26

x1,x3,x4,x5,
x8

6

0.303

91.22

x1,x3,x4,x5,
x6,x8

7

0.287

90.63

x1,x3,x4,x5,
x6,x7,x8

8

0.272

90.21

x1,x2,x3,x4,
x5,x6,x7,x8

説明変数

yn  16.3230  0.1611xn 3  8.5340 xn 4  1.7072xn8   nが最もよいことがわかります．
20
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

説明変数の選択（続）前⾴のAICの値をみると，2次の多項式モデルに⾼齢化率を加えたモデルが最もよいことになります．
また，変数を1個だけ使うとすると⼈⼝密度の2乗（xn8）が最もよく，2個を使うとすると⼈⼝密度の2乗と⾼齢化率（xn3）がよいこともわかります．
ちなみに，このときのモデルはyn  0.0849 xn 3  0.4573xn8   nとなり，⾼齢化率の回帰係数がマイナスですが，
これは⾼齢化率が上がると患者数が減少するという因果関係を意味しているのではなく，図のように⾼齢化率が⾼い県では患者数が少ないという相関関係を表しているにすぎないことに注意しましょう．
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

21

２．ロジスティック回帰分析

ロジスティック回帰分析は{0, 1}のようなカテゴリー(離散)変数やその観測値を集約した割合などを説明変数で表現する⽅法です．形式的に普通の回帰モ
デルをあてはめると不都合なことがおこるので，0, 1が発⽣する確率を考え,その⾮線形変換したものの回帰モデルを考えます．
本稿では最初に，観測値の⾮線形変換に対して普通の回帰モデルを適⽤する⽅法を考え，次に本来のロジスティック回帰モデルを考えます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

2222割合やカテゴリ変数のデータ
ロジスティック回帰分析は {0, 1} のような離散変数やそれらを集約した割合などを説明変数で表現する⽅法です．
例として，携帯電話の普及率の推移と湿度と降⾬の関係を考えてみましょう．
左図：1993年から2019年までの携帯電話（スマホを含む）の世帯普及率(%)の推移を⽰します．20世紀に急激に⽴ち上がってその後は100%近くに張り付いています．
⻄暦を説明変数，普及率を⽬的変数とします．
右図：東京の1年間の気象データ．横軸に1⽇の最低湿度(%)，縦軸は降⾬で⾬が降った
場合は y=1, 降らなかった場合は y=0 の2値データ．湿度を説明変数，降⽔を⽬的変数とします．
携帯電話普及率

降⾬事象

100
80
60
40
20
0
1993 1996 1999 2002 2005 2008 2011 2014 2017

最低湿度（%）
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

2323

回帰分析における問題点
割合データや予測誤差を定義できない離散変数の場合でもその事象が発⽣する確率を想定することによって，形式的には通常の回帰モデルを適⽤することは可能ですが，いくつかの問題があることを理解しておく必要があります.
1.

単純な単回帰モデルではその確率が０以下や1以上になるという不都合なことが起こります（左図）．

2.

割合データでは，特に発⽣確率が０と１の近辺で誤差分布が正規分布から⼤きく
外れます（右図のように場所によって分布の形が変化します）．

3.

変動の分散が⼤きく変化します．

6
降⽔事象︵●︶と降⽔確率︵⾚︶

誤差分布
4
2
0
0.0
最⼩相対湿度(%)

東京⼤学 数理・情報教育研究センター

0.2

0.4

0.6

0.8

1.0

発⽣確率

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

24

24

ロジット変換
まず，携帯電話普及率のように，説明変数（年）に対する割合 pn が⼤量な観測数に基づいて得られている場合には，pnをそのまま回帰分析するのではなく， p のロジット変換
ロジット変換の逆変換
1

 p 
q  log 

1

p



を定義しそのモデル化を考えます．逆に q から
qep1  eq

0.8
0.6

p 0.4
0.2
0
-6

-4

-2

0

q

2

4

6

で確率 p を求めると，q がどんな値をとっても p は常に 0 と 1 の間の値をとる
ので q の空間では⾃由なモデリングが可能となります．
【参考】p/(1p) はオッズと呼ばれ，競⾺などのギャンブルでも利⽤されます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

25

25

ロジット変換（例）
携帯電話普及率（上図）の場合には，ロジット
変換によって，下図が得られます．このロジット変換 q は−∞から∞までの値をとれるので，説明変数 xn を⽤いてqn  a0  a1 xn    am xnm

0.8
0.6
0.4
0.2
0.0
1993 1996 1999 2002 2005 2008 2011 2014 2017と回帰できます．このとき，q の逆変換によって y が 1 をとる確率がexp( a0  a1 xn    am xnm )pn 
1  exp( a0  a1 xn    am xnm )で求められます．

携帯電話普及率

1.0

携帯電話普及率のロジット変換
2
1
0
-1
-2
1993 1996 1999 2002 2005 2008 2011 2014 2017東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

26

26

ロジット変換（例2）

0.2

致死率
0.15

コロナウイルスの致死率（感染者のうち亡くなった⼈の割合）を考えます．致死率は年齢に⼤きく依存することが知られているので，横軸に20代以下，30代などの世代をとり，縦軸に致死率をプロットしてみると上図のようになります．
世代が増えるに従って致死率が上昇する傾向があきらかですが，このデータのロジット変換は中図のようになり，直線でかなりよく近似できます．
各世代を真ん中の年齢で代表させて単回帰を⾏うと

qn  13.108  0.129 xnという直線(⻩⾊）が得られます．したがって，

exp(13.108  0.129 xn )pn 
1  exp( 13.108  0.129 xn )で確率 pn の予測式（⻩⾊）が得られます．90歳以上ではデータよりもかなり⼤きな値が予測値になっています.
東京⼤学 数理・情報教育研究センター

0.1
0.05
0

0

ロジット変換
ロジット変換

0
‐2
‐2
‐4
‐4
‐6
‐6
‐8
‐8
‐10
‐10
‐12
‐12

ロジット変換
ロジット変換
予測値
予測値

0.3

致死率
0.2
致死率

0.1

予測値

0

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

27

27

ロジスティック回帰モデル
N 組のデータ ( x1 , y1 ), ,( xN , y N ) があり，xn は説明変数，yn は病気に罹患している場合は0，いない場合は１となる離散変数とします．
0 か1の2値をとる変数に対して予測誤差は考えにくいので，1という事象がおきる確率 p を考え，そのロジット変換 q を説明変数 xn で表現するモデル
qn  a  bxn

を考えることにします．このとき
 pn log   a  bxn
 1  pn 

なので，事象１および０が起きる確率はそれぞれ次のように表されます．
pn 

exp(a  bxn ),1  exp( a  bxn )1  pn 1
1  exp( a  bxn )このように事象１と０が起きる確率を表現するモデルをロジスティック回帰モ
デルと呼びます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

28

28

ロジスティック回帰モデルの尤度
ロジスティック回帰モデルを想定すると，説明変数 xn が与えられているとき，
yn の確率は11  exp( a  bx )n
P ( yn | xn , a, b)   exp( a  bxn )1  exp( a  bxn )exp{ yn ( a  bxn )}1  exp( a  bxn )yn  0 のとき
yn  1 のとき【注意】 yn は 0 か１なので上の式はこのように⼀つの式で書ける．

で与えられます．したがって，N 組の独⽴なデータ ( x1 , y1 ), ,( xN , y N ) が与えられているとき尤度は以下のように定義されます．
NN

exp{ yn ( a  bxn )}L ( a, b) p ( yn | xn , a, b) n 1n 1 1  exp  a  bxn 

次⾴ではこの尤度を⽤いてパラメータ a, b を推定する⽅法を考えます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

29

29

最尤法
尤度の最⼤化によって，未知数   (a, b) を推定する⽅法は最尤法と呼ばれます.
ただし，通常は尤度の対数をとった以下の対数尤度を最⼤化します．対数関数は単調増加関数なので，どちらでも同じ推定値が得られます．
ロジスティック回帰モデルの場合，対数尤度は次のようになります．
( a, b)  log L( a, b)N  y (a  bx )  log 1  exp (a  bx )
n

n

n

n 1

したがって，この対数尤度を最⼤化，すなわちmax ( a, b)  ( aˆ , bˆ)
a ,bとなる aˆ と bˆ を求めることによって，係数の推定値が得られます．このように，最尤法で推定された係数を最尤推定値と呼びます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA3030

例：降⽔確率降⽔

東京の1年間（2020年）の気象データを考えます．横軸に1⽇の最低湿度(%)，縦軸は
降⾬で⾬が降った場合は y=1, 降らなかった
場合は y=0 の2値データ．湿度を説明変数,
降⽔事象を⽬的変数として，湿度からどの程度降⾬を予測できるかを⾒てみます．

最低湿度（%）

ロジスティック回帰モデルを推定してみるとaˆ  5.2652bˆ  0.0892という推定値が得られます．したがって，最低湿度が x の時の降⽔確率はp( x) exp( 5.2652  0.0892 x )1  exp(5.2652  0.0892 x)となります．
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

31

31

例：降⾬確率

降⽔確率
最低湿度

⾚い曲線はロジスティック回帰モデルで最低湿度に基づき推定した東京の降⽔確率です．湿度が40%，60％，80％の時，降⽔確率がそれぞれおよそ
15%，50％，85%となることがわかります．この⾚い曲線が縦に⽴っているほど予測精度が⾼いことになります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

32

32

３．時系列分析
時間と共に変動する現象の記録が時系列です．ここでは，時系列データの簡単な分析の仕⽅を学びます．より詳しい説明がリテラシーレベル教材の4-4にあります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

33

33

時系列とは3000
新規感染者数︵東京︶

図は横軸に⽉⽇，縦軸に東京の新型コ
ロナウイルス新規感染者数を⽰します．
このように時間と共に変化する現象を記録したものが時系列データで，それを図⽰したものは 時系列グラフと呼ば
れます．

2500
2000
1500
1000
500
0
6/1 6/29 7/27 8/24 9/21 10/1911/1612/14 1/11 2/8

3/8

4/5

データが何番⽬に観測されたかを n で表し，そのときの観測値を ynとします．
N 個の観測値 y1,…,yN が得られたとき，N は時系列の⻑さあるいはデータ数と呼ばれます．
時系列分析の⽬的は，時系列の特徴を捉えたり，モデルを推定することによって，今後の変化を予測したり，意思決定に必要な情報を得ることです．そのための基本的なツールとして，以下では時系列の周期を捉える⽅法と時系列の傾向を捉える⽅法を考えます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

34

34

時系列を変換してみる
4.0
新規感染者数の対数

⼈数や⾦額などを数えたデータの場合は対数をとった⽅が特徴を捉えやすいこともあります．右図は新規感染者数の常⽤
対数をとって表⽰したものです．

3.0

2.0

1.0
1

29

57

85

113

141

169

197

225

253

281

309

前ページの原データでは，時刻によって変動の幅が⼤きく変化していましたが，対数変換の結果ほぼ同じような変動幅となり，同じようなパターンで周期的に上下を繰り返していることが顕著になります．この周期性は，検査体制の影響で⽇曜⽇や⽉曜⽇の新規感染者数が少なくなるからです．
ただし，時系列の場合は完全な周期的変動を⽰すことは少なく，多くの場合周期的な変動パターン⾃体が徐々に変化します．このような時系列の周期的変動を可視化するツールとして，ピリオドグラムとパワースペクトルがあります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

35

＊リテラシーレベル教材4.4から流⽤

パワースペクトル

光をプリズムで沢⼭の⾊に分解できるように，時系列はサイン・コサインの和に分解できます．その強さを表⽰するとどの周期が強いかが分ります.

プリズム

同じようにランダムな波形も周期関数の和で表現できます。
 スペクトルアナライザや⾳声のグラフィック・イコライザが例です

スペクトル
アナライザ

2020
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

36

＊リテラシーレベル教材4.4から流⽤

周期関数の場合• サイン・コサインのようなひとつの周期関数で表現できる
場合にはスペクトルは1点に集中します．
p( f )

y n  sin ( 2  f 1 n )

n

f

f1

• ２つの周期関数の和の場合はスペクトルは２つになります．
p( f )

y n  sin(2  f1 n )  0.8 sin(4  f 1 n )

n

f1

f2

f

( f 2  2 f1 )

縦棒の⾼さは振幅の2乗に⽐例します．
2020
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

37

例：

Covid-19データの周期性
1週間の周期

Covid-19 時系列のピリオドグラムの例です．横軸 f は周期の逆数で周波数と呼ばれます．また，縦軸は各周波数におけるピリオドグラムの対数を表します． 0  f  0.5 の範囲が図⽰されていますが，f = 0.5
は 2⽇に 1回の周期成分，f = 0 は周期無限⼤の成分です．ピリオドグ
ラムが極⼤値をとる f=1/7 = 0.1429は1週間の周期成分に対応します．
ピリオドグラムからこのデータには顕著な1週間周期が存在することがわかります．また，f=2/7 と 3/7 にもピークが⾒られますが，これは⾼調波と呼ばれるもので1週間の周期成分が単純なサイン・コサイン
ではなく，複雑な波形であることを⽰しています．

Log-periodogram

周波数 f

ピリオドグラムに似たものとしてパワースペクトルがあります．パワースペクトルはピリオドグラムと違って連続関数になりますが，時系列モデルを推定することによって求められます．パワースペクトルでも同様な周波数にピークが⾒られます．
周波数 f

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA38周期性の除去 周期性を除去する理由
Covid-19の場合について考えてみると，検査体制の事情から毎週⽇曜⽇や⽉曜⽇は⼩さな値になりますが，この事情を考慮しないと，患者数が減少したと誤った判断をしかねません．このような場合には，時系列の周期成分を除去してしまう⽅が簡単です．

 周期性を除去する⽅法
周期の⻑さが分っている周期成分を除去する簡単な⽅法として，移動平均があります．
⼀般に，時系列 y1 ,, yn が与えられるとき(2k+1)項の移動平均はtn 

1( ynk    yn    ynk )2k  1で定義されます．ここで，k=3 とすると項数は 2k+1=7 となり，カッコ内には⽇曜から⼟曜までの各曜⽇が1回ずつ含まれるので，曜⽇の影響を除去できます．この⽅法では直近の tn2, tn1, tn は計算できないので，上記の式で求めたものを現時点の推定値とみなして tnk と定義することもあります．ただし，この⽚側移動平均の場合は，k 時点だ
け波形が遅れることに注意する必要があります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

39

例：Covid-19データの場合
上図の⾚線は k=3 として求めた7項移動
平均です．1週間の変動パターンがほとんど除去されて滑らかになり，感染者の増減の傾向がよく捉えられているよ
うにみえます．

【参考】移動平均において k=10 とすると，前後の1週間を加え，合計3週間の平均をとることになり，さらに滑らかな移動平均が得られます(下図)．ただし，その反⾯ 急激な構造変化を捉えにくくなるという副作⽤があり, n=230 付近の急激な変化も滑らかになっている事には注意する必要があります.

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA40季節調整法：より⾼度な⽅法
擬似周期的な変動を含む時系列の分析法として季節調整法があります．この⽅法ではyn  Tn  Sn  wnと時系列をトレンド Tn，季節成分 Sn，ノイズ wnの3成分に分解します．（季節調整法は年周期の除去を主要な⽬的とするため，1週間の周期の場合でも季節成分と呼ばれています．詳細はリテラシー教材の4.4節)
周期を p=7，トレンド次数を2として，この⽅法を適⽤すると，右図のような結果が得られます．
⾮常に滑らかなトレンドが得られ，1週間周期の変動が適切に除去されています．また,
季節成分が時間の推移にともなって少しずつ変化していることもわかります．

東京⼤学 数理・情報教育研究センター

log10(covid-19)
トレンド

季節成分

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

41

⾃⼰相関関数
周期的変動に限らず時系列の変動の特徴は⾃⼰相関関数で表現することができます．ある時刻と k だけ離れた時刻との相関係数を計算しR ( k )  Cor( yn , yn  k )と表します．R(k) は yn と ynｰk の相関係数で k はラグ，R(k)は ⾃⼰相関関数と呼ばれます．
•••

R(k) は k だけ離れた地点との相関の⼤きさを表します．
yn とynkに正の相関があるとき R(k) > 0 となります．
ynとynkに負の相関があるとき R(k) < 0 となります．
時間差

k

nk
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SAn42

例：Covid-19 データの⾃⼰相関関数
Covid-19データのラグ32までの⾃⼰相関関数を⽰します．
• 左上は原データの⾃⼰相関関数で，7⽇周期で変動しながらゆっくり減衰しています．
• 左下は季節調整で得られたトレンドの場合で，⾮常にゆっくり減衰しています．また原データの変動の主要部分がトレンドであることもわかります．
• 右上は季節成分で7⽇周期でほぼ規則的に変動しています．
• 右下はノイズ成分で，ほとんど⾃⼰相関がない系列であることがわかります．

Covid-19
新規感染者数

季節成分

k

k

ノイズ

トレンド

k
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SAk43

時系列の予測
時系列分析の重要な応⽤が予測です．時刻nまでの時系列 y1,…,yn が与えられたとき，将来の値 yn+1,…,yn+k を推定する問題が予測です．とくに１時点先のyn+1 を推定する問題を1期先予測，yn+j ( j は2以上)を推定するときは⻑期予測と呼びます．
 ARモデルによる予測
代表的な予測⽅法はARモデルを⽤いる⽅法で，時刻 nまでのデータを使ってARモデルyn  a1 yn 1    am yn  m   n
観測値が推定できると，yn+1の予測値はyn 1|n  a1 yn    am yn  m 1となります．以下，逐次代⼊により2期先以上の予測値も求められます．

東京⼤学 数理・情報教育研究センター

予測誤差
過去の値 (m個)n m4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

予測値

yn|n1n 1 n44

例：Covid-19データの予測
AIC

まずn=260までのCovid-19データからARモデルを推定し,そのモデルを使って予測した場合を⽰します．
上図のAICで選択した結果15次のARモデルが得られます.
真ん中の図は，そのモデルで推定したパワースペクトルです. 7⽇とその1/2および1/3の周期が⽬⽴っています．

パワースペクトル

下図はこのモデルによる⻑期予測結果で，予測値及び標準誤差を⽰します．⿊⾊で⽰す実際の値を⽐較的よく予測していますが，誤差幅が広くなっています．
ARモデルによる予測は便利ですが，
このデータのように明らかなトレンドがある場合には，トレンドを分離する⽅がよい結果が得られます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

45季節調整モデルによる予測の例
季節調整法によるトレンド，季節成分，ノイズへの分解を利⽤した⻑期予測の結果を⽰します．
トレンドは⾮常に滑らかなので⻑期予測も安定して実施でき，予測誤差も⽐較的⼩さくな
ります．また，季節成分も安定しているので，その予測誤差も⽐較的⼩さくなります．
結果として，この⼆つを加算した時系列の予測値の予測誤差もARモデ
ルの場合よりも⼩さく，また周期パ
ターンもよく再現されています．
コロナウイルス感染者数の予測のためにはこの予測値 pnから10 pn を計算すればよく，1週間周期の変動も考慮した予測値が得られます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

46

４．クラスター分析
説明変数の情報を利⽤せず，データを似たもの同⼠に分類するのがクラスター分析です．どんなデータを似ていると考えるかは，２点間の距離やクラスター間の距離（類似度）の定義に反映されます．距離の定義によって結果がかなり異なることがあることは注意する必要があります．
⼆つの代表的なクラスター分析の⽅法として，階層的クラスター分析と⾮階層的クラスター分析があります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

47

47

クラスターとは
k 次元のベクトルで表されるデータを考えます．

xn  ( xn1 , , xnk )例えば，n 番⽬の学⽣の数学, 物理, ⽣物, ... , 英語の成績などです．学⽣が N ⼈いるとN 個のベクトル x1 , , xNが得られます．
k 次元空間の似た点の集まりをクラスターと呼びます．クラスター分析では，データから図の様に近いものを⼀つのクラスターにして，いくつかのグループに分けます．ただし，いくつのクラスターが存在するか前もってわかっているわけではないので，2点間の近さやクラスター間の近さをどのように定義するかによって結果が変わってきます．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

48

クラスター間の距離
⼆点間の距離としては通常，次に⽰すユークリッド距離を使います．

xn  ( xn1 , , xnk ),yn  ( yn1 , , ynk )d ( x, y )  ( x1  y1 ) 2    ( xk  yk ) 2

・ただし，距離の定義としては，このほかマンハッタン距離，ミンコフスキー距離，チェビ
シェフ距離などがあります．

クラスター間の距離（類似度）
クラスター間の距離の定義としては，最短距離法，最⻑距離法，重⼼法，群平均法，ウォード法などいろいろ提案されており，どれを選ぶかによってクラス
ター分析の結果が異なります．


最短距離法

東京⼤学 数理・情報教育研究センター

最⻑距離法

重⼼法

群平均法

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SAn A nBn A  nB

ウォード法

49

階層的クラスタリング
階層的クラスタリングでは、最初は全てのデータが異なるクラスターに属するものとし,（クラスターの数）=（データ数）とする状態を考えます．
次にすべてのクラスター間で最も距離が⼩さな２つを探し，それらを⼀つのクラスターに統合し，クラスター間の距離を更新します．以下，このステップを繰り返すことによってボトムアップにひとつずつクラスター数を減らしていきます．
このクラスター統合のときにΠ型にふたつのクラスターを繋ぎ，縦の⾼さをクラスター
間の距離とします．この操作を続けていくと，最終的にクラスターがひとつになると，
樹形図（デンドログラム）ができます．
既に樹形図が与えられている場合には，横線を引いて上からだんだん下げていくと2クラスタ
ー，3クラスターの順に分類ができます.
右図の場合は２クラスターの場合は，{A,B},{C,D,E},３クラスターの場合は, {A}, {B},{C,D,E}
であることがわかります．

東京⼤学 数理・情報教育研究センター

Cと(D,E)の距離
AとBの距離DとEの距離AB

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

C

D

E

50

移動距離データの例
デンドログラム作成の説明のために東北・北海
道の県庁所在地間の移動時間が距離として与えられている場合を想定し，最⻑距離法で階層的クラスタリングを⾏ってみます．
① 福島と宮城を併合：距離42分② ⼭形と１の２県を併合：距離76.5分③ ⻘森と岩⼿を併合；距離88分④ 秋⽥と３を併合：距離170分⑤ ２と４を併合：距離251分⑥ 北海道と５を併合：距離373分

北海道

⻘森

岩⼿

宮城

秋⽥

⼭形

北海道

0

⻘森

313.5

0

岩⼿

290.5

88

0

宮城

253.5

146.5

61.5

0

秋⽥

373

170

97

153.5

0

⼭形

319

247.5

161

76.5

251

0

福島

283

193.5

115,5

42

205.5

72.5

福島

0

（移動距離：県庁所在地のJR駅を朝8時に出発した時の到達時間）

結果として，右下の樹形図が得られます．
①〜⑥は併合していった順番を⽰します．
逆に樹形図が与えられると，横線を上からおろしていくと，次のように 2つ〜５つのクラスターが順番に得られます．
2
3
4
5

⑥
⑤
④
②

(北海道） (⻘森，岩⼿，宮城，秋⽥，⼭形，福島）(北海道） (⻘森，秋⽥，岩⼿）（宮城，⼭形，福島）(北海道）（秋⽥）（⻘森，岩⼿）（宮城，⼭形，福島）(北海道）（秋⽥）（⻘森，岩⼿）（⼭形）（宮城，福島）

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

①

③

51

例：多変量データの場合
重回帰モデルの説明変数として⽤いたデータを拡⼤した９変量（下記）のデータに対して階層的クラスタリングを適⽤した計算⼿順と結果を⽰します．
１．９変量のデータをそれぞれ，平均0，分散1に変換して規準化します．
２．９つの基準化したデータの距離を計算します．
３．クラスタリングを実施します．

１．⼈⼝密度
２．⾼齢化率３．県庁所在地の⼈⼝密度
４．交通便利度（4時間）
５．平均気温６．最⾼気温７．最低気温８．平均湿度
９．快晴⽇数東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

52

⾮階層的クラスタリング
階層的クラスタリングは全データバラバラの状態から，⼀つずつクラスターを形成していき，また最終的に樹形図ができるので，直感的にもわかりやすいです．
しかしながら，この⽅法ではまず全データ間の距離を計算する必要があり，またク
ラスター間の距離の定義によっては，毎回距離を更新していく必要があるために，巨⼤なデータに適⽤することは困難です．
⾮階層的クラスタリングの代表的な k-means法では、樹形図（階層的な構造）を想定せず，予め定めたクラスター数について，適当に定めた初期クラスターから繰り返し計算によって，与えられた距離を最⼩にするようにクラスターを求めます．
３クラスターの場合

収束するまで繰り返す

3つの点（⻩⾊）をランダムに配置し，すべての点（⿊）との距離を計算する．
東京⼤学 数理・情報教育研究センター

最短距離を基準に3クラスターに分類し，それぞれのクラスターの中⼼との距離を計算し直す．

同じプロセスを繰り返す．

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

重⼼の位置が変化しなくなったら終了．

53

⾮階層的クラスタリングの例
k =2の場合

k =4の場合

k =3の場合

k =5の場合

階層的クラスタリングでも⽤いた多変量データのうち1, 3, 4, 5, 6番⽬の5変量データにk-means法を適⽤した結果です. kとしては2, 3, 4, 5 の4つの場合を⽰しています．
右図では5次元のうち，第1主成分と第2主成分だけを表⽰しています．
⼤雑把にいうと沖縄を除き, k=2 の場合は⾸都圏
・近畿・中部の⼈⼝密度が⾼い県とそれ以外にグループ化されています．
k = 3ではそれ以外が２つに分けられ，k = 4の場合
は，さらに⼤都市圏が２分割されています．
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

54

54

⾮階層的クラスタリングのまとめ
⾮階層的クラスタリングのメリット• 全サンプル間の距離を計算する必要がなく，仮定した重⼼と全サンプル間の距離だけを計算すればよいので，階層的クラスタリングよりも計算量が少な
く，ビッグデータにも適⽤できます。

⾮階層的クラスタリングのデメリット• 最初にランダムに指定する重⼼の位置によって，最終的結果が変わることがあります（初期値依存性）• 何個のクラスターに分けるのがよいかの明確な判断基準がありません。

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

55

５．パターン発⾒POSデータなどの販売・取引に関するデータから意思決定に役に⽴つ規則や注⽬すべきパターンを発⾒する⽅法としてアソシエーション分析があります．⽀持度，確信度，リフト値などの概念と相関ルール抽出や頻出パターン発⾒の⽅法としてのAprioriアルゴリズムを紹介します.

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

56

56

POSデータからのパターン発⾒
コンビニなどの店舗で買い物をするとレジで様々な情報が⼊⼒されており
POS（Point of Sale）データが得られます．POSデータはチェーン店などを中⼼に販売戦略に活かされています．
POSデータのうち，どのような商品を⼀緒に購⼊したかの部分だけからも重要な情報が得られることがあります．例えば，下記のような販売データがあったとします．このような販売データが⼤量にあったとき，その中から販売促進に役に⽴ちそうな情報を取り出す⽅法がアソシエーション分析
です．バスケット分析と呼ばれることもあります．
ID

バスケット

1

{A, B, C, D}

2{A, B, E}

3{B, C, F, G}

4{D, E, F, G, H}

5{B, C, E, H}

東京⼤学 数理・情報教育研究センター

パン，ジュース，⽜乳などの商品はアイテムと呼ばれます．
⼀⼈のお客はいくつかのアイテムを購⼊するので，それを{ A, B, C, D} のように表してバスケットと呼びます．また，I
I ={A, B}のようにアイテムを集めたものをアイテム集合と呼びます．
多くのIDのバスケットの中に同じアイテム集合（例えば IとJ ）が現れる場合，I とJ は⼀緒に買われる可能性が⾼いことがわかるので，並べて陳列するなどの販売促進策が考えられます．アソシエーション分析によって，紙オムツとビールが頻出アイテムとして⾒出されたというのが有名な例です．

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

57

57

アソシエーション分析
アソシエーション分析では，厖⼤な販売データに⾒られる商品Bを購⼊すると商品Cも買う傾向があるというような，関連性や規則性，相関ルール（association rule）を⾒出すことを⽬的とします．
販売データは以下のような情報を持ちます．

ID

通し番号，何番⽬の購⼊(バスケット)かを n で表す

j (n)n 番⽬の客が購⼊した商品数

In{I1 ,, I j ( n ) } n 番⽬のバスケット内のアイテム

アソシエーション分析では，全バスケットの中でアイテム集合 I が購⼊されたときに，アイテム集合 J が購⼊される傾向があるときに I  Jと表現します．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

58

58

アイテムとアイテム集合
相関ルールではアイテム集合 I と J について， I  J という関係を考えます．
ただし，I と J は ｛A} のような単独のアイテムの場合も {A, B, E} のような複数のアイテムの集合の場合もあります．
相関ルールは I ならば必ず J が起こるというルールではないので，次⾴のように⽀持度，確信度，リフト値などの指標が考慮されます．
これらの定義には以下の記号が使われます．
N
I, J
k(I )

東京⼤学 数理・情報教育研究センター

全データ数（バスケットの数）
アイテム集合
全データ中にアイテム集合 I が現れる回数

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

ID

バスケット

1

{A, B, C, D}

2{A, B, E}

3{B, C, F, G}



・・・

N

{B, C, E, H}

59

59

⽀持度，確信度，リフト値
⽀持度：

全データ中に I と J が同時に現れる割合
S (I  J ) 

確信度：

k (I  J )N

I が出現した時に J も同時に現れる割合
C (I  J ) k (I  J )k (I )

リフト値：⽀持度と全データの中でアイテム 集合 J が出現する割合の⽐
L( I  J ) 

C(I  J ) k (I  J )Nk(J ) / Nk ( I )k ( J ) ⽀持度は全体の中で I と J が同時に現れる割合で，因果関係を意味しているわけではな
いが，シンプルで分かりやすい指標です．
 確信度は I と J の関係以上に J の出現割合に依存する指標なので注意が必要です．
 リフト値が１より⼤きければアイテム集合 J の販売に I の関連性があると考えられます.

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

60

60

例：簡単な仮想データの場合
右図の販売データの場合，N (データ数 )  5, k ( B )  4, k (C )  3, k ( B  C )  3

なので，相関ルール B  C についてk (B  C) 3
S (B  C)   0.6
5Nk (B  C) 3
C(B  C)   0.75
k ( B)
4
L( B  C ) 

ID

アイテム集合1{A, B, C, D}

2{A, B, E}

3{B, C, F, G}

4{D, E, F, G, H}

5{B, C, E, H}

k (B  C)N 3 5

 1.2543
k ( B ) k (C )となり，⽀持度，確信度，リフト値はそれぞれ，0.6，0.75，1.25 となります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

61

61

Aprioriアルゴリズム
販売データの中に頻繁に表れる頻出パターンや相関ルールを求めるためには，すべてのアイテム集合がデータの中に何回現れるかを調べることにな
ります．これは⼀⾒すると簡単そうに⾒えますが，ビッグデータの場合には組み合わせ爆発によって計算が困難になります．Apriori アルゴリズム
はアソシエーション分析において⾼速に相関ルールや頻出アイテム集合を検出するために1990年代前半にIBM研究所で開発されたアルゴリズムで,その後開発されたアルゴリズムの基礎ともなっています．
Aprioriアルゴリズムでは，⽀持度(support)と確信度(confidence) の閾値を指定するとそれ以上の値の相関ルールを列挙します．⽀持度を⼩さく設定しすぎると⼤量のルールを列挙するので注意が必要です．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

62

62

Aprioriアルゴリズムの適⽤例（１）
Rのライブラリーに標準データとして付属している Groceries（⾷料品データ N=9835）に対してAprioriアルゴリズムを適⽤し，リフト値が⼤きな順に並べた結果を⽰します．ただし，support=0.005, confidence=0.6と設定した場合の結果です．
rhs（結論）

⽀持度

確信度

リフト値

⇒ {他の野菜｝

0.00579

0.633

3.273

２

｛果物, 根菜, ⽜乳｝⇒｛他の野菜｝

0.00549

0.614

3.171

3

｛果物, クリーム｝

0.00559

0.604

3.124

4

{根菜, ⽟葱｝

0.00569

0.602

3.112

ID

lhs（条件）⇒１{柑橘類, 根菜, ⽜乳｝⇒｛他の野菜｝⇒ {他の野菜｝

5｛南国果物, 根菜, ヨーグルト｝⇒｛⽜乳｝

0.00569

0.700

2.740

6

｛バター，クリーム｝

0.00549

0.675

2.642

⇒｛⽜乳｝

この結果は，⾷料品だけの間の相関分析で⽬新しいものではありませんが，アメリカにおける相関分析の有名な例として紙オムツとビールの関連性の発⾒がよく知られています．
このように⼈間では気づきにくい関連性の発⾒が，ビッグデータ分析の醍醐味で，それらを並べて陳列することによって販売促進に役に⽴つ可能性があります．
東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

63

63

Aprioriアルゴリズムの適⽤例（２）同じデータでも⽀持度と確信度の閾値を変えると全く違った結果が得られることがあります．
下の結果はsupport=0.001, confidence=0.5と設定した場合です．

ID
１

lhs（条件）rhs（結論）

⽀持度

確信度

リフト値

⇒ {ハンバーグ⾁｝

0.001220.632

18.996⇒{インスタント⾷品，ソーダ｝

２｛ソーダ，ポップコーン｝⇒｛塩味スナック｝

0.001220.632

16.6983｛⼩⻨粉，ベーキングパウダー｝⇒｛砂糖｝

0.001020.556

16.4084｛ハム，プロセスチーズ｝⇒｛⽩パン｝

0.00193

0.633

15.045

5{⽜乳，インスタント⾷品}⇒｛ハンバーグ⾁｝

0.001530.500

15.038

6｛野菜，凝乳，ヨーグルト，クリーム｝⇒｛クリームチーズ｝

0.001020.588

14.834

⽀持度の閾値を下げるとリフト値が15を超えるパターンが検出されるようになります．また前⾴の結果と違って，野菜や乳製品以外のアイテムも出現します．ただし，⽀持度は0.1％とかなり⼩さな値になります．

東京⼤学 数理・情報教育研究センター

4-4. 時系列データの解析
東京⼤学 数理・情報教育研究センター
北川源四郎 2021 CC BY-NC-SA

64

64

1-4-6

次元削減

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

65

動機：多次元のデータの取り扱い• 多次元のデータは，さまざまな場面で現れます．
• たとえば，100 人の「身長，体重，体温，血圧，脈拍，酸素濃度，血糖値」のデータは，7 次元空間の中の 100 個の点の集合とみなせます．
• 高次元のデータは，低次元のデータに（近似的に）変換してから扱うと，さまざまな点で便利です．
• データの可視化が，容易になります．
• 少ない数の因子で，データを説明できるようになります．
• データの容量が小さくなり，計算機で扱いやすくなります．
• このように，与えられたデータをより低次元のデータに近似的に変換することを，次元削減とよびます．
• 元のデータがもつ情報を完全に保ったまま次元削減を行うことは，一般に，不可能です．したがって，元のデータがもつ情報をなるべく保持したまま次元を減らすことが，次元削減の目標です．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

66

2 次元のデータを 1 次元に圧縮する例 (1/4)• 説明の都合上，まず，2 次元のデータを 1 次元に変換する状況を考えましょう．

➔• 次の図のようなデータ集合が与えられたとします：➔• ある特定の視点（方向）からこのデータを見るとき，情報をなるべく失わないためにはどの視点（方向）から見ればよいでしょうか．
• たとえば， の方向から見た場合は，次のように見えます：• これは，𝑥2 の情報を捨てて，𝑥1 のみの情報をとり出したことに相当します．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA67

2 次元のデータを 1 次元に圧縮する例 (2/4)

➙

➔
➔

• 一方，

の方向からみると，次のように見えます：➙• データをの方向からみると，次のように見えます：• 後者では，ほとんどの情報が失われています（データ点どうしの区別が，ほとんどつきません）．
• 前者のように，ばらつき（つまり，分散）が大きくなる方向が，望ましい方向であることがわかります．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

68

2 次元のデータを 1 次元に圧縮する例 (3/4)• 先ほどは極端な例でしたが，このデータ集合ではどうでしょうか．

➙

➔
➔

• 一方，

の方向からみると，次のように見えます：➙• データをの方向からみると，次のように見えます：• やはり，分散が大きくなる前者の方向が，望ましい方向であると考えられます．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

69

2 次元のデータを 1 次元に圧縮する例 (4/4)• 「ある方向からデータを見る」ことは，「データ点の座標の重み和をとる（ただし，重みベクトルのノルムは 1 とする）」ことに相当します．

➙➔➔• 𝑥2 軸から 𝜋/3 回転した方向

から見ることは， 𝑥 1

計算することに相当します．

➙• 𝑥2 軸から −𝜋/6 回転した方向

から見ることは， 𝑥 1を計算することに相当します．







𝑤1
を用いて 𝑥 1• 一般に， k𝒘k = 1 を満たす 𝒘 =𝑤2 cos(𝜋/3)𝑥2をsin(𝜋/3)
 cos(−𝜋/6)𝑥2sin(−𝜋/6)



 
 𝑤1
𝑥2
と表せます．
𝑤2

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA70主成分分析：分散最大化としての導入• 𝑛 次元のデータを 1 次元に圧縮することを考えましょう．
• 𝑚 個のデータ点

𝒙 𝑖 = 𝑥𝑖1 𝑥𝑖2 · · ·が与えられているとします．

𝑥𝑖𝑛

(𝑖 = 1, . . . , 𝑚)

• 以降の表記を簡潔にするため，データ点は行ベクトルで表します．
• 1 次元への圧縮は， k𝒘k = 1 を満たす列ベクトル 𝒘 ∈ R𝑛 を用いて 𝒙 1 𝒘,
𝒙 2 𝒘, . . . , 𝒙 𝑚 𝒘 と表せます．
• 以下では，データの情報をできるだけ失わないため，𝒙 1 𝒘, 𝒙 2 𝒘, . . . , 𝒙 𝑚 𝒘の分散が最大になるように 𝒘 を選ぶことを考えます．
• このような次元削減の手法を，主成分分析とよびます．
• 主成分分析の解説の前に，予備知識として，行列の固有値分解と特異値分解についてまとめておきましょう．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

71

予備知識：対称行列の固有値分解 (1/2)• 𝑛 次の対称行列 𝐴 に対して，𝑛 次の直交行列 𝑄 をうまく選ぶことで𝜆 1 0 · · · 0  0 𝜆2 · · · 0  ⊤
𝐴=𝑄 .
..  𝑄.. . .

 ..
.
.
.


 0 0 · · · 𝜆𝑛


|
{z
}
対角行列

と表すことができます．これを， 𝐴 の固有値分解とよびます．

• ただし，𝑛 次の正方行列 𝑄 が 𝑄 ⊤ 𝑄 = 𝐼 （つまり，𝑄 −1 = 𝑄 ）を満たすと
き，𝑄 を直交行列とよびます．
• また，𝜆 1 , 𝜆 2 , . . . , 𝜆 𝑛 は 𝐴 の固有値です（➪ 1-6-5 線形代数）．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

72

予備知識：対称行列の固有値分解 (2/2)• 𝐴 の固有値分解において，𝑄 を𝜆 1 0 · · · 0 𝜆2 · · ·
𝐴=𝑄 .
.. . .
 ..
.
.

0 0 ···


0 
0 
..  𝑄 ⊤. 
𝜆 𝑛 
𝑄 =  𝒒 1 𝒒 2 · · · 𝒒 𝑛 と表すと，𝑄 の列ベクトル 𝒒 𝑗 ( 𝑗 = 1, 2, . . . , 𝑛) は固有値 𝜆 𝑗 に対応する 𝐴の固有ベクトルです（➪ 1-6-5 線形代数）．
• 𝑄 が直交行列であることから， 𝒒 𝑗 ( 𝑗 = 1, 2, . . . , 𝑛) は， k𝒒 𝑗 k = 1 , 𝒒 ⊤𝑗 𝒒 𝑙 = 0( 𝑗 ≠ 𝑙) を満たします．
• このように，ノルムが 1 で，かつ，互いに直交するベクトルの組を，正規直交系とよびます．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

73

予備知識：行列のランク• ベクトル 𝒙 1 , 𝒙 2 , . . . , 𝒙 𝑛 が与えられたとき，定数 𝑎 1 , 𝑎 2 , . . . , 𝑎 𝑛 に対して
𝑎1 𝒙1 + 𝑎2 𝒙2 + · · · + 𝑎 𝑛 𝒙 𝑛 = 0が成り立つのが 𝑎 1 = 𝑎 2 = · · · = 𝑎 𝑛 = 0 である場合に限られるとき，𝒙 1 ,𝒙 2 , . . . , 𝒙 𝑛 は線形独立（または，1 次独立）であるといいます．
• 𝑚 × 𝑛 型の実行列 𝐴 を考えます．
• 𝐴 の 𝑛 本の列ベクトルのうち線形独立なものの最大個数を， 𝐴 のランク
（または，階数）とよびます．
• 𝐴 のランクは，また， 𝐴 の 𝑚 本の行ベクトルのうち線形独立なものの最大個数でもあります．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

74

予備知識：特異値分解 (1/2)• 𝑚 × 𝑛 型の実行列 𝐴 のランクを 𝑟 で表します．𝑚 次の直交行列 𝑈 と 𝑛 次の直交行列 𝑉 をうまく選ぶことで， 𝐴 を
𝐴 = 𝑈 𝛴𝑉 ⊤と表すことができます．ただし， 𝛴 は 𝑚 × 𝑛 型の実行列で，
𝜎1 0 · · · 0 0 · · · 0
 0 𝜎2 · · · 0 0 · · · 0 ..
.. .. . .
.. .. . .
.
. . .
. .
.



𝛴 =  0 0 · · · 𝜎𝑟 0 · · · 0 , 𝜎1 ≧ 𝜎2 ≧ · · · ≧ 𝜎𝑟 > 0
 0 0 · · · 0 0 · · · 0


 ..

.
.
.
.
.
.
.
.
.
.
.
.
.
. . 
. . .
.

 0 0 · · · 0 0 · · · 0
の形です．これを， 𝐴 の特異値分解とよびます．
• ここで，𝜎1 , 𝜎2 , . . . , 𝜎𝑟 を 𝐴 の特異値とよびます．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

75

予備知識：特異値分解 (2/2)• 𝐴 の特異値分解を考えます： 𝐴 = 𝑈 𝛴𝑉 ⊤ .
• 𝐴 の特異値は， 𝐴⊤ 𝐴 の（ゼロでない）固有値の平方根に等しいです．
• 𝑉 = 𝒗 1 𝒗 2 · · · 𝒗 𝑛 と表すと，𝑉 の列ベクトル 𝒗 𝑗 は 𝐴⊤ 𝐴 の固有ベクトルです（この 𝒗 𝑗 を， 𝐴 の右特異ベクトルとよびます）．
• 𝐴 の特異値は， 𝐴𝐴⊤ の（ゼロでない）固有値の平方根に等しいです．
• 𝑈 = 𝒖 1 𝒖 2 · · · 𝒖 𝑚 と表すと，𝑈 の列ベクトル 𝒖 𝑖 は 𝐴𝐴⊤ の固有ベクトルです（この 𝒖 𝑖 を， 𝐴 の左特異ベクトルとよびます）．
• 行列 𝐴 のノルム（フロベニウスノルム）は，
ppk 𝐴k F = tr( 𝐴⊤ 𝐴) = tr( 𝐴𝐴⊤ )で定義されます（正方行列 𝐵 に対して，対角成分の和をトレースとよび，tr 𝐵 で表します）．
• 𝐴 の特異値を用いると， pk 𝐴k F = 𝜎1 2 + 𝜎2 2 + · · · + 𝜎𝑟 2と表せます．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

76主成分分析：1 次元への圧縮 (1/2)• データ行列 𝑋 ∈ R𝑚×𝑛 を次式で定義します：


𝒙1




𝒙2


𝑋=
 .
..


.


 𝒙𝑚




• 各 𝑗 = 1, 2, . . . , 𝑛 に対して，𝑥1 𝑗 , 𝑥2 𝑗 , . . . , 𝑥 𝑚 𝑗 の平均を 𝑥¯ 𝑗 で表します：𝑚1 Õ𝑥¯ 𝑗 =
𝑥𝑖 𝑗 .
𝑚 𝑖=1• 平均 𝑥¯ 𝑗 を並べた行ベクトルを 𝒙¯ で表します：𝒙¯ = 𝑥¯1 𝑥¯2 · · ·𝑥¯ 𝑛 =

1 ⊤1𝑚 𝑋 ,𝑚• 𝑋 の中心化 𝑋0 ∈ R𝑚×𝑛 を次式で定義します：
𝑋0 = 𝑋 − 1𝑚 𝒙¯ .

1
 
 
1𝑚 =  ...  .
 
1
 

1 ⊤
𝑋0 は，データ行列の平均を 0 としたものです： 1𝑚 𝑋0 = 𝑂 1,𝑛 .
𝑚
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

77主成分分析：1 次元への圧縮 (2/2)• 𝒙 1 𝒘, . . . , 𝒙 𝑚 𝒘 の平均は，
𝑚
𝑚

1 Õ
1 Õ
𝒙𝑖 𝒘 =
𝒙 𝑖 𝒘 = 𝒙¯ 𝒘 .
𝑚 𝑖=1
𝑚 𝑖=1

• 𝒙 1 𝒘, . . . , 𝒙 𝑚 𝒘 の分散は，
𝑚
𝑚
1

Õ

1
1 Õ
2
(𝒙 𝑖 − 𝒙¯ )𝒘 = 𝒘 ⊤𝑋0⊤ 𝑋0 𝒘 .
(𝒙 𝑖 𝒘 − 𝒙¯ 𝒘) 2 =
𝑚 𝑖=1
𝑚 𝑖=1
𝑚

1 ⊤
• したがって，𝑆 = 𝑋0 𝑋0 とおくと，求めたいものは，条件 k𝒘k = 1 を満𝑚⊤
たす 𝒘 のうち 𝒘 𝑆𝒘 を最大にするものです．
• 𝑆 は，標本分散共分散行列とよばれる 𝑛 次の対称行列です．

• そのような 𝒘 は，𝑆 の最大固有値に対応する単位固有ベクトルであることがわかります（スライド80で説明します）．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

78主成分分析： 𝑝 次元への圧縮 (1/3)• 以上では，𝑛 次元のデータの 1 次元への圧縮 𝒙 1 𝒘, . . . , 𝒙 𝑚 𝒘 について，圧縮後の分散が最大になるように単位ベクトル 𝒘 を選ぶことを考えました．
• 以下では，これを一般化して，データの 𝑝 次元 ( 𝑝 < 𝑛) への圧縮を考えます．これは，分散ができるだけ大きくなるように正規直交系 𝒘 1 , . . . , 𝒘 𝑝
を選ぶことに相当します．
• 正規直交系は，スライド 73 で定義しました．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

79主成分分析： 𝑝 次元への圧縮 (2/3)• 標本分散共分散行列 𝑆 は 𝑛 次の対称行列ですから，その固有ベクトル𝒗 1 , . . . , 𝒗 𝑛 を正規直交系としてとると𝒗⊤ 𝜆 1 0 · · · 0  

1



⊤



𝒗2
 𝒗 1 𝒗 2 · · · 𝒗 𝑛   0 𝜆2 · · · 0  
. . .
𝑆 = 



..
..
. .
.



. .
.
. .




⊤



  0 0 · · · 𝜆𝑛 𝒗𝑛と表せます．ただし，実数 𝜆 1 ≧ · · · ≧ 𝜆 𝑛 は 𝑆 の固有値です．

• 𝒙 1 𝒘, . . . , 𝒙 𝑚 𝒘 の分散は 𝒘 ⊤ 𝑆𝒘 でした．

• 𝒘 = 𝑐 1 𝒗 1 + · · · + 𝑐 𝑛 𝒗 𝑛 とおけば，𝒘 ⊤ 𝑆𝒘 = 𝜆 1 𝑐 1 2 + · · · + 𝜆 𝑛 𝑐 𝑛 2 です．また，
条件 k𝒘k = 1 は 𝑐 1 2 + · · · + 𝑐 𝑛 2 = 1 と書けます．
• つまり，分散 𝒘 ⊤ 𝑆𝒘 は，𝜆 1 , . . . , 𝜆 𝑛 の非負の重み付き和（重みの和をとると 1）とみなせますから，𝒘 = 𝒗 1 で最大値をとることがわかります．
• 次に，𝒗 1 に直交するベクトルは 𝒘 = 𝑐 2 𝒗 2 + · · · + 𝑐 𝑛 𝒗 𝑛 と表せます．このとき，𝒘 ⊤ 𝑆𝒘 = 𝜆 2 𝑐 2 2 + · · · + 𝜆 𝑛 𝑐 𝑛 2 は 𝒘 = 𝒗 2 において最大値をとります．
• つまり， 𝑝 次元への圧縮には 𝒗 1 , . . . , 𝒗 𝑝 を用いればよいことがわかります．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

80主成分分析： 𝑝 次元への圧縮 (3/3)• このようにして，データ行列 𝑋 の 1 次元への圧縮は𝑋𝒗 1で与えられ，2 次元への圧縮は





𝑋 𝒗 1 𝒗 2 で与えられます．一般に， 𝑝 次元への圧縮は𝑋 𝒗 1 𝒗 2 · · · 𝒗 𝑝 で与えられます．

• 𝑋𝒗 1 の第 𝑖 行 𝒙 𝑖 𝒗 1 を 𝒙 𝑖 の第 1 主成分とよび， 𝑋𝒗 2 の第 𝑖 行 𝒙 𝑖 𝒗 2 を 𝒙 𝑖 の
第 2 主成分とよびます．一般に，𝒙 𝑖 𝒗 𝑝 を 𝒙 𝑖 の第 𝑝 主成分とよびます．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

81

低ランク行列近似としての主成分分析 (1/3)
1 ⊤• 主成分分析で用いるベクトル 𝒗 1 , . . . , 𝒗 𝑛 は，行列 𝑆 = 𝑋0 𝑋0 の正規直交𝑚な固有ベクトルでした．
• したがって，𝒗 1 , . . . , 𝒗 𝑛 は中心化されたデータ行列 𝑋0 ∈ R𝑚×𝑛 の右特異ベクトルでもあります．特異値分解により 𝑋0 は






  𝒗⊤
1




.
..
𝑋0 = 𝜎1 𝒖 1 · · · 𝜎𝑟 𝒖𝑟 𝑂 𝑚,𝑛−𝑟  





⊤


  𝒗𝑛と表せます．ただし，𝑟 は 𝑋0 のランクであり，𝜎1 ≧ · · · ≧ 𝜎𝑟 > 0 は 𝑋0 の特異値であり，𝒖 1 , . . . , 𝒖𝑟 ∈ R𝑚 はこれに対応する（正規直交な）左特異ベ
クトルです．
• ここで 𝑋0 𝒗 𝑗 = 𝜎 𝑗 𝒖 𝑗 が成り立ちますから，𝜎 𝑗 𝒖 𝑗 ( 𝑗 = 1, . . . , 𝑟) は中心化されたデータ 𝒙 1 − 𝒙¯ , . . . , 𝒙 𝑚 − 𝒙¯ の第 𝑗 主成分を縦に並べたベクトルである
ことがわかります．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

82

低ランク行列近似としての主成分分析 (2/3)• スライド82の結果から， 𝑋0 ∈ R𝑚×𝑛 は







 
  𝒗⊤


1
 




.



..
𝑋0 = 𝜎1 𝒖 1 · · · 𝜎𝑝 𝒖 𝑝   + 𝜎𝑝+1 𝒖 𝑝+1 · · · 𝜎𝑟 𝒖𝑟  

 
 


⊤
 
  𝒗𝑝
 










|
{z

= 𝐸 とおきます．

𝒗 ⊤𝑝+1
..
.
𝒗 𝑟⊤









}

と表せます．ここで，もし行列 𝐸 のノルムが小さければ，𝑛 次元のデータ
𝑋0 はランクが 𝑝 の行列で近似できることになります：𝑚 × 𝑛 型で，ランクは 𝑝 (< 𝑛)

z
}|


 

 

 


𝑋0 ≃ 𝜎1 𝒖 1 · · · 𝜎𝑝 𝒖 𝑝  
 

 





|
{z
}
𝑋0 の 𝑝 次元への圧縮

𝒗⊤
1
..
.
𝒗 ⊤𝑝

{




 .




これを，𝑋0 の低ランク行列近似とよびます．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

83

低ランク行列近似としての主成分分析 (3/3)• スライド83で，行列  𝒗⊤




𝑝+1




..
𝐸 = 𝜎𝑝+1 𝒖 𝑝+1 · · · 𝜎𝑟 𝒖𝑟  

.




⊤



𝒗𝑟





のノルムが小さいほど，低ランク行列近似の精度はよいと考えられます．
𝐸 のフロベニウスノルム k𝐸 k F の 2 乗はk𝐸 k 2F = tr(𝐸 ⊤ 𝐸) = 𝜎𝑝+1 2 + · · · + 𝜎𝑟 2
です．
• 同様に， k 𝑋0 k 2F = 𝜎1 2 + · · · + 𝜎𝑟 2 が成り立ちます．
• したがって，近似の精度の尺度としてk 𝑋0 k 2F − k𝐸 k 2F
𝜎1 2 + · · · + 𝜎𝑝 2
= 22
𝜎1 + · · · + 𝜎𝑝 2 + · · · + 𝜎𝑟 2
k 𝑋0 k Fが得られますが，この値を累積寄与率とよびます．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

84非線形の次元削減：導入 (1/2)• 主成分分析は，「データを見る方向」を上手に選ぶことで，分散をなるべく保ったまま次元を削減する手法でした．
「見る方向」を定めることは，主成分分析が本質的に線形変換を行っていることを意味します．
• このため，主成分分析は，たとえば左図のデータならうまく扱うことができます．
• 右図のデータ（スイスロールとよばれます）は，どうでしょうか．
• これは 3 次元のデータですが，データ点は（2 次元的な）曲面上に載っているように見えます．

x3
➔

x2
x1
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

85非線形の次元削減：導入 (2/2)• 右図の「スイスロール」の曲面を（ある方法で）2 次元平面上に広げると，
左図のようになります．これで，次元が 1 つ削減できました．
• 曲面を平面に広げるのは非線形な変換ですので，主成分分析では扱うことができません．
• 以下では，非線形な次元削減の手法として，カーネル主成分分析をとりあげます．
• その解説の前に，予備知識として，カーネル法の概要をみておきましょう．

x3

←
x2
x1
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

86

カーネル法：回帰分析への適用例 (1/3)• 線形回帰（左下の図）では，目的変数 𝑦 ∈ R を説明変数 𝒙 ∈ R𝑛 の 1 次関
数として近似します．
• 与えられたデータ点を (𝒙 1 , 𝑦 1 ), (𝒙 2 , 𝑦 2 ), . . . , (𝒙 𝑚 , 𝑦 𝑚 ) とします．
• 𝒙 1 , 𝒙 2 , . . . , 𝒙 𝑚 は，これまでと同様に，行ベクトルとします．
• 回帰式を 𝑦 = 𝒙𝜶 + 𝛽 とします（𝜶 = (𝛼1 , 𝛼2 . . . , 𝛼𝑛 ) ⊤ です）．
• 最小 2 乗法では，回帰式の予測とデータ 𝑦 𝑖 の差の 2 乗和𝑚Õ2(𝒙𝜶 + 𝛽) − 𝑦 𝑖 = k 𝑋𝜶 + 𝛽1𝑚 − 𝒚k 2 が最小になるように 𝜶 ∈ R𝑛 と𝑖=1

𝛽 ∈ R を選びます．

20
10
0
-10
-20
-2

0

2

• 𝑋 と 𝒚 の定義は次のとおりです：
 𝑦1 


𝒙1
 


 𝑦2 


𝒙2
 


,
𝒚
=
𝑋=
 ..  .

.
.
 . 


.
 


𝑦 𝑚

 𝒙𝑚
 



東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

87

カーネル法：回帰分析への適用例 (2/3)• 右下の図のデータは，1 次関数では当てはまりがとても悪そうです．このような場合に用いられる方法の 1 つがカーネル法です．
• カーネル法では，カーネル関数とよばれる関数を用います．
• 代表例に，ガウスカーネルがあります（𝛾 > 0 はパラメータです）：
𝑘 (𝒙, 𝒔) = exp(−𝛾k𝒙 − 𝒔k 2 ) .
𝑚
Õ• 回帰式は 𝑦 =𝜏𝑙 𝑘 (𝒙 𝑙 , 𝒙) とします．そして，回帰式の予測とデータの差の 2 乗和

𝑙=1
𝑚
𝑚
Õ hÕ

𝜏𝑙 𝑘 (𝒙 𝑙 , 𝒙 𝑖 ) − 𝑦 𝑖𝑖=1 𝑙=1

を選びます．
20i2が最小になるように 𝜏1 , 𝜏2 , . . . , 𝜏𝑚 ∈ R

20

10

0

0
-20
-10
-40
-20
-2

0

2

-2

0

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

2

88

カーネル法：回帰分析への適用例 (3/3)• 関数 𝑓 (𝝉) =

𝑚 hÕ
𝑚
Õ
𝑖=1 𝑙=1

𝜏𝑙 𝑘 (𝒙 𝑙 , 𝒙 𝑖 ) − 𝑦 𝑖i2を最小にする 𝝉 = (𝜏1 , 𝜏2 , . . . , 𝜏𝑚 ) ⊤は，次のようにして求められます．

• まず，𝑚 次の対称行列 𝐾 を
 𝑘 (𝒙 1 , 𝒙 1 ) 𝑘 (𝒙 1 , 𝒙 2 ) · · · 𝑘 (𝒙 1 , 𝒙 𝑚 )  𝑘 (𝒙 2 , 𝒙 1 ) 𝑘 (𝒙 2 , 𝒙 2 ) · · · 𝑘 (𝒙 2 , 𝒙 𝑚 ) 


𝐾=

..
..
..
..


.
.
.
.


 𝑘 (𝒙 𝑚 , 𝒙 1 ) 𝑘 (𝒙 𝑚 , 𝒙 2 ) · · · 𝑘 (𝒙 𝑚 , 𝒙 𝑚 ) 
で定めます．この 𝐾 を，グラム行列とよびます．

• カーネル関数 𝑘 は，グラム行列が常に正則であるように定めます（たとえばガウスカーネルはこの条件を満たしています）．

• この 𝐾 を用いると， 𝑓 (𝝉) = k𝐾𝝉 − 𝒚k 2 (≧ 0) と書けます．
• したがって， 𝑓 (𝝉) は 𝝉 = 𝐾 −1 𝒚 において最小値 0 をとることがわかります．
• 実際には過学習を防ぐために何らかの正則化を施しますが，以上がカーネル法を回帰分析へ適用する際の概要です．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

89

カーネル主成分分析 (1/2)• それでは，カーネル主成分分析を説明します．主成分分析では，中心化されたデータ行列 𝑋0 の特異値分解を考えました（スライド82）
：






  𝒗⊤
1




.


.
𝜎
𝒖
·
·
·
𝜎
𝒖
𝑂
𝑋0 = 𝑈 𝛴𝑉 =  1 1
.
𝑟 𝑟 𝑚,𝑛−𝑟  
.
(∗1)






  𝒗⊤𝑛
ここで，𝜎 𝑗 𝒖 𝑗 ( 𝑗 = 1, . . . , 𝑟) は 𝒙 1 − 𝒙¯ , . . . , 𝒙 𝑚 − 𝒙¯ の第 𝑗 主成分を縦に並べたベクトルでした．
• また， 𝑋0 の定義（スライド 77）は1 ⊤
𝑋0 = 𝑋 − 1𝑚 𝒙¯ , 𝒙¯ = 1𝑚 𝑋
𝑚
1でしたから，𝐻𝑚 = 𝐼 𝑚 − 1𝑚 1⊤
𝑚 とおくと 𝑋0 = 𝐻 𝑚 𝑋 と書けます．
𝑚• 次に，式 (∗1) より 𝑋0 𝑋0⊤ = 𝑈 𝛴 2𝑈 ⊤ が成り立ちますから，行列⊤ の（ゼロでない）固有値は 𝜎 2 , . . . , 𝜎 2 であり，対𝑋0 𝑋0⊤ = 𝐻𝑚 (𝑋 𝑋 ⊤ )𝐻𝑚1𝑟
応する固有ベクトルは 𝒖 1 , . . . , 𝒖𝑟 です．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

90

カーネル主成分分析 (2/2)⊤ を書き下すと• ここで行列 𝐻𝑚 (𝑋 𝑋 ⊤ )𝐻𝑚
 h𝒙 1 , 𝒙 1 i h𝒙 1 , 𝒙 2 i · · · h𝒙 2 , 𝒙 1 i h𝒙 2 , 𝒙 2 i · · ·⊤ = 𝐻 𝐻𝑚 (𝑋 𝑋 ⊤ )𝐻𝑚
..
..
𝑚 
..

.
.
.

 h𝒙 𝑚 , 𝒙 1 i h𝒙 𝑚 , 𝒙 2 i · · ·
です．

h𝒙 1 , 𝒙 𝑚 i 
h𝒙 2 , 𝒙 𝑚 i 
⊤
 𝐻𝑚
..

.

h𝒙 𝑚 , 𝒙 𝑚 i 

• カーネル主成分分析では，データ点どうしの内積 h𝒙 𝑖 , 𝒙 𝑗 i をカーネル関数⊤ は 𝐻 𝐾 𝐻 ⊤ に置きの値 𝑘 (𝒙 𝑖 , 𝒙 𝑗 ) で置き換えます．すると，𝐻𝑚 (𝑋 𝑋 ⊤ )𝐻𝑚𝑚𝑚
換わります（𝐾 は，グラム行列です）．
⊤ の（ゼロでない）固有値を 𝜆˜ , . . . , 𝜆˜ で表し，対応• したがって，𝐻𝑚 𝐾 𝐻𝑚
1𝑟
する固有ベクトルを 𝒖˜ 1 , . . . , 𝒖˜ 𝑟 で表すと，カーネル主成分分析による 𝑋の 𝑝 次元への圧縮は


p

q
 ˜ p˜ 𝜆 1 𝒖˜ 1 𝜆 2 𝒖˜ 2 · · · 𝜆˜ 𝑝 𝒖˜ 𝑝 






であることがわかります．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA91

1-4-7

最適化概論

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

92最適化とは：最適化問題の定義• 最適化問題とは，与えられた条件の下である関数の値を最小にする解（または最大にする解）を求める問題のこと（その解を最適解とよびます）．
• 最適化とは，何らかの意思決定を支援するために，最適化問題を定式化し，それを解く方法論のこと．
• 最適化問題の構成要素• （決定）変数：意思決定において値を決めるべき量• 目的関数：最小（または最大）にしたい評価尺度となる関数• 制約条件：決定変数が満たすべき条件
𝑌 = 𝑓 (𝑥)

𝑌

𝑎

𝑐

𝑏

𝑥

1 変数の最適化問題の例• 左の最小化問題の例では，決定変数は 𝑥 ∈ R , 目的関数は 𝑓 : R → R , 制
約条件は 𝑎 ≦ 𝑥 ≦ 𝑏 , 最適解は 𝑥 = 𝑐です．
• 最適解における目的関数の値を，最適値とよびます．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

93最適化とは：局所最適解と大域的最適解• ( 𝑓 (𝑥) を最大化する問題の最適解) = (− 𝑓 (𝑥) を最小化する問題の最適解) .
• したがって，最小化問題に限って議論してもよいことがわかります．
• 制約条件をもたない最適化問題は，無制約最適化問題とよばれます．
• 左の最小化問題の例で，10• 点 × のように，その点の近傍には目的関数値がより小さい解が存在しない点を，局所最適解とよびます．

5
0
40

×

∗

30

40

20

30
10

20
10

• 最適解 ∗ のことを他の局所最適解と特に区別したい場合は，大域的最適解とよびます．

2 変数の無制約最適化問題の例
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

94最適化とは：変分問題と最適化問題• 最速降下線問題：• 鉛直面内に曲線 𝑋 = 𝑥(𝑧) があり，端点 (𝑎, 0),(𝑏, ℎ) の位置は指定されています．質点を，重力に従って曲線上を降下させます．

(𝑎, 0)• 質点が降下に要する時間が最小の曲線は？

𝑥(𝑧)• これは，変分問題とよばれる問題の例です．

𝑋

𝑧

𝑧

(𝑏, ℎ)• 下図のように，折れ線で考えてみましょう．
• 各折れ点の高さを決めればよいので，決定変数の数は折れ点の数（図では 9 個）です．
• 逆に言うと，上図の曲線の問題は，無限個の点の高さを決める最適化問題とみなせます．

𝑥1
𝑥2
𝑥9

• このように，変分問題は，無限個の決定変数をもつ最適化問題とみなせます．理工学におけるさまざまな支配式が，変分問題として記述できます．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

95

離散最適化問題：最短路問題を例として• 左図において，点 𝑠 から点 𝑡 に移動する経路のうち，経路上の数字（重みとよびます）の和が最小のものを求める問題を考えます．
s

s
52

20
30

30

51

51

12

58

58

52

20

12

18

18
7

16

17

34

34

19

17

16

7

25
23

19

25

t

辺が 14 本ある例題

23

t

最適解

• この最適化問題は，2 つの点を結ぶ矢印（辺とよびます）それぞれを経路として採用するかしないかを選択する問題とみることができます．この選択は離散的ですので，このような問題を離散最適化問題とよびます．
• 点 𝑠 から点 𝑡 に移動する経路は，有限個しかありません．このように，制約条件を満たす解の個数が有限個の問題を組合せ最適化問題とよびます．
• 離散最適化問題と組合せ最適化問題は，厳密に区別しないことが多いです．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

96

連続最適化問題：線形計画問題を例として• 制約条件
5𝑥1 + 4𝑥2 ≦ 80 ,の下で，目的関数

2𝑥1 + 4𝑥2 ≦ 40 ,2𝑥1 + 8𝑥2 ≦ 64

20𝑥1 + 60𝑥2を最大化する最適化問題を考えます．
• この問題は，数式で次のように記述します：
Maximize

20𝑥1 + 60𝑥2

𝑥2

subject to 5𝑥1 + 4𝑥2 ≦ 80 ,2𝑥1 + 4𝑥2 ≦ 40 ,2𝑥1 + 8𝑥2 ≦ 64 .

𝑥1

• 決定変数 𝑥1 および 𝑥2 は連続的に変化する量ですので，このような問題を連続最適化問題とよびます．
• 特に，この問題は，目的関数も制約条件もすべて 1 次式で表されています．
このような問題を線形計画問題とよびます．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

97最適化の具体例：回帰分析 (1/3)• 回帰分析の 1 つである線形回帰（左下の図）では，目的変数 𝑦 ∈ R を説明
変数 𝑥 ∈ R の 1 次関数として近似します．
• 与えられたデータ点を (𝑥1 , 𝑦 1 ), (𝑥2 , 𝑦 2 ), . . . , (𝑥 𝑚 , 𝑦 𝑚 ) とします．
• 回帰式を 𝑦 = 𝛽0 + 𝛽1 𝑥 とします．
• 最小 2 乗法では，回帰式による予測とデータの差の 2 乗和𝑚Õ2
(𝛽0 + 𝛽1 𝑥𝑖 ) − 𝑦 𝑖 を最小にする 𝛽0 と 𝛽1 の値を求めます．
𝑖=1

20
10

0
-10
-20
-2

0

2

• 𝑋, 𝜷 𝒚を
 𝑦1 
1 𝑥1    𝑦2 
1 𝑥2 𝛽0 , 𝒚= . 
𝑋 = . .  , 𝜷 =𝛽1 .. 
 .. .. 
 


𝑦 𝑚
1 𝑥 𝑚  で定義すると，目的関数は k 𝑋 𝜷 − 𝒚k 2 と書けます（ 𝜷 が，最適化の決定変数です）．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA98最適化の具体例：回帰分析 (2/3)• 重回帰分析では，説明変数の数を 𝑝 とすると，各 𝑖 = 1, 2, . . . , 𝑚 に対して，
データ点は (𝑥 1𝑖 , 𝑥 2𝑖 , . . . , 𝑥 𝑝𝑖 ) と 𝑦 𝑖 の組として与えられます．
• 回帰式を 𝑦 = 𝛽0 + 𝛽1 𝑥1 + 𝛽2 𝑥2 + · · · + 𝛽 𝑝 𝑥 𝑝 とします．
• 最小 2 乗法では，回帰式による予測とデータとの差の 2 乗和𝑚Õ2
(𝛽0 + 𝛽1 𝑥1𝑖 + 𝛽2 𝑥2𝑖 + · · · + 𝛽 𝑝 𝑥 𝑝𝑖 ) − 𝑦 𝑖𝑖=1を最小にする 𝛽0 , 𝛽1 , 𝛽2 , . . . , 𝛽 𝑝 の値を求めます．

• 𝑋, 𝜷 𝒚を
 𝛽0 
 𝑦1   
 𝛽1 
 𝑦2 
  
 𝛽2 
𝜷=  , 𝒚= . 
 .. 
 .. 
 
 . 
𝑦 𝑚
 
 
𝛽 𝑝 
 
で定義すると，目的関数は k 𝑋 𝜷 − 𝒚k 2 と書けます．さらにこれを展開すると，2 次関数 𝜷⊤ (𝑋 ⊤ 𝑋) 𝜷 − 2(𝑋 ⊤ 𝒚) ⊤ 𝜷 を最小化する 𝜷 を求めればよいことがわかります（定数項 k 𝒚k 2 は最適化に影響しないことに注意します）．
1 𝑥11 𝑥21 · · ·1 𝑥12 𝑥22 · · ·
𝑋 = .
..
..
..
 ..
.
.
.

1 𝑥1𝑚 𝑥2𝑚 · · ·𝑥 𝑝1 
𝑥 𝑝2 
..  ,. 𝑥 𝑝𝑚 

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA99最適化の具体例：回帰分析 (3/3)• データへの過学習（過度なあてはまり）を防ぐための種々の方法が知られています．
• ティコノフ正則化付き最小 2 乗法（リッジ回帰）とよばれる方法では，誤算の 2 乗和に 𝛾k 𝜷k 2 を加えます（𝛾 > 0 はパラメータです）．目的関数を整理すると
𝜷⊤ (𝑋 ⊤ 𝑋 + 𝛾𝐼) 𝜷 − 2(𝑋 ⊤ 𝒚) ⊤ 𝜷となりますので，やはり 2 次関数の最小化問題であることがわかります．
• 一方，ℓ1 ノルム正則化付き最小 2 乗法とよばれる方法では，
𝛾(|𝛽0 | + |𝛽1 | + |𝛽2 | + · · · + |𝛽 𝑝 |) を加えます．この問題は，補助的な変数として 𝑠 𝑗 ( 𝑗 = 0, 1, 2, . . . , 𝑝) を導入すると，1 次不等式で表される制約条件
𝑠 𝑗 ≧ 𝛽 𝑗 , 𝑠 𝑗 ≧ −𝛽 𝑗 ( 𝑗 = 0, 1, 2, . . . , 𝑝)の下で 2 次関数𝑝
Õ𝜷⊤ (𝑋 ⊤ 𝑋) 𝜷 − 2(𝑋 ⊤ 𝒚) ⊤ 𝜷 +𝑠𝑗
𝑗=0を最小化する問題に帰着できます（最適解では各 𝑗 = 0, 1, 2, . . . , 𝑝 に対して 𝑠 𝑗 = |𝛽 𝑗 | が成り立つことに注意します）．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

100

ナップサック問題とは• 品物が 𝑛 個あり，そのうちのいくつかを選んで 1 つのナップサックに詰め込むことを考えます．
• それぞれの品物には，重量と価格が与えられています．ナップサックには重量制限が与えられており，詰め込んだ品物の重量の和がその制限を超えてはいけません．
• このとき，ナップサックに詰め込んだ品物の価格の和が最大になるように品物を選ぶ問題を，ナップサック問題とよびます．
• ナップサック問題は，それぞれの品物をナップサックに入れるか入れないかを選択する問題とみることができます．この選択は離散的ですので，
ナップサック問題は離散最適化問題です．

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

101

整数計画問題としての定式化• 品物の数が 𝑛 = 5 個の場合の具体例を考えます：
品物 1 品物 2 品物 3 品物 4
重量
価格

品物 5

3
4
1
2
33
32
10
28
ナップサックの重量制限：7545• 品物 𝑗 ( 𝑗 = 1, 2,(. . . , 5) に対して，決定変数 𝑥 𝑗 を
1 (品物 𝑗 をナップサックに入れるとき),
𝑥𝑗 =0 (品物 𝑗 をナップサックに入れないとき)と対応させます．すると，ナップサック問題は次のように定式化できます：
Maximize

33𝑥1 + 32𝑥2 + 10𝑥3 + 28𝑥4 + 45𝑥5

subject to 3𝑥1 + 4𝑥2 + 𝑥3 + 2𝑥4 + 5𝑥2 ≦ 7 ,𝑥 𝑗 = 0 または 1 ,𝑗 = 1, 2, . . . , 5 .

• 決定変数が整数値のみをとるため，このような最適化問題は整数計画問題とよばれます．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

102

貪欲算法
• 重量が軽くて価格が高い品物と，重量が重くて価格が安い品物とが，あったとします．ナップサック問題では，前者を詰め込む方が，後者を詰め込
むよりも，得策のように思えます．
• そこで，単位重量あたりの価格が大きい順に，できるだけ多くの品物をナップサックに詰め込むことにすると，(𝑥1 , 𝑥2 , 𝑥3 , 𝑥4 , 𝑥5 ) = (1, 0, 1, 1, 0)が得られます．
品物 1 品物 2 品物 3 品物 4 品物 5
重量
価格
価格/重量

3
4
1
2
33
32
10
28
11
8
10
14
ナップサックの重量制限：7545
9• このような解法を，貪欲算法とよびます．
• 貪欲算法では，必ずしも最適解は得られません．実際，この問題では(𝑥1 , 𝑥2 , 𝑥3 , 𝑥4 , 𝑥5 ) = (0, 0, 0, 1, 1)が最適解です．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

103

緩和問題
• ナップサック問題の「𝑥 𝑗 は 0 または 1 をとる」という制約条件を，不等式「0 ≦ 𝑥 𝑗 ≦ 1」に置き換えてみます：
Maximize

33𝑥1 + 32𝑥2 + 10𝑥3 + 28𝑥4 + 45𝑥5

subject to 3𝑥1 + 4𝑥2 + 𝑥3 + 2𝑥4 + 5𝑥2 ≦ 7 ,0 ≦ 𝑥𝑗 ≦ 1,(∗1)

𝑗 = 1, 2, . . . , 5 .

• このように制約条件を緩めることで得られる問題のことを，緩和問題とよびます．
• 問題 (∗1) の最適解は，貪欲算法で得られることを示すことができます：(𝑥1 , 𝑥2 , 𝑥3 , 𝑥4 , 𝑥5 ) = (1, 0, 1, 1, 1/5)
品物 1
重量
価格
価格/重量

品物 2

品物 3

品物 4

品物 5

3
4
1
2
33
32
10
28
11
8
10
14
ナップサックの重量制限：75459東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

104

ナップサック問題の 1/2-近似解法• 緩和問題はもとの問題の制約条件を緩めて得られた問題ですから，(ナップサック問題の最適値) ≦ (緩和問題の最適値)が成り立ちます．
• ナップサック問題に貪欲算法を適用したときの目的関数の値を 𝑓 g で表し，品物 1, . . . , 𝑛 の価格の最大値を 𝑐 max で表すと，(緩和問題の最適値) ≦ 𝑓 g + 𝑐max
が成り立ちます．
• したがって，max{ 𝑓 g , 𝑐max }max{ 𝑓 g , 𝑐max } 1
≧
≧
g
max
𝑓 +𝑐
2
(ナップサック問題の最適値)が成り立つことがわかります．つまり，「貪欲算法の解」と「価格が 𝑐 maxの品物のみを詰め込む解」のうち目的関数の値が大きい方を出力する解法を考えると，この解法で得られる解の目的関数の値は最適値の 1/2 倍以上であることが保証されます．
• このときの比 1/2 のことを，近似比とよびます．また，近似比が 𝛼 以上である解法を，𝛼-近似解法とよびます．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

105

巡回セールスマン問題とは• いくつかの点と，任意の 2 点の間の距離（正の実数）が与えられているとします．2 点間を結ぶ線分を辺とよび，距離をその辺の重みとよびます．
• ある点から出発し，それぞれの点をちょうど 1 回ずつ訪れてから，もとの点に戻る経路を考えます（このような経路を，巡回路とよびます）．
• 巡回路のうち重みの和が最小のものを求める問題を，巡回セールスマン問
題とよびます．
• 下の例では，辺の重みは 2 点間のユークリッド距離としています：

点が 5 個ある例題

最適解• 巡回セールスマン問題は，それぞれの辺を巡回路に用いるか用いないかという離散的な選択をする問題とみることができますので，離散最適化問題
です．
東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

106

巡回路の個数• 点が 𝑛 個の巡回セールスマン問題を考えます．
• 出発する点を 1 つ決めると，巡回路は「残りの 𝑛 − 1 個の点をどの順番に訪れるか」で決まります．したがって，巡回路は (𝑛 − 1)! 個だけあります．
• 例として 𝑛 = 100 の場合，巡回路は 99! ≃ 10156 個あります．

𝑛 = 100 の例題

最適解• 制約条件を満たす解の個数は有限個ですが，非常に多いので，それらを単純に列挙することで最適解をみつけることは事実上不可能です．
• 巡回セールスマン問題は，解くことが難しい最適化問題の代表例として知られています，東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA
107貪欲算法の例• 巡回セールスマン問題に対して，たとえば次のような貪欲算法が考えられます：• ある点 𝑠 を出発点とし，まず，𝑠 から最も近い点（対応する辺の重みが最小の点）に移動します．
• 「まだ訪れていない点のうち，今いる点から最も近い点（対応する辺の重みが最小の点）に移動する」ことを繰り返します．
• 貪欲算法では，必ずしも巡回セールスマン問題の最適解は得られません：s

貪欲算法により得られる解

最適解

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA108

最小木問題とその貪欲算法• ある点から出発していくつかの辺をたどって同じ点に戻ることができるとき，それらの辺の集合を閉路とよびます．
• 巡回セールスマン問題の問題設定で，閉路ができないように辺をいくつか選び，すべての点が辺でつながるようにします．点と選んだ辺との組を，
全域木とよびます．辺の重みの和が最小の全域木を，最小木とよびます．
• 最小木は，次のような貪欲算法で得られることが知られています（クラス
カルのアルゴリズムとよばれています）：• 重みが小さい順に，辺を採用します．ただし，その辺を選ぶことで閉路ができる場合には，その辺は採用しません．

𝑛 = 6 の例題

最小木

東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

109

巡回セールスマン問題の 1/2-近似解法 (1/2)• ここでは，辺の重みは三角不等式を満たしていることを仮定します．
• 下記は，巡回セールスマン問題の 1/2-近似解法です：• 最小木を求めます．
• その辺を 2 本ずつにし，相異なる向きを付けます．
• ある点 𝑠 を出発点とし，向き付きの辺を順にたどります．その際に，既に訪れた点はスキップします．

s

最小木

各辺を 2 本ずつに得られた巡回路東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

110

巡回セールスマン問題の 1/2-近似解法 (2/2)• ここでは，辺の重みは三角不等式を満たしていることを仮定します．
• 近似比が 1/2 であることは，次のようにして示せます：• 最適解から辺を 1 本とり除くと全域木が得られますから，
1(最適値) ≧ (最小木の重み) = (真ん中の図の重み) .
2• 辺の重みに関する三角不等式より，(真ん中の図の重み) ≧ (得られた巡回路の重み) .

s

最小木

各辺を 2 本ずつに得られた巡回路東京大学 数理・情報教育研究センター 寒野善博 2021 CC BY-NC-SA

111