主成分分析

基本的な考え方

(Press ? for help, n and p for next and previous slide)

村田 昇

講義概要

  • 第1日 : 主成分分析の考え方
  • 第2日 : 分析の評価と視覚化

主成分分析の例

県毎の生活環境の違いの分析

県名 年少人口比 老年人口比 婚姻率 離婚率 高校数／面積 交通事故 犯罪件数 食料費 住居費 貯蓄率
北海道 11.7 26.0 4.86 2.12 1.34 274.2 8.98 22.5 5.9 19.4
青森県 12.1 27.0 4.33 1.78 2.63 386.7 6.12 24.6 5.4 16.2
岩手県 12.4 27.9 4.32 1.52 2.19 261.6 4.83 24.7 7.1 19.1
宮城県 13.0 22.9 5.30 1.70 3.18 447.7 8.85 23.7 5.0 9.5
秋田県 11.1 30.7 3.78 1.41 1.85 266.2 4.12 22.9 7.9 15.2
山形県 12.6 28.3 4.24 1.46 2.24 614.9 5.54 22.2 6.0 5.1
福島県 12.9 26.1 4.73 1.64 2.65 498.9 8.13 22.7 5.7 25.1
茨城県 13.2 23.8 4.92 1.79 3.09 500.6 13.00 21.9 6.8 20.0
栃木県 13.2 23.2 5.13 1.85 2.68 404.3 11.53 20.4 6.9 17.4
群馬県 13.4 24.9 4.64 1.77 3.56 925.2 10.49 23.9 4.4 1.6
埼玉県 13.0 22.0 5.10 1.86 7.81 493.6 13.91 23.0 7.0 21.6
千葉県 12.8 23.2 5.19 1.86 5.24 370.2 13.36 26.2 5.2 8.9
東京都 11.3 21.3 6.75 1.91 31.03 358.5 14.13 25.1 7.7 18.0
神奈川県 13.0 21.5 5.68 1.85 16.09 408.6 9.46 24.5 7.4 15.9
新潟県 12.5 27.2 4.35 1.37 2.35 357.2 8.71 23.5 5.7 13.0

ja-pairs1.png

Figure 1: 県別の生活環境(教育・労働などに関連する項目)

ja-pairs2.png

Figure 2: 県別の生活環境(教育・労働などに関連する項目)

ja-pairs3.png

Figure 3: 県別の生活環境(教育・労働などに関連する項目)

ja-biplot.png

Figure 4: 県別の生活環境の主成分分析

主成分分析の考え方

主成分分析

  • 多数の変量のもつ情報の分析・視覚化
    • 変量を効率的に縮約して少数の特徴量を構成する
    • 特徴量に関与する変量間の関係を明らかにする
  • PCA (Principal Component Analysis)
    • 構成する特徴量 : 主成分 (princial component)

分析の枠組み

  • \(x_{1},\dotsc,x_{p}\) : 変数
  • \(z_{1},\dotsc,z_{d}\) : 特徴量 ( \(d\leq p\) )

  • 変数と特徴量の関係 (線形結合)

    \begin{equation} z_k=a_{1k}x_{1}+\cdots+a_{pk}x_{p}\quad(k=1,\dotsc,d) \end{equation}

    • 特徴量は定数倍の任意性があるので以下を仮定

      \begin{equation} \|\boldsymbol{a}_k\|^2=\sum_{j=1}^pa_{jk}^2=1 \end{equation}

主成分分析の用語

  • 特徴量 \(z_k\)
    • 第 \(k\) 主成分得点 (principal component score)
    • 第 \(k\) 主成分
  • 係数ベクトル \(\boldsymbol{a}_k\)
    • 第 \(k\) 主成分負荷量 (principal component loading)
    • 第 \(k\) 主成分方向 (principal component direction)

分析の目的

  • 目的

    主成分得点 \(z_{1},\dots,z_{d}\) が変数 \(x_{1},\dotsc,x_{p}\) の情報を効率よく反映するように主成分負荷量 \(\boldsymbol{a}_{1},\dotsc,\boldsymbol{a}_{d}\) を観測データから決定する

  • 分析の方針 (以下は同値)
    • データの情報を最も保持する変量の 線形結合を構成
    • データの情報を最も反映する 座標軸を探索
  • 教師なし学習 の代表的手法の1つ
    • 特徴抽出 : 情報処理に重要な特性を変数に凝集
    • 次元縮約 : 入力をできるだけ少ない変数で表現

実習

R : 主成分分析を実行する関数

  • Rの標準的な関数
    • stats::prcomp()
    • stats::princomp()
  • 計算法に若干の違いがある
    • 数値計算の観点からみると prcomp() が優位
    • princomp() はS言語(商用)との互換性を重視した実装
  • 本講義では prcomp() を利用

R : 関数 prcomp() の使い方

  • データフレームの全ての列を用いる場合 

    prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE,
       tol = NULL, rank. = NULL, ...)
#' x: 必要な変数のみからなるデータフレーム
#' center: 中心化(平均0)を行って処理するか否か
#' scale.: 規格化(分散1)を行って処理するか否か

  • 列名を指定する(formulaを用いる)場合 

    prcomp(formula, data = NULL, subset, na.action, ...)
#' formula: ~ 変数名 (解析の対象を + で並べる) 左辺はないので注意
#' data: 必要な変数を含むデータフレーム
#' 詳細は '?stats::prcomp' を参照

R : 分析結果を取得する関数

  • 主成分得点を計算する 

    predict(object, newdata, ...)
#' object: prcomp が出力したオブジェクト
#' newdata: 主成分得点を計算するデータフレーム
#' 詳細は '?stats::prcomp' または '?stats::predict.prcomp' を参照
    • newdata を省略すると分析に用いたデータフレームの得点が計算される

  • 主成分分析の結果を取得する 

    tidy(x, matrix = "u", ...)
#' x: prcomp が出力したオブジェクト
#' matrix: 結果として取り出す行列 u:scores, v:loadings, d:eigenvalues
#' 詳細は '?broom::tidy.prcomp' を参照

  • 主成分得点を計算する 

    augment(x, data = NULL, newdata, ...)
#' x: prcomp が出力したオブジェクト
#' data: 元のデータ (通常不要)
#' newdata: 主成分得点を計算するデータフレーム
#' 詳細は '?broom::augment.prcomp' を参照

練習問題

  • 数値実験により主成分分析の考え方を確認しなさい

    • 以下のモデルに従う人工データを生成する 

      #' 観測データ (2次元) の作成 (aのスカラー倍に正規乱数を重畳)
a <- c(1, 2)/sqrt(5) # 主成分負荷量 (単位ベクトル)
n <- 100 # データ数
toy_data <- # aのスカラー倍に正規乱数を重畳(行列として作成)
  runif(n,-1,1) %o% a + rnorm(2*n, sd=0.3) 
colnames(toy_data) <- paste0("x", 1:2) # 行列に列名(x1,x2)を付与
toy_data <- as_tibble(toy_data)
    • 観測データの散布図を作成

    • 観測データから第1主成分負荷量を推定 

      prcomp(toy_data) # 全ての主成分を計算する
a_hat <- prcomp(toy_data)$rotation[,1] # 負荷量(rotation)の1列目が第1主成分

    • 散布図上に主成分負荷量を描画 

      geom_abline(slope = 傾き, intercept = 切片) # 指定の直線を追加できる

第1主成分の計算

記号の準備

  • 変数 : \(x_{1},\dotsc,x_{p}\) (\(p\)次元)

  • 観測データ : \(n\) 個の \((x_{1},\dotsc,x_{p})\) の組

    \begin{equation} \{(x_{i1},\dots,x_{ip})\}_{i=1}^n \end{equation}
  • ベクトル表現
    • \(\boldsymbol{x}_{i}=(x_{i1},\dots,x_{ip})^{\mathsf{T}}\) : \(i\) 番目の観測データ (\(p\) 次元空間内の1点)
    • \(\boldsymbol{a}=(a_{1},\dots,a_{p})^{\mathsf{T}}\) : 長さ1の \(p\) 次元ベクトル

係数ベクトルによる射影

  • データ \(\boldsymbol{x}_{i}\) の \(\boldsymbol{a}\) 方向成分の長さ

    \begin{equation} \boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_{i} \quad\text{(スカラー)} \end{equation}

  • 方向ベクトル \(\boldsymbol{a}\) をもつ直線上への点 \(\boldsymbol{x}_{i}\) の直交射影

    \begin{equation} (\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_{i})\,\boldsymbol{a} \quad\text{(スカラー \(\times\) ベクトル)} \end{equation}

幾何学的描像

pca-figure.png

Figure 5: 観測データの直交射影 (\(p=2,n=2\) の場合)

ベクトル \(\boldsymbol{a}\) の選択の指針

  • 射影による特徴量の構成

    ベクトル \(\boldsymbol{a}\) を うまく 選んで 観測データ \(\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{n}\) の情報を最も保持する1変量データ \(z_{1},\cdots,z_{n}\)を構成

    \begin{equation} z_{1}=\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_{1}, z_{2}=\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_2, \dotsc, z_{n}=\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_n \end{equation}

  • 特徴量のばらつきの最大化

    観測データの ばらつき を最も反映するベクトル \(\boldsymbol{a}\) を選択

    \begin{equation} \arg\max_{\boldsymbol{a}} \sum_{i=1}^n(\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_{i} -\boldsymbol{a}^{\mathsf{T}}\bar{\boldsymbol{x}})^2, \quad \bar{\boldsymbol{x}} = \frac{1}{n}\sum_{i=1}^n\boldsymbol{x}_{i}, \end{equation}

ベクトル \(\boldsymbol{a}\) の最適化

  • 最適化問題

    制約条件 \(\|\boldsymbol{a}\|=1\) の下で 以下の関数を最大化せよ

    \begin{equation} f(\boldsymbol{a}) = \sum_{i=1}^n(\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_{i} -\boldsymbol{a}^{\mathsf{T}}\bar{\boldsymbol{x}})^2 \end{equation}
  • この最大化問題は必ず解をもつ
    • \(f(\boldsymbol{a})\) は連続関数
    • 集合 \(\{\boldsymbol{a}\in\mathbb{R}^p:\|\boldsymbol{a}\|=1\}\) はコンパクト(有界閉集合)

第1主成分の解

行列による表現

  • 中心化したデータ行列

    \begin{equation} X = \begin{pmatrix} \boldsymbol{x}_{1}^{\mathsf{T}}-\bar{\boldsymbol{x}}^{\mathsf{T}} \\ \vdots \\ \boldsymbol{x}_{n}^{\mathsf{T}}-\bar{\boldsymbol{x}}^{\mathsf{T}} \end{pmatrix} = \begin{pmatrix} x_{11}-\bar{x}_{1} & \cdots & x_{1p}-\bar{x}_{p}\\ \vdots & & \vdots \\ x_{n1}-\bar{x}_{1} & \cdots & x_{np}-\bar{x}_{p} \end{pmatrix} \end{equation}

  • 評価関数 \(f(\boldsymbol{a})\) は行列 \(X^{\mathsf{T}}X\) の二次形式

    \begin{equation} f(\boldsymbol{a}) = \boldsymbol{a}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{a} \end{equation}

ベクトル \(\boldsymbol{a}\) の解

  • 最適化問題

    \begin{equation} \text{maximize}\quad f(\boldsymbol{a}) = \boldsymbol{a}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{a} \quad\text{s.t.}\quad \boldsymbol{a}^{\mathsf{T}}\boldsymbol{a}=1 \end{equation}

  • 制約付き最適化なので未定係数法を用いればよい

    \begin{equation} L(\boldsymbol{a},\lambda) =f(\boldsymbol{a})+\lambda(1-\boldsymbol{a}^{\mathsf{T}}\boldsymbol{a}) \end{equation}

    の鞍点

    \begin{equation} \frac{\partial}{\partial\boldsymbol{a}}L(\boldsymbol{a},\lambda) =0 \end{equation}

    を求めればよいので

    \begin{align} 2X^{\mathsf{T}}X\boldsymbol{a}-2\lambda\boldsymbol{a} &=0\\ X^{\mathsf{T}}X\boldsymbol{a} &=\lambda\boldsymbol{a} \quad\text{(固有値問題)} \end{align}

  • 解の条件

    \(f(\boldsymbol{a})\) の極大値を与える \(\boldsymbol{a}\) は \(X^{\mathsf{T}}X\) の固有ベクトルとなる

    \begin{equation} X^{\mathsf{T}}X\boldsymbol{a} = \lambda\boldsymbol{a} \end{equation}

第1主成分

  • 固有ベクトル\(\boldsymbol{a}\)に対する\(f(\boldsymbol{a})\) は行列 \(X^{\mathsf{T}}X\) の固有値

    \begin{equation} f(\boldsymbol{a}) =\boldsymbol{a}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{a} =\boldsymbol{a}^{\mathsf{T}}\lambda\boldsymbol{a} =\lambda \end{equation}
  • 求める \(\boldsymbol{a}\) は行列 \(X^{\mathsf{T}}X\) の最大固有ベクトル (長さ1)
  • 第1主成分負荷量 : 最大(第一)固有ベクトル \(\boldsymbol{a}\)

  • 第1主成分得点

    \begin{equation} z_{i1} =a_{1}x_{i1}+\cdots+a_{p}x_{ip} =\boldsymbol{a}^{\mathsf{T}}\boldsymbol{x}_{i}, \quad(i=1,\dots,n) \end{equation}

実習

練習問題

  • 第1主成分とGram行列の固有ベクトルの関係を調べなさい
    • 人工データを生成する
    • 主成分分析を実行する

    • Gram 行列を計算し固有値・固有ベクトルを求める 

      #' 中心化を行う
X <- scale(toy_data, scale = FALSE)
#' 詳細は '?base::scale' を参照
#' Gram 行列を計算する
G <- crossprod(X)
#' 固有値・固有ベクトルを求める
eigen(G) # 返り値 'values, vectors' を確認
#' 詳細は '?base::eigen' を参照

Gram 行列の性質

Gram 行列の固有値

  • \(X^{\mathsf{T}}X\) は非負定値対称行列
  • \(X^{\mathsf{T}}X\) の固有値は0以上の実数

    • 固有値を重複を許して降順に並べる

      \begin{equation} \lambda_{1}\geq\dotsb\geq\lambda_{p}\quad(\geq0) \end{equation}

    • 固有値 \(\lambda_{k}\) に対する固有ベクトルを \(\boldsymbol{a}_{k}\)(長さ1)とする

      \begin{equation} \|\boldsymbol{a}_{k}\|=1, \quad (k=1,\dotsc,p) \end{equation}

Gram 行列のスペクトル分解

  • \(\boldsymbol{a}_{1},\dotsc,\boldsymbol{a}_{p}\) は 互いに直交 するようとることができる

    \begin{equation} j\neq k \quad\Rightarrow\quad \boldsymbol{a}_{j}^{\mathsf{T}}\boldsymbol{a}_k=0 \end{equation}

  • 行列 \(X^{\mathsf{T}}X\) (非負定値対称行列) のスペクトル分解

    \begin{align} X^{\mathsf{T}}X &=\lambda_{1}\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}}+ \lambda_{2}\boldsymbol{a}_{2}\boldsymbol{a}_{2}^{\mathsf{T}}+ \dotsb+\lambda_{p}\boldsymbol{a}_{p}\boldsymbol{a}_{p}^{\mathsf{T}}\\ &=\sum_{k=1}^{p}\lambda_{k}\boldsymbol{a}_{k}\boldsymbol{a}_{k}^{\mathsf{T}} \end{align}
    • 固有値と固有ベクトルによる行列の表現

第2主成分以降の計算

第2主成分の考え方

  • 第1主成分
    • 主成分負荷量 : ベクトル \(\boldsymbol{a}_{1}\)
    • 主成分得点 : \(\boldsymbol{a}_{1}^{\mathsf{T}}\boldsymbol{x}_{i}\) (\(i=1,\dotsc,n\))

  • 第1主成分負荷量に関してデータが有する情報

    \begin{equation} (\boldsymbol{a}_{1}^{\mathsf{T}}\boldsymbol{x}_{i})\,\boldsymbol{a}_{1} \quad(i=1,\dotsc,n) \end{equation}

  • 第1主成分を取り除いた観測データ (分析対象)

    \begin{equation} \tilde{\boldsymbol{x}}_{i} = \boldsymbol{x}_{i} -(\boldsymbol{a}_{1}^{\mathsf{T}}\boldsymbol{x}_{i})\,\boldsymbol{a}_{1} \quad(i=1,\dotsc,n) \end{equation}

第2主成分の最適化

  • 最適化問題

    制約条件 \(\|\boldsymbol{a}\|=1\) の下で 以下の関数を最大化せよ

    \begin{equation} \tilde{f}(\boldsymbol{a}) = \sum_{i=1}^n(\boldsymbol{a}^{\mathsf{T}}\tilde{\boldsymbol{x}}_{i} -\boldsymbol{a}^{\mathsf{T}}\bar{\tilde{\boldsymbol{x}}})^2 \quad\text{ただし}\quad \bar{\tilde{\boldsymbol{x}}} = \frac{1}{n}\sum_{i=1}^n\tilde{\boldsymbol{x}}_{i} \end{equation}

第2主成分以降の解

行列による表現

  • 中心化したデータ行列

    \begin{equation} \tilde{X} = \begin{pmatrix} \tilde{\boldsymbol{x}}_{1}^{\mathsf{T}}-\bar{\tilde{\boldsymbol{x}}}^{\mathsf{T}} \\ \vdots \\ \tilde{\boldsymbol{x}}_{n}^{\mathsf{T}}-\bar{\tilde{\boldsymbol{x}}}^{\mathsf{T}} \end{pmatrix} = X-X\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}} \end{equation}

  • Gram 行列

    \begin{align} \tilde{X}^{\mathsf{T}}\tilde{X} &= (X-X\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}})^{\mathsf{T}} (X-X\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}})\\ &= X^{\mathsf{T}}X - X^{\mathsf{T}}X\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}} - \boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}}X^{\mathsf{T}}X + \boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}}X^{\mathsf{T}}X \boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}}\\ &= X^{\mathsf{T}}X - \lambda_{1}\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}} - \lambda_{1}\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}} + \lambda_{1}\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}} \boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}}\\ &= X^{\mathsf{T}}X - \lambda_{1}\boldsymbol{a}_{1}\boldsymbol{a}_{1}^{\mathsf{T}}\\ &= \sum_{k=2}^{p}\lambda_{k}\boldsymbol{a}_{k}\boldsymbol{a}_{k}^{\mathsf{T}} \end{align}

第2主成分

  • Gram 行列 \(\tilde{X}^{\mathsf{T}}\tilde{X}\) の固有ベクトル \(\boldsymbol{a}_{1}\) の固有値は 0

    \begin{equation} \tilde{X}^{\mathsf{T}}\tilde{X}\boldsymbol{a}_{1} = 0 \end{equation}
  • Gram 行列 \(\tilde{X}^{\mathsf{T}}\tilde{X}\) の最大固有値は \(\lambda_2\)
  • 解は第2固有値 \(\lambda_2\) に対応する固有ベクトル \(\boldsymbol{a}_2\)
  • 以下同様に 第 \(k\) 主成分負荷量は \(X^{\mathsf{T}}X\) の第 \(k\) 固有値 \(\lambda_k\) に対応する固有ベクトル \(\boldsymbol{a}_k\)

実習

データセットの準備

  • 主成分分析では以下のデータセットを使用する

    • japan_social.csv (配付)

      総務省統計局より取得した都道府県別の社会生活統計指標の一部

      • Pref : 都道府県名
      • Forest : 森林面積割合 (%) 2014年
      • Agri : 就業者１人当たり農業産出額(販売農家）(万円) 2014年
      • Ratio : 全国総人口に占める人口割合 (%) 2015年
      • Land : 土地生産性（耕地面積１ヘクタール当たり）(万円) 2014年
      • Goods : 商業年間商品販売額［卸売業＋小売業］（事業所当たり）(百万円) 2013年
      • Area : 地方区分

練習問題

  • 前掲のデータを用いて主成分分析を行いなさい

    • 都道府県名を行名としてデータを読み込む 

      js_data <- read_csv("data/japan_social.csv")
    • データの散布図行列を描く
    • 各データの箱ひげ図を描き，変数の大きさを確認する

    • 主成分負荷量を計算する 

      js_pca <- prcomp(js_data[-c(1,7)], scale. = TRUE)
#' '-c(1,7)' は都道府県名・地方区分を除く．関数 select() を利用することもできる
#' 'scale.=TRUE' とすると変数を正規化してから解析する

次回の予定

  • 第1日 : 主成分分析の考え方
  • 第2日 : 分析の評価と視覚化