Processing math: 100%

R演習 : 主成分分析

潜在的意味分析

テキストデータからBoW形式データへ

library(tm) #パッケージtmをインストールの上読み込む

データを読み込んで中身を確認する．これは単なるテキストデータである．

rawtext <- read.csv("test_sentences.csv",stringsAsFactors=F,header=F)[,1]
rawtext

1. 'Luke I am your father.'
2. 'When I was your age we didn\'t have Star Destroyers...'
3. 'Luke may I have some privacy please'
4. 'No Luke you\'re not old enough'
5. 'Hey Dad look at me! I\'m jedi!'
6. 'Harry Potter was a highly unusual boy in many ways.'
7. 'Get out of the way Harry Potter.'
8. 'Listen to me Harry.'

テキストデータからコーパスを作成する．

sentences <- VCorpus(VectorSource(rawtext))

各文章へのアクセスにはcontentを使う．

sentences[[1]]$content

'Luke I am your father.'

意味のない単語を削除する．大文字を小文字へ，数字を削除，ストップワードを削除など．

sent.clean <- tm_map(sentences, content_transformer(tolower))
sent.clean <- tm_map(sent.clean, removeNumbers)
sent.clean <- tm_map(sent.clean, removeWords, stopwords("en"))
sent.clean <- tm_map(sent.clean, removePunctuation)
sent.clean <- tm_map(sent.clean, stemDocument)
sent.clean <- tm_map(sent.clean, stripWhitespace)

削除後の文章は以下のようになる．

sent.clean[[1]]$content; sentences[[1]]$content

'luke father'

'Luke I am your father.'

BoW形式へのデータへ変換する．wordLengthsを指定することで単語の長さを制限できる．下記の場合は単語長が2~20のものに限定している．

bow <- DocumentTermMatrix(sent.clean,control=list(wordLengths=c(2,20)))

中身を確認すると，各行(文章)における出現単語(列)の有無が0 or 1で与えられている．

as.matrix(bow)

A matrix: 8 × 23 of type dbl

	age	boy	dad	destroy	enough	father	get	harri	hey	high	⋯	luke	mani	may	old	pleas	potter	privaci	star	unusu	way
1	0	0	0	0	0	1	0	0	0	0	⋯	1	0	0	0	0	0	0	0	0	0
2	1	0	0	1	0	0	0	0	0	0	⋯	0	0	0	0	0	0	0	1	0	0
3	0	0	0	0	0	0	0	0	0	0	⋯	1	0	1	0	1	0	1	0	0	0
4	0	0	0	0	1	0	0	0	0	0	⋯	1	0	0	1	0	0	0	0	0	0
5	0	0	1	0	0	0	0	0	1	0	⋯	0	0	0	0	0	0	0	0	0	0
6	0	1	0	0	0	0	0	1	0	1	⋯	0	1	0	0	0	1	0	0	1	1
7	0	0	0	0	0	0	1	1	0	0	⋯	0	0	0	0	0	1	0	0	0	1
8	0	0	0	0	0	0	0	1	0	0	⋯	0	0	0	0	0	0	0	0	0	0

2つ以上の文で出現する単語のみ取得することもoptionで可能．

rbow <- bow[,findFreqTerms(bow,2)]

as.matrix(rbow)

A matrix: 8 × 4 of type dbl

	harri	luke	potter	way
1	0	1	0	0
2	0	0	0	0
3	0	1	0	0
4	0	1	0	0
5	0	0	0	0
6	1	0	1	1
7	1	0	1	1
8	1	0	0	0

特異値分解(SVD)による行列の近似

任意の実$n\times p$行列$X$は$X=UDV^{T}$のように分解が可能である．ここで，$D={\rm diag}(d_{1},\dots, d_{r}), d_{j} > 0, r={\rm rank}(X)$であり， $U^{T}U=V^{T}V=I$である．

特異値分解をR上で実行する．

result <- svd(bow)
str(result) #中身を確認

List of 3
 $ d: num [1:8] 3.02 2.28 2 1.73 1.57 ...
 $ u: num [1:8, 1:8] 3.95e-18 -8.58e-17 3.33e-17 -4.12e-17 -5.72e-16 ...
 $ v: num [1:23, 1:8] 1.72e-16 -2.75e-01 -4.29e-16 1.40e-16 2.80e-18 ...

例えば特異値$d_{j}$を取得するためには以下のように行えばよい．

result$d

1. 3.02005186652227
2. 2.28348850300971
3. 2
4. 1.73205080756888
5. 1.56869727072788
6. 1.52289849543483
7. 1.24902637927215
8. 1.15102959537683

対角行列$D$のサイズを$r$から$k (<r)$に変更することで，元の行列$X$を近似することができる．

U <- as.matrix(result$u)
D <- diag(result$d)
V <- as.matrix(result$v)

ap.bow <- U[,1:2] %*% D[1:2,1:2] %*% t(V[,1:2])  #k=2とした近似
sum(abs(bow-ap.bow))  ##近似の精度

23.7238335706983

ap.bow <- U[,1:5] %*% D[1:5,1:5] %*% t(V[,1:5])  #k=5とした近似
sum(abs(bow-ap.bow))  ##近似の精度

12.4068964307878

$k=r$ならばもちろん等しい．

ap.bow <- U[,1:8] %*% D[1:8,1:8] %*% t(V[,1:8])  #k=rとした近似
sum(abs(bow-ap.bow))  ##近似の精度

4.0626071112701e-14

行列$V$の中身を解釈することで各文章がどのトピックに属しているかを解釈することができる．

rownames(V) <- colnames(bow)  #名前を付ける

round(V[,1],3)  ##harriとpotterの絶対値が大きい
round(V[,2],3)  ##lukeの絶対値が大きい

age

0

boy

-0.275

dad

0

destroy

0

enough

0

father

0

get

-0.173

harri

-0.511

hey

0

high

-0.275

jedi

0

listen

-0.063

look

0

luke

0

mani

-0.275

may

0

old

0

pleas

0

potter

-0.448

privaci

0

star

0

unusu

-0.275

way

-0.448

age

0

boy

0

dad

0

destroy

0

enough

-0.228

father

-0.174

get

0

harri

0

hey

0

high

0

jedi

0

listen

0

look

0

luke

-0.733

mani

0

may

-0.331

old

-0.228

pleas

-0.331

potter

0

privaci

-0.331

star

0

unusu

0

way

0

おそらく2つのトピックが存在しているので，$UD$の第1列と第2列で散布図を描く．

plot(D[1,1]*U[,1], D[2,2]*U[,2])

主成分分析

データサイエンスコンソーシアムの民力データ(http://datascience.jp/tutorial.html) を読み込む．

Ppower <- read.csv("Ppower.csv",header=T,row.names=1)
head(Ppower)

A data.frame: 6 × 24

	population	households	office	income	nationaltax	localtax	agriculturalout	forestryprod	landing	factory	⋯	deposit	publicexpense	houses	cars	educationalamount	booksales	newspapersales	tv	phone	posting
	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<dbl>	<int>	<int>	⋯	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<int>
Hokkaido	5650573	2522295	270504	145293	1379098	546638	10579	463.7	761331	10668	⋯	336487	1611254	49183	3720320	945495	177154	2235101	1652081	6154045	650837
Aomori	1479358	551806	74341	32498	242730	131442	2402	88.2	107727	3188	⋯	70174	324170	8971	999912	264053	40213	531167	470546	1384375	114626
Iwate	1405060	488354	72456	34152	224269	123476	2587	193.4	108752	4126	⋯	76438	339411	8906	983234	277808	35185	471524	443942	1300483	114618
Miyagi	2350026	856527	115297	61092	735121	245372	1870	76.1	373037	5782	⋯	125693	391923	19382	1557230	365365	70995	812314	719845	2525483	351549
Akita	1173722	410308	65300	27294	186373	97512	2208	123.2	5599	4130	⋯	56790	328282	6664	831589	214608	30773	442938	385189	1087529	93834
Yamagata	1225990	387732	70523	29851	202270	107894	2349	62.5	4116	5971	⋯	67083	363756	7088	920699	215490	29966	462638	371689	1127332	105540

24変数を合成して新しい変数を作っていく．

xx <- scale(Ppower,scale=T) #scaling

PCAの実行.

sigma <- t(xx)%*%xx/47
result <- eigen(sigma)
evalue <- result$values
evector <- result$vectors

rownames(evector) <- colnames(xx) #解釈をしやすくするために名前をつける．

round(evector[,1],3)
round(evector[,2],3)

population

-0.225

households

-0.227

office

-0.228

income

-0.229

nationaltax

-0.207

localtax

-0.212

agriculturalout

-0.023

forestryprod

0.017

landing

-0.043

factory

-0.207

industryamount

-0.156

employer

-0.227

storesales

-0.205

powerconsum

-0.228

deposit

-0.223

publicexpense

-0.211

houses

-0.227

cars

-0.208

educationalamount

-0.223

booksales

-0.228

newspapersales

-0.228

tv

-0.226

phone

-0.229

posting

-0.209

population

0.053

households

0.03

office

0.002

income

-0.02

nationaltax

-0.163

localtax

-0.124

agriculturalout

0.582

forestryprod

0.445

landing

0.504

factory

-0.05

industryamount

0.083

employer

0.039

storesales

-0.149

powerconsum

0.021

deposit

-0.106

publicexpense

0.192

houses

-0.044

cars

0.186

educationalamount

0.119

booksales

-0.046

newspapersales

-0.002

tv

0.046

phone

-0.027

posting

-0.161

第2主成分スコアまで計算してプロットする．

z1 <- xx %*% evector[,1]
z2 <- xx %*% evector[,2]
plot(z1,z2,col="blue",pch=3)
text(z1,z2,row.names(xx),cex=0.9) #点にサンプル名を付与する．

主成分スコアの性質 : 主成分スコアは平均0で互いに相関しない．

Z <- xx %*% evector[,1:5] #第5主成分スコアまで例えば計算する．
round(apply(Z,2,mean),3) #各主成分スコアの平均値
round(t(Z)%*%Z/47,3) #主成分スコア間の共分散

1. 0
2. 0
3. 0
4. 0
5. 0

A matrix: 5 × 5 of type dbl

18.531	0.000	0.000	0.00	0.000
0.000	2.512	0.000	0.00	0.000
0.000	0.000	1.164	0.00	0.000
0.000	0.000	0.000	0.63	0.000
0.000	0.000	0.000	0.00	0.266

主成分スコアの分散，すなわち上記行列の対角成分は固有値である．

round(evalue,3)[1:5]

1. 18.531
2. 2.512
3. 1.164
4. 0.63
5. 0.266

biplotによる主成分スコアの解釈

biplot(Z[,1:2],evector[,1:2],ylim=c(-5,9),xlim=c(-25,7),cex=0.85)

主成分数の選択については，寄与率がひとつの尺度となる．累積寄与率が一定の値以上になるように主成分数を選択する．

round(cumsum(evalue)/sum(evalue),3)  #累積寄与率
plot(1:24,round(cumsum(evalue)/sum(evalue),3))  #そのプロット

1. 0.789
2. 0.896
3. 0.945
4. 0.972
5. 0.984
6. 0.991
7. 0.996
8. 0.998
9. 0.998
10. 0.999
11. 0.999
12. 0.999
13. 1
14. 1
15. 1
16. 1
17. 1
18. 1
19. 1
20. 1
21. 1
22. 1
23. 1
24. 1