2018-07-10

R で主成分分析 (1)

PCA, Principal Component Analysis(主成分分析)は、線形独立であるか自明でない多数の変数を、ばらつきが大きい順から線形独立な主成分 (PC, Principal Component ) と呼ばれる変数に変換(合成)します。データの次元を削減するために用いられます。

次元が削減されてなにが嬉しいかというと、それは、多数の変数=多次元の情報を、少ない次元で端的に表現できるということです。主成分はデータのばらつき=変動の大きい順に PC1, PC2, ... と変換されています。そのため、主成分のいくつかだけでデータのふるまい(特徴)を概ね説明できてしまう場合があります。

最近、PCA を駆使して何百、何千という変数を持つデータの解析をする業務が多くなりました。社内とは言え、他人が開発をしたシステムを使っているので、いざとなれば自力である程度検証できるようにしておかなければなりません。

ということで、簡単なデータを使って R で主成分分析をしてみます。

iris データの解析

R に組み込まれているデータセット iris を使って主成分分析をしてみます。sepal は花びらの「がく」、petal は「花弁」のことで、それぞれ長さと幅を測定したデータが並んでいます。また species は「種」で、アヤメ属 (iris) の花 setosa, versicolor, virginica の三種類あります。

> head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
> unique(iris$Species)
[1] setosa     versicolor virginica 
Levels: setosa versicolor virginica
> 

そこで、この三種類の花がどのような特徴を持っているか調べてみることにします。「種」を除けば変数が4つしかないので「がく」と「花弁」の長さと幅をボックスプロットで比較する方法と、同じ方法で主成分分析を使うとどのようになるかを比較してみます。

R スクリプト

本記事で使用したスクリプトを以下に示しました。このスクリプトでは ggplot2 と qcc の二つのパッケージを利用しています。お使いの環境にインストールされていなければ、install.packages('ggplot2') のようにしてインストールしてください。

library(ggplot2)
library(qcc)
# GLOBAL
theme_default <- theme(panel.background = element_blank(),
                       panel.border = element_rect(fill = NA),
                       panel.grid.major = element_blank(),
                       panel.grid.minor = element_blank(),
                       strip.background = element_blank(),
                       axis.text.x = element_text(colour = "black"),
                       axis.text.y = element_text(colour = "black"),
                       axis.ticks = element_line(colour = "black"),
                       plot.margin = unit(c(1, 1, 1, 1), "line"))
# -----------------------------------------------------------------------------
#  myboxplot
# -----------------------------------------------------------------------------
myboxplot <- function(pp, ymin, ymax) {
  pp <- pp + geom_boxplot()
  pp <- pp + scale_y_continuous(limits = c(ymin, ymax))
  pp <- pp + theme_default
  return(pp)
}
# -----------------------------------------------------------------------------
#  BOXPLOT for RAW DATA
# -----------------------------------------------------------------------------
df <- as.data.frame(iris)
#
y.min <- min(df[, 1:4])
y.max <- max(df[, 1:4])
# -----------------------------------------------------------------------------
# Sepal.Length
p <- ggplot(df, aes(x = Species, y = Sepal.Length))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Sepal.Width
p <- ggplot(df, aes(x = Species, y = Sepal.Width))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Petal.Length
p <- ggplot(df, aes(x = Species, y = Petal.Length))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Petal.Width
p <- ggplot(df, aes(x = Species, y = Petal.Width))
p <- myboxplot(p, y.min, y.max)
plot(p)
# -----------------------------------------------------------------------------
#  BOXPLOT for PCA DATA
# -----------------------------------------------------------------------------
df <- as.data.frame(iris)
row.names(df) <- paste(df$Species, row.names(df), sep="_") 
df$Species <- NULL
# PCA
df_pca <- prcomp(df)
df_out <- as.data.frame(df_pca$x)
df_out$Species <- sapply(strsplit(as.character(row.names(df)), "_"), "[[", 1)
#
y.min <- min(df_out[, 1:4])
y.max <- max(df_out[, 1:4])
# -----------------------------------------------------------------------------
# PC1
p <- ggplot(df_out, aes(x = Species, y = PC1))
p <- myboxplot(p, y.min, y.max)
plot(p)
# PC2
p <- ggplot(df_out, aes(x = Species, y = PC2))
p <- myboxplot(p, y.min, y.max)
plot(p)
# PC3
p <- ggplot(df_out, aes(x = Species, y = PC3))
p <- myboxplot(p, y.min, y.max)
plot(p)
# PC4
p <- ggplot(df_out, aes(x = Species, y = PC4))
p <- myboxplot(p, y.min, y.max)
plot(p)
# Pareto Chart
PCA <- df_pca$sdev
names(PCA) <- c("PC1", "PC2", "PC3", "PC4")
pareto <- pareto.chart(PCA, ylab = "Variation")
plot(pareto)

生データのボックスプロット

以下に、「がく」と「花弁」の長さと幅について、「種」別に描画したボックスプロットを示しました。議論が分かれるところですが、ここでは全体を比較するため縦軸のスケールを揃えてあります。

これらのボックスプロットをじ〜っと見較べて判る差異を以下にまとめました。

  • 「がく」と「花弁」の大きさは、概ね setosa, versicolor, virsinica の順に大きい。
  • 「花弁 (petal)」の長さの差異が、他の数値に比べて際立っている。
  • 「がく (sepal)」の幅の差異が、他の差異に比べて小さい。

主成分分析後のボックスプロット

「がく」と「花弁」の長さと幅について、prcomp で主成分分析をして、主成分を「種」別に描画したボックスプロットを示しました。生データのボックスプロットの時と同様に、全体を比較するため縦軸のスケールを揃えてあります。

主成分のばらつき(変動)の大きさをパレート図にしたものを下記に示しました。

パレート図の内訳をみると、主成分 PC1 と PC2 で全体の変動の 80% 以上を説明できている事がわかります。

> pareto
     
Pareto chart analysis for PCA
        Frequency   Cum.Freq.  Percentage Cum.Percent.
  PC1   2.0562689   2.0562689  68.9345126   68.9345126
  PC2   0.4926162   2.5488851  16.5145035   85.4490161
  PC3   0.2796596   2.8285447   9.3753300   94.8243460
  PC4   0.1543862   2.9829309   5.1756540  100.0000000
> 

参考サイト

  1. irisの正体 (R Advent Calendar 2012 6日目) - どんな鳥も
  2. R plot PCA using ggplot2 - Boqiang Hu

 

ブログランキング・にほんブログ村へ
にほんブログ村

2018-07-07

才能の枯渇はまだ、いや、あればだけど…

最近忙しいから疲れているのかな?独り言を書きたくなりました。

転職して職場が変わっても、その業界にいる限り、昔の同僚や先輩後輩に出くわすことがある。

「何だかんだ言ってもまだ半導体業界にいます。」と、かつてお世話になった先輩に某顧客の工場の喫煙室で挨拶をしたとき、「半導体業界以外に移る必要がなぜある?」と返され、返事ができなかった。

週の大半が出張の毎日。ホテルから出張先へ向う時間調整で、毎日ではないのだけれど NHK の朝の連続テレビ小説を観るようになってしまった。今放映している「半分、青い。」で、主人公が伸びない才能に諦めをつけて漫画家をやめる決意をする段階にさしかっている(第13週 仕事が欲しい 6月25日(月)~6月30日(土)放送)。この場面に心がチクチクする。いや主人公にではない、自分にである。

漫画家に限らず、プロの仕事であればなんでも「才能」が必要とされる。自分は転職を重ねたけど、半導体業界の片隅でエンジニアとしてまだなんとかやっている。でも50歳半ばになり、年寄り臭く過去を振り返ることが多くなってしまった。一体自分が一番輝いていたのはいつだったのだろうかと。いや、過去を振り返っても仕方がないことなのだが、最近では才能ある若い世代をとても羨しく思うようになっている自分に気付く。

若い頃は負けるものかと、自分も切磋琢磨を続けて、自分ができないことをする同僚に追いつくことを目指したものだ。今は聞き分けが良い年嵩の同僚に成り下がってしまったようだ。もちろん少しでもキャッチアップしようと、知らなかった技術を仕入れて悪あがきをしているのではあるが…。

前述の番組でナレーターは、「才能が枯渇する」というような表現で主人公の状況を言い表したが、それは創造力を必要とする分野の仕事では致命的な事なのだろう。若い時に自分の才能の限界に気がついてしまったら、漫画家のような職業だともう廃業するしかない。そんな漫画家の「作品」はもはや読者に受け入れられるわけがないからだ。

幸い、エンジニアという職業はそうでもない。それは技術の蓄積の上に成り立っている仕事であるからだ。もちろん創造力を必要とする部分はあるが全てではない。きらきらと光り輝く才能が無くともなんとかやっていける。だから諦めきれずに悪あがきができる。

大学を卒業後、半導体プロセスエンジニアとしてキャリアをスタートした自分は、担当のドライエッチングプロセスの実験に熱中し、夜遅くまで SEM 室でエッチングしたサンプルの断面を観察していたものだ。日常的にプロセスの改善に明け暮れる中で DOESPC を駆使するにつれて統計解析を扱う仕事が増え、昨今 FDCVM に取り組むに至っては、半端ない量の多変量変数を扱うために機械学習の知識が必要になってきた。同じ場所、仕事に留まるとテクノロジーのトレンドから取り残される。だから新しい知識を習得するための努力は続く。

人生百年と言えど、職業人生はとっくに折り返し地点を超えてゴールを意識する歳になってしまった。そろそろラストスパートだ。

参考サイト

  1. NHK連続テレビ小説『半分、青い。』

 

ブログランキング・にほんブログ村へ
にほんブログ村