とらりもん - 主成分分析 Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
2012/02/08 田中 健太郎
!1. はじめに
[[前回|教師なし分類による土地被覆分類図の作成]], 教師なし分類による土地被覆分類について学んだ. 教師なし分類に用いる衛星データとして, LandsatのBand1〜4のデータを使用しているが, 一般にマルチスペクトルのデータでは, 隣接するBand同士が相関を持っていると考えられる. つまりBand1〜4のデータを用いたデータセットは, 助長性を持つことになる. そのためマルチスペクトルデータを用いて, 土地被覆分類を行う場合, 主成分分析と呼ばれる手法がよく使われる. 本章では, 主成分分析について学ぶことにする.
!2. 主成分分析とは?
主成分分析は, 統計手法の一つであり, 土地被覆分類に限らず様々な分野で一般的に使用されている手法である. 簡単に説明すると, 「主成分分析とは, 相関関係を持つ観測値のデータから, 互いに無関係な因子を取り出し, 観測値をその因子の線型結合で表す手法である.」 またこの互いに無関係な因子のことを, 主成分と呼ぶ. 主成分分析の理解を助けるために, 簡単な演習問題を行おう.
!3. 主成分分析の演習
次の表は, とある学校の期末テストの結果である (もちろんすべてフィクションである). 生徒数は20人で, テスト科目は, 国語, 英語, 数学, 理科とする.
{{br}}
{{attach_view(pca1.png)}}
{{br}}
単純なケースとして, 生徒は理系向きか文系向きの生徒がいるとすると (もちろんオールラウンド型や体育会系, その他例外もいるのだが) , 国語と英語もしくは数学と理科の点数には, 相関関係があると考えられる. こうしたデータに対して主成分分析を行い, 主成分分析の意味と, 主成分がどのような意味を持つかを調べてみる.
※衛星データを用いるケースを考えた場合, 各教科が衛星データの各バンドデータやNDVI, 出席番号が位置情報に相当する.
観測データを20×4の行列Xとして考えると, 主成分分析の最初のステップは, 4×4の共分散行列Vを作成することである.
共分散行列Vは, 次の式で作成される.
{{br}}
{{attach_view(pca2.png)}}
{{br}}
次に共分散行列Vを対角化する. 行列Vの対角化とは, 行列Vの複数の固有ベクトルを列とする行列Uを用いて,
{{br}}
{{attach_view(pca3.png)}}
{{br}}
とすることである. 新たに定義した行列V'は, 行列Vの固有値を対角成分に持つ対角行列である. ここで値がもっとも大きい固有値に対応する固有ベクトルのことを, 第1主成分と呼び, n番目に大きい固有値に対応する固有ベクトルを第n主成分と呼ぶ. つまり主成分分析とは, 共分散行列の固有値問題であるといえる.固有ベクトルは, 定数倍しても固有ベクトルであることを思い出そう. つまり任意の定数を固有ベクトルに掛けることによって, 大きさ (ノルム) が1の単位固有ベクトルを作ることができる. 固有ベクトルは, 単位固有ベクトルに変換しておくと都合が良い.
実際に共分散行列Vから第一主成分の固有ベクトルおよび固有値を求めてみる.
手計算では困難なので, Rと呼ばれる無料の統計解析ソフトを使おう. Rは非常に優秀であり, 様々な統計解析をカバーしている. GRASSデータをRで処理することも可能である.
!!Rの入手とインストール
Ubuntuの場合
sudo apt-get update
sudo apt-get install r-base
sudo apt-get install r-base-dev
sudo apt-get install r-cran-*
その他OSの場合は, [[こちら|http://www.okada.jp.org/RWiki/?R%20%A4%CE%A5%A4%A5%F3%A5%B9%A5%C8%A1%BC%A5%EB]]を参考に,ダウンロードおよびインストールする.
インストールしたら
$ R
とするとRが起動する.
!!Rによる主成分分析
主成分分析のやり方はいたって簡単である.
まずテキストデータを適当な変数 (tmp) に代入する
tmp<-c(53,80,29,90,37,52,31,42,53,71,70,40,\
84,25,34,36,22,40,34,48,44,68,36,63,25,51,17,\
41,50,54,64,34,98,31,38,19,15,44,35,52,44,49,\
32,87,60,50,73,49,41,38,53,44,56,41,27,61,83,\
81,44,83,56,36,37,84,39,38,65,47,34,40,41,45,49,25,23,71,79,94,42,69 )
つぎにテキストデータを行列に変換し, 適当な変数 (mat) に代入する.
> mat<-matrix(tmp, 20,4, byrow=F)
行列の各列に名前付けを行う.
> colnames(mat)<-c("J", "E", "M", "S")
最後にprcompと呼ばれる関数によって, 共分散行列の作成, 行列の固有値, 固有ベクトルが求まる.結果を変数result.pcに出力する.
> result.pc<-prcomp(mat)
さっそく結果を表示してみよう.
> summary(result.pc)
PC1 PC2 PC3 PC4
Standard deviation 27.5571 26.0666 6.6951 6.39585
Proportion of Variance 0.4981 0.4457 0.0294 0.02683
Cumulative Proportion 0.4981 0.9438 0.9732 1.00000
PC1〜4は, 第1主成分〜第4主成分を意味している. Standard deviationは, 主成分方向における標準偏差であり, 2乗すると固有値になる. Proportion of Varianceは, 寄与率と呼ばれる値であり, 各主成分の標準偏差を標準偏差の合計 (27.5571+26.0666+6.6951+6.39585)で除算したものである. Cumulative Proportionは, 累積寄与率である.結果より第二主成分までで, 累積寄与率が94%となっている. 第3主成分以下は, 科目ごとの得点の大小による影響がほとんど無いことを意味している.そのためテストの点数は4科目分あるが, 実際には第1主成分の値と第2主成分の値の2つを考慮すれば, データの解析には十分であるといえる.
次に各主成分が何を意味しているかを見るために, 固有ベクトルの結果をみる.
> result.pc$rotation
PC1 PC2 PC3 PC4
J -0.70456635 0.141688348 0.69041127 0.08272215
E -0.70067461 0.009512981 -0.70131341 -0.13085907
M 0.04752647 0.660990134 -0.17441084 0.72829536
S 0.10189501 0.736835095 0.03275114 -0.66754686
結果より第1主成分は, 文系科目 (国語と英語) の点数が低い度合いを示していると解釈できる. つまり第一主成分の値が大きい場合, 文系科目の点数が低いことを意味している. 一方第二主成分は, 文系科目に関係なく, 理系科目 (数学と理科) の点数の出来を示す指標であると解釈できる. また主成分分析から推察される結果として, 今回のテストでは, 文系科目の点数が高い場合, 理系科目の点数が低い (もしくはその逆) といった相関はみられなかった.
Rを終了する場合,
> q()
と入力する.
以上がRを用いた主成分分析の方法である.続いて衛星データに対して, 主成分分析を行う方法について解説する.
!GRASSで行う主成分分析
GRASSでは, i.pcaというコマンドで主成分分析が実装されている.
[[前回|高分解能衛星Landsat: 植生指標NDVI, 温度画像]]と同じマップセットを使用し, i.pcaを使用してみる.
$ grass
LOCATION: UTM54
MAPSET: PERMANENT
$ g.region rast=B1
実行するのはたった一行である. (さすがGRASS!!)
$ i.pca in=B1L,B2L,B3L,B4L out=PCA rescale=0,0
実行すると,
$ g.list rast | cat | grep PCA
PCA.1
PCA.3PCA.2
PCA.3
PCA.4
というラスターデータが作成される. それぞれ第1主成分から第4主成分に対応している.
それぞれの寄与率を見るためには,r.infoを実行する.
r.info PCA.1
PC1 1712.60 (-0.5562,-0.4776,-0.4857,-0.4761) [81.82%]
PC2 313.03 ( 0.3810, 0.1851, 0.2383,-0.8739) [14.96%]
PC3 64.45 ( 0.6447,-0.0419,-0.7605, 0.0648) [ 3.08%]
PC4 3.04 ( 0.3603,-0.8578, 0.3590, 0.0733) [ 0.15%]
結果より, 第1と第2主成分で寄与率が96%を越えていることがわかる. また第1主成分は, デジタルナンバーの低さ度合いを示している. つまり第1主成分の値が大きい場合, どの観測バンドにおいてもデジタルナンバーが低いことを意味する. 第2主成分は, 主に観測バンド1〜3のデジタルナンバーと観測バンド4のデジタルナンバーとの差を意味している. つまり第2主成分が大きいほど, 観測バンド1〜3に比べ, 観測バンド4の値が低い.
主成分分析結果を用いて, 教師なし分類を行ってみる. 主成分分析の結果より, 第3主成分以下は寄与率が低いため, 解析に使用しない.
やり方がわからないという人は, [[前回|教師なし分類による土地被覆分類図の作成]]の内容を参照.
i.group group=landsat sub=pca input=PCA.1,PCA.2
i.cluster group=landsat subgroup=pca sigfile=sig_pca classes=10 report=class_report.txt
i.maxlik group=landsat subgroup=pca sigfile=sig_pca class=LC_map_pca reject=LC_map_pca_re
前回と同様に分類結果の分類項目に意味付けを行う.
cat << EOF > reclass_pca.txt
10 9 = 11 water
2 3 5 6 7 8 = 12 vegetation
4 = 13 urban area
1 = 14 bare
EOF
r.reclass input=LC_map_pca out=LC_map_pca_reclass rules=reclass_pca.txt title=LC_map_pca_reclass
d.rast.leg LC_map_pca_reclass
{{attach_view(pca4.png)}}
最後に色合いを調節して, 表示してみよう.
cat << EOF > color.txt
11 0:0:255
12 0:100:0
13 255:0:0
14 200:200:0
EOF
r.colors LC_map_pca_reclass color=rules < color.txt
g.region rast=LC_map_pca_reclass
d.rast.leg LC_map_pca_reclass
{{attach_view(pca5.png)}}
以上が主成分分析を用いた教師なし分類の方法である. 主成分分析を行うことで, 分類に用いるラスターデータが2つになっていることがわかる (特徴空間が2次元になる).
!1. はじめに
[[前回|教師なし分類による土地被覆分類図の作成]], 教師なし分類による土地被覆分類について学んだ. 教師なし分類に用いる衛星データとして, LandsatのBand1〜4のデータを使用しているが, 一般にマルチスペクトルのデータでは, 隣接するBand同士が相関を持っていると考えられる. つまりBand1〜4のデータを用いたデータセットは, 助長性を持つことになる. そのためマルチスペクトルデータを用いて, 土地被覆分類を行う場合, 主成分分析と呼ばれる手法がよく使われる. 本章では, 主成分分析について学ぶことにする.
!2. 主成分分析とは?
主成分分析は, 統計手法の一つであり, 土地被覆分類に限らず様々な分野で一般的に使用されている手法である. 簡単に説明すると, 「主成分分析とは, 相関関係を持つ観測値のデータから, 互いに無関係な因子を取り出し, 観測値をその因子の線型結合で表す手法である.」 またこの互いに無関係な因子のことを, 主成分と呼ぶ. 主成分分析の理解を助けるために, 簡単な演習問題を行おう.
!3. 主成分分析の演習
次の表は, とある学校の期末テストの結果である (もちろんすべてフィクションである). 生徒数は20人で, テスト科目は, 国語, 英語, 数学, 理科とする.
{{br}}
{{attach_view(pca1.png)}}
{{br}}
単純なケースとして, 生徒は理系向きか文系向きの生徒がいるとすると (もちろんオールラウンド型や体育会系, その他例外もいるのだが) , 国語と英語もしくは数学と理科の点数には, 相関関係があると考えられる. こうしたデータに対して主成分分析を行い, 主成分分析の意味と, 主成分がどのような意味を持つかを調べてみる.
※衛星データを用いるケースを考えた場合, 各教科が衛星データの各バンドデータやNDVI, 出席番号が位置情報に相当する.
観測データを20×4の行列Xとして考えると, 主成分分析の最初のステップは, 4×4の共分散行列Vを作成することである.
共分散行列Vは, 次の式で作成される.
{{br}}
{{attach_view(pca2.png)}}
{{br}}
次に共分散行列Vを対角化する. 行列Vの対角化とは, 行列Vの複数の固有ベクトルを列とする行列Uを用いて,
{{br}}
{{attach_view(pca3.png)}}
{{br}}
とすることである. 新たに定義した行列V'は, 行列Vの固有値を対角成分に持つ対角行列である. ここで値がもっとも大きい固有値に対応する固有ベクトルのことを, 第1主成分と呼び, n番目に大きい固有値に対応する固有ベクトルを第n主成分と呼ぶ. つまり主成分分析とは, 共分散行列の固有値問題であるといえる.固有ベクトルは, 定数倍しても固有ベクトルであることを思い出そう. つまり任意の定数を固有ベクトルに掛けることによって, 大きさ (ノルム) が1の単位固有ベクトルを作ることができる. 固有ベクトルは, 単位固有ベクトルに変換しておくと都合が良い.
実際に共分散行列Vから第一主成分の固有ベクトルおよび固有値を求めてみる.
手計算では困難なので, Rと呼ばれる無料の統計解析ソフトを使おう. Rは非常に優秀であり, 様々な統計解析をカバーしている. GRASSデータをRで処理することも可能である.
!!Rの入手とインストール
Ubuntuの場合
sudo apt-get update
sudo apt-get install r-base
sudo apt-get install r-base-dev
sudo apt-get install r-cran-*
その他OSの場合は, [[こちら|http://www.okada.jp.org/RWiki/?R%20%A4%CE%A5%A4%A5%F3%A5%B9%A5%C8%A1%BC%A5%EB]]を参考に,ダウンロードおよびインストールする.
インストールしたら
$ R
とするとRが起動する.
!!Rによる主成分分析
主成分分析のやり方はいたって簡単である.
まずテキストデータを適当な変数 (tmp) に代入する
tmp<-c(53,80,29,90,37,52,31,42,53,71,70,40,\
84,25,34,36,22,40,34,48,44,68,36,63,25,51,17,\
41,50,54,64,34,98,31,38,19,15,44,35,52,44,49,\
32,87,60,50,73,49,41,38,53,44,56,41,27,61,83,\
81,44,83,56,36,37,84,39,38,65,47,34,40,41,45,49,25,23,71,79,94,42,69 )
つぎにテキストデータを行列に変換し, 適当な変数 (mat) に代入する.
> mat<-matrix(tmp, 20,4, byrow=F)
行列の各列に名前付けを行う.
> colnames(mat)<-c("J", "E", "M", "S")
最後にprcompと呼ばれる関数によって, 共分散行列の作成, 行列の固有値, 固有ベクトルが求まる.結果を変数result.pcに出力する.
> result.pc<-prcomp(mat)
さっそく結果を表示してみよう.
> summary(result.pc)
PC1 PC2 PC3 PC4
Standard deviation 27.5571 26.0666 6.6951 6.39585
Proportion of Variance 0.4981 0.4457 0.0294 0.02683
Cumulative Proportion 0.4981 0.9438 0.9732 1.00000
PC1〜4は, 第1主成分〜第4主成分を意味している. Standard deviationは, 主成分方向における標準偏差であり, 2乗すると固有値になる. Proportion of Varianceは, 寄与率と呼ばれる値であり, 各主成分の標準偏差を標準偏差の合計 (27.5571+26.0666+6.6951+6.39585)で除算したものである. Cumulative Proportionは, 累積寄与率である.結果より第二主成分までで, 累積寄与率が94%となっている. 第3主成分以下は, 科目ごとの得点の大小による影響がほとんど無いことを意味している.そのためテストの点数は4科目分あるが, 実際には第1主成分の値と第2主成分の値の2つを考慮すれば, データの解析には十分であるといえる.
次に各主成分が何を意味しているかを見るために, 固有ベクトルの結果をみる.
> result.pc$rotation
PC1 PC2 PC3 PC4
J -0.70456635 0.141688348 0.69041127 0.08272215
E -0.70067461 0.009512981 -0.70131341 -0.13085907
M 0.04752647 0.660990134 -0.17441084 0.72829536
S 0.10189501 0.736835095 0.03275114 -0.66754686
結果より第1主成分は, 文系科目 (国語と英語) の点数が低い度合いを示していると解釈できる. つまり第一主成分の値が大きい場合, 文系科目の点数が低いことを意味している. 一方第二主成分は, 文系科目に関係なく, 理系科目 (数学と理科) の点数の出来を示す指標であると解釈できる. また主成分分析から推察される結果として, 今回のテストでは, 文系科目の点数が高い場合, 理系科目の点数が低い (もしくはその逆) といった相関はみられなかった.
Rを終了する場合,
> q()
と入力する.
以上がRを用いた主成分分析の方法である.続いて衛星データに対して, 主成分分析を行う方法について解説する.
!GRASSで行う主成分分析
GRASSでは, i.pcaというコマンドで主成分分析が実装されている.
[[前回|高分解能衛星Landsat: 植生指標NDVI, 温度画像]]と同じマップセットを使用し, i.pcaを使用してみる.
$ grass
LOCATION: UTM54
MAPSET: PERMANENT
$ g.region rast=B1
実行するのはたった一行である. (さすがGRASS!!)
$ i.pca in=B1L,B2L,B3L,B4L out=PCA rescale=0,0
実行すると,
$ g.list rast | cat | grep PCA
PCA.1
PCA.3
PCA.4
というラスターデータが作成される. それぞれ第1主成分から第4主成分に対応している.
それぞれの寄与率を見るためには,r.infoを実行する.
r.info PCA.1
PC1 1712.60 (-0.5562,-0.4776,-0.4857,-0.4761) [81.82%]
PC2 313.03 ( 0.3810, 0.1851, 0.2383,-0.8739) [14.96%]
PC3 64.45 ( 0.6447,-0.0419,-0.7605, 0.0648) [ 3.08%]
PC4 3.04 ( 0.3603,-0.8578, 0.3590, 0.0733) [ 0.15%]
結果より, 第1と第2主成分で寄与率が96%を越えていることがわかる. また第1主成分は, デジタルナンバーの低さ度合いを示している. つまり第1主成分の値が大きい場合, どの観測バンドにおいてもデジタルナンバーが低いことを意味する. 第2主成分は, 主に観測バンド1〜3のデジタルナンバーと観測バンド4のデジタルナンバーとの差を意味している. つまり第2主成分が大きいほど, 観測バンド1〜3に比べ, 観測バンド4の値が低い.
主成分分析結果を用いて, 教師なし分類を行ってみる. 主成分分析の結果より, 第3主成分以下は寄与率が低いため, 解析に使用しない.
やり方がわからないという人は, [[前回|教師なし分類による土地被覆分類図の作成]]の内容を参照.
i.group group=landsat sub=pca input=PCA.1,PCA.2
i.cluster group=landsat subgroup=pca sigfile=sig_pca classes=10 report=class_report.txt
i.maxlik group=landsat subgroup=pca sigfile=sig_pca class=LC_map_pca reject=LC_map_pca_re
前回と同様に分類結果の分類項目に意味付けを行う.
cat << EOF > reclass_pca.txt
10 9 = 11 water
2 3 5 6 7 8 = 12 vegetation
4 = 13 urban area
1 = 14 bare
EOF
r.reclass input=LC_map_pca out=LC_map_pca_reclass rules=reclass_pca.txt title=LC_map_pca_reclass
d.rast.leg LC_map_pca_reclass
{{attach_view(pca4.png)}}
最後に色合いを調節して, 表示してみよう.
cat << EOF > color.txt
11 0:0:255
12 0:100:0
13 255:0:0
14 200:200:0
EOF
r.colors LC_map_pca_reclass color=rules < color.txt
g.region rast=LC_map_pca_reclass
d.rast.leg LC_map_pca_reclass
{{attach_view(pca5.png)}}
以上が主成分分析を用いた教師なし分類の方法である. 主成分分析を行うことで, 分類に用いるラスターデータが2つになっていることがわかる (特徴空間が2次元になる).