ブログ一覧

ブログ始めました

テンプレートにblogが入っていたので、blogの項目を作りました。気が向いたら研究について何かしら書こうかと思います。

2017年08月29日ブログのカテゴリー:blog, その他

多元素データの生成

Rを使ってデータの可視化・解析を行うにあたって,自分たちが良く解析するデータを練習に使いたいと常々思っていました.まだ,論文にしていないものは公開するわけにもいかず… 必ずしも,解析したいようなデータがそろっていない…ということで,まずは実際の分析データに似たデータを生成するところから始めたいと思います.

ここでは,島根県内の水試料中に含まれる67元素の分析データの代表値から,多元素データを生成してみます.

一般的に,環境中の化学物質濃度は対数正規分布していることが多いので,幾何平均値と幾何標準偏差を計算したファイルを用意しました.さらに相関係数は,対数変換したのちにPearsonの相関係数を計算しました.

それぞれのファイルは以下の通りです (現在は非公開)

  • 幾何平均値:
  • 幾何標準偏差:
  • 相関係数行列:

正規分布と対数正規分布

library(mvtnorm)
# {MASS}の()では負の相関を扱えない

m = 67 # number of variable
g = 5
n = 60 # sample size
N = g *n

d.cor <- read.csv("cor.csv") 
m.cor <- as.matrix(d.cor[,-1])

GM <- read.csv("GM.csv")
GM <- GM[,-1]

GSD <- read.csv("GSD.csv")
GSD <- GSD[,-1]

#
df <- rmvnorm(n=N, mean=rep(0, m), sigma=m.cor)
# もう一度標準化しておく
df <- scale(
  df, 
  center = TRUE, # 平均値を0にする
  scale = TRUE # 標準偏差を1にする
  )

# for文を使ってセル[i, j]にlog(GSD[,j])を乗じて,log(GM[,j])を足して
# logから元に戻すために10^ の計算をする
for (i in 1:N)  {
  for (j in 1:m) {
    df[i, j] <- 10 ^ ( df[i, j]* log10(GSD[1, j]) + log10(GM[1, j])  )
  }
}

# 変数名を元素名に変更
colnames(df) <- colnames(GM)
d.ele <- as.data.frame(df)

# 後々比較を行うことを考えて,5つの地域と月の列を追加
area <- rep(1:5, 60)
month <- rep(1:12, 25)
df.ele <- cbind(area, month, d.ele)

summary(df.ele)
# データが確認出来たらcsvファイルとしてアウトプット
write.csv(df.ele, "67ele-sim.csv")
# Boxplotの作成
# ggplotでの作図はコードも長くなり,とっつきにくいですがより高品質なグラフを作ることができます
library(ggplot2)
library(reshape2)
library(scales)

d.ele.hist <- melt(d.ele) # 一つの変数として扱うために,データを列方向にまとめる
## No id variables; using all as measure variables
ggplot(d.ele.hist, aes(x=variable, y=value)) + theme_bw() +
  annotation_logticks(sides = "l") + # l で左だけ,trblだと上下左右にtickが入る
  scale_y_log10( # y軸を対数軸に
    breaks = 10^(-7:7),
    labels = trans_format("log10", math_format(10^.x)) # y軸の表記を10^x
  ) +
  geom_boxplot(outlier.alpha = 0.6) + # boxplotを追加
  stat_summary(fun.y="mean", geom="point", shape=21, size=3, fill="white") + # 平均値を追加する
  xlab("") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=45, vjust=1, hjust=1)
  )
# 目的のグラフが作れたらファイルを保存する
# とくにWindows環境は透明度を上手く扱えないので,alphaの値を調製した場合はpngファイルで保存する
ggsave("boxplot.png", width=40, height=25, unit="cm", dpi=600)

# {ggplot2}の場合は,ggsave()で画像を保存できる.

環境分析を行った時にまず確認すべき事の一つは,測定値の妥当性です. 認証物質を測定して測定方法に問題が無いか確認しますが,正解が分からない環境試料の場合には,計算間違い一つで重大な影響を与える可能性もあり,より注意を払う必要があります.

  • 濃度レベルがこれまでの知見と大きな違いが無いか?
  • 外れ値はどうか?多すぎたり,偏っていたりしないか?
  • Oddo-Harkins則のように,(地球)化学的な知見から外れていないか?

以上の確認をして,場合によってはスペクトル干渉の有無を再確認し,他の分析条件を採用すべきか検討すべきでしょう.

適切なデータがあれば,以下のようなBoxplotが作成できます.

2017年09月29日ブログのカテゴリー:R, blog, ggplot, 微量元素

平均値の比較

Rとggplot2を使って平均値を比較してみましょう.

使うデータは前回作成した67元素のシミュレーションデータです.⇒ データ

 

d <- read.csv("67ele-sim.csv")

library(ggplot2)

ggplot(
  d, # 使用するデータフレームを指定
  aes(x=factor(area), y=Li) #指定したデータフレームからxとyを指定
  ) +
  geom_boxplot() # Boxplotを描画すると指定

# 背景や軸の色などを設定する
ggplot(d, aes(x=factor(area), y=Li)) + 
  geom_boxplot() +
  xlab("Area") + ylab("Concentration (μg/L)") + #x,y軸名を設定
  theme_bw() + # 基本的なテーマを白黒にする
  theme(
    panel.grid.minor = element_line(linetype = "blank"), #補助線を無し 
    axis.text = element_text(colour = "black"), # 軸のテキストを黒にする
    legend.position="bottom", # 凡例は下に表示
    text=element_text(size=14), # テキストのサイズを14に
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5) #x軸のテキストの位置を調整
  )

# 何のデータか分からないのでラベルを追加する
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_boxplot(aes(fill=factor(area))) +
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]), #ラベルの位置を表示するためのデータフレームを作成
    aes(x=x, y=y, label=label), # 使用するデータを指定
    hjust=0, # ラベルの左端を水平の位置に合わせる
    vjust=1 # ラベルの上端を垂直位置に合わせる
    ) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
# boxplotに平均値を追加する
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_boxplot() +
  stat_summary( # 平均値を追加する
    fun.y="mean", # 算術平均を指定
    geom="point", # 点をプロット
    shape=21, # 縁のある丸
    size=2.5, 
    fill="white"
    ) + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
# バイオリンプロット
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_violin(adjust = 1.2) +
  geom_boxplot(
    width=0.15, #ボックスの幅
    fill="Black", #ボックスの塗りつぶし
    outlier.alpha = 0 #外れ値の透明度を0にして表示しないように
  ) +
  stat_summary(fun.y="mean", geom="point", shape=21, size=2.5, fill="white") + # 平均値を追加する
     xlab("Area") + ylab("Concentration (μg/L)") +
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
# 蜂群図
library(ggbeeswarm) # {ggplot2}のgeom_dotplot()でも描画できるが,{ggbeeswarm}を使った方がy値の違いを表現できる
## Warning: package 'ggbeeswarm' was built under R version 3.3.3
ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_violin(fill="grey90", adjust = 1.5) +
  geom_boxplot(width=0.15, fill="White", outlier.alpha = 0) +
  geom_beeswarm(
    colour="black",
    priority="density",
    cex=1.5,
    size=1
    ) +
    geom_boxplot(width=0.15, fill="White", outlier.alpha = 0, alpha=0.5) +
  stat_summary(fun.y="mean", geom="point", shape=21, size=3, fill="white") + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )

Excelでこれらよりももっと簡単な棒グラフを作成して,x軸を地域間から採取月差を確認したいと思っても,12組のxとyの範囲を設定する必要が出てかなり面倒です.

Rを使って作図するようになると,x軸を他のパラメーターに変更するのはとても簡単に実行できます.

# xのデータをmonthに変更するだけで簡単に月別のデータを表示できる
ggplot(d, aes(x=factor(month), y=Li)) + theme_bw() +
  geom_boxplot() +
  stat_summary( # 平均値を追加する
    fun.y="mean", # 算術平均を指定
    geom="point", # 点をプロット
    shape=21, # 縁のある丸
    size=2.5, 
    fill="white"
    ) + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  xlab("Month") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )

さらに,地域×採取月のようなパターンの比較も簡単に作図できます.

# 地点ごと,採取月ごとの比較を行う

ggplot(d, aes(x=factor(area), y=Li)) + theme_bw() +
  geom_boxplot() +
  stat_summary(fun.y="mean", geom="point", shape=21, size=2.5, fill="white") + 
  geom_label(
    data=data.frame(x=-Inf, y=Inf, label=colnames(d)[4]),
    aes(x=x, y=y, label=label), hjust=0, vjust=1
    ) +
  facet_wrap(~month) +
  xlab("Area") + ylab("Concentration (μg/L)") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )

最後に,Boxplotで地域間の比較を67元素で一度に比較してアウトプットしてみます.

その際に,イオン価数,イオン半径,イオン化傾向など周期表の配置は元素間の関係を把握するために非常に重要です.そこで,周期表を模した配置してみます.これは{geofacet}というパッケージを使うと非常に簡単に作図することができます.

{geofacet}でしようする周期表上のグリッド( Grid_Period_Tbl_Ele (-ARI).csv )をまず作成・インポートします.


library(tidyr)
library(dplyr)
library(scales)

d.ele.melt <- tidyr::gather(d, "area", "month", 4:70) # Longフォーマットに変形
colnames(d.ele.melt) <- c("ID", "area", "month", "name", "Conc") # 列名を変更,元素は"name"にしておくこと
d.ele.melt$area <- as.factor(d.ele.melt$area)
  
library(geofacet)
PTEGrid <- read.csv("Grid_Period_Tbl_Ele (-ARI).csv") # カスタムグリッドは'code'と'name'にする必要がある

ggplot(d.ele.melt, aes(x=area, y=Conc)) + theme_bw() +
  annotation_logticks(sides = "l") + # l で左だけ,trblだと上下左右にtickが入る
  scale_y_log10( # y軸を対数軸に
    breaks = 10^(-7:7),
    labels = trans_format("log10", math_format(10^.x)) # y軸の表記を10^x
  ) +
  geom_boxplot(aes(fill=area)) +
  stat_summary(fun.y="mean", geom="point", shape=21, size=2.5, fill="white") +
  facet_geo(facets = ~ name, grid = PTEGrid, scales="free_y") +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=14), 
    axis.text.x=element_text(angle=0, vjust=1, hjust=0.5)
  )
## You provided a user-specified grid. If this is a generally-useful
##   grid, please consider submitting it to become a part of the
##   geofacet package. You can do this easily by calling:
##   grid_submit(__grid_df_name__)

以下に示すように,各元素の地域別濃度を周期表上に表示することができます.

カーソルを移動させると,右側のウインドウで各元素のデータを確認することができます.

2017年10月30日

シンタックスハイライト

データサイエンス関係のブログを拝見していると,普段Rstudioで表示されるようにコードの視認性が高くなるようなスタイルが設定されています.HPBではRmarkdownで生成したhtmlをそのまま埋め込むことができないので少し放置していました.

このあたりの知識は全くなかったので,よく参考にしているHPのソースコードを見ると○○syntax highlightのような記述が・・・

このキーワードを検索してみたところ以下のサイトが比較的簡単に導入できそうな例を紹介しています.

JavaScriptでシンタックスハイライトする入力エリアを提供するライブラリを探した

このサイトの内容を参考にして,Rmarkdownで生成したhtmlファイルのソースコードからjavascriptとcssファイルを作成して,多元素データの生成の前半部分のコードに適用してみました.

 

library(mvtnorm)
# {MASS}の()では負の相関を扱えない

m = 67 # number of variable
g = 5
n = 60 # sample size
N = g *n

d.cor <- read.csv("cor.csv") 
m.cor <- as.matrix(d.cor[,-1])

GM <- read.csv("GM.csv")
GM <- GM[,-1]

GSD <- read.csv("GSD.csv")
GSD <- GSD[,-1]

#
df <- rmvnorm(n=N, mean=rep(0, m), sigma=m.cor)
# もう一度標準化しておく
df <- scale(
  df, 
  center = TRUE, # 平均値を0にする
  scale = TRUE # 標準偏差を0にする
  )

# for文を使ってセル[i, j]にlog(GSD[,j])を乗じて,log(GM[,j])を足して
# logから元に戻すために10^ の計算をする
for (i in 1:N)  {
  for (j in 1:m) {
    df[i, j] <- 10 ^ ( df[i, j]* log10(GSD[1, j]) + log10(GM[1, j])  )
  }
}

# 変数名を元素名に変更
colnames(df) <- colnames(GM)
d.ele <- as.data.frame(df)

# 後々比較を行うことを考えて,5つの地域と月の列を追加
area <- rep(1:5, 60)
month <- rep(1:12, 25)
df.ele <- cbind(area, month, d.ele)

HPBだけではjavascriptやcssファイルの場所を指定・保存できていないので,FFFTPで別途jsフォルダとcssフォルダを作成してその中に保存しました.

これはHPBのプレビューではコードにカラーが付いていることを確認できないので,何度か試行錯誤しましたが,なんとか上手くいったようです.

これまでの内容もすべて更新しました.モヤモヤしていたのがスッキリしました.これでもう少しブログの更新頻度を上げたいところです.

2017年11月09日

機械学習1 Random forest

我々も研究に機械学習を取り入れ始めています.機械学習(あるいは昨今のAIと思われているもの)というとDeep Learningのイメージが強いですね.Deep LearningはA Neural Network Playgroundというサイトで体験することが出来ます.このサイトではパラメーターを変化させてその挙動を見て学ぶことが出来ます.よく参考にさせていただいている有名なデータサイエンティストであるTJO氏のブログでも紹介されていますし,Deep Learning以外の機械学習の手法との比較をしています.

今回はその記事をなぞった内容で,上記NN Playgroundの二重らせんデータに対して何回かに分けて機械学習を試してみたいと思います.まずは二重らせんデータ生成し,用いました.本来はデータの一部のみを学習用データとするべきでしょうが,今回は全データを学習用データとして,TJO氏のブログにあるようにグリッド内の全点をテストデータとしています.また,私はggplot脳なので,{ggplot2}でアウトプットしていきたいと思います.

学習用データと試験データの作成

n <- 500
x <- seq(0.02, 4, by = 4/n)

noise1 <- rnorm(n, mean=0, sd=1)
noise2 <- rnorm(n, mean=0, sd=1)
noise3 <- rnorm(n, mean=0, sd=1)
noise4 <- rnorm(n, mean=0, sd=1)

L1 <- x + 0.05 * noise1 +0.02
L2 <- x + 0.05 * noise2 +0.02

r <- 0.05

X1 <- (1+r*noise3)*L1 * cos(10*L1/pi)
Y1 <- (1+r*noise3)*L1 * sin(10*L1/pi)

X2 <- -(1+r*noise4)*L2 * cos(10*L2/pi)
Y2 <- -(1+r*noise4)*L2 * sin(10*L2/pi)

df <- rbind(cbind(0, X1, Y1),cbind(1, X2, Y2)) 
df <- as.data.frame(df)
colnames(df) <- c("type", "x", "y")
df$type <- as.factor(df$type)

write.csv(df, "double spiral.csv")

px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")

NN PlaygroudでいうNoiseを加えています.実際はcsvファイルにアウトプットして毎度読み込んでます.{ggplot2}でNN Playgroud風にプロットするとこんな感じになります.

また,テストデータとしてgridという0.02刻みのデータを作成しました.

library(ggplot2)
theme_set(
  theme_bw() +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=12)
  )
)

ggplot(df, aes(x=x, y=y, colour=type)) + 
  geom_point(alpha=0.5) +
  scale_color_manual(values = c("darkorange", "darkblue")) 

Random forestを試してみる

我々は環境分析の研究室なので,環境汚染の度合いとその影響,発生源と環境汚染など目的変数と説明変数の関係性を評価したいというのが解析の目的になることが多いです.

そのため,変数の重要度を評価できるRandom forestという方法をよく使用しています.Random forestについての説明はWikipedia等見てください.

他にもその発展形のXgboostも変数重要度が計算できるので今後取り入れたいと思っているところです.

d <- read.csv("double spiral.csv")
d <- d[,2:4]
d$type <- as.factor(d$type)

px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")
library(randomForest)
## Warning: package 'randomForest' was built under R version 3.3.3
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
fit.rf <- randomForest(type~., d)
set.seed(8800)
tuneRF(d[-1], d$type, doBest=T)
## mtry = 1  OOB error = 0.7% 
## Searching left ...
## Searching right ...
## mtry = 2     OOB error = 1% 
## -0.4285714 0.05
## 
## Call:
##  randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1]) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 0.5%
## Confusion matrix:
##     0   1 class.error
## 0 496   4       0.008
## 1   1 499       0.002
fit.rf <- randomForest(type~., d, mtry=1)

pred.rf <- predict(fit.rf, newdata=grid)
summary(pred.rf)
##     0     1 
## 94210 94661
d.pred.rf <- as.data.frame(cbind(grid, as.numeric(pred.rf)))
colnames(d.pred.rf) <- c("x", "y", "pred")
d.pred.rf$pred <- d.pred.rf$pred - 1 # 0,1にする

predict()関数でテストデータであるgridをd.rfというモデルで予測しています.予測した結果を対応する"x","y"に対して"pred"という変数で保存しています.

 

library(ggplot2)
theme_set(
  theme_bw() +
  theme(
    panel.grid.minor = element_line(linetype = "blank"), 
    axis.text = element_text(colour = "black"),
    legend.position="bottom",
    text=element_text(size=12)
  )
)

ggplot(NULL) + # 異なるデータ数のデータフレームを重ねる時は初めにggplot(NULL) としておく
  geom_raster(aes(x=d.pred.rf$x, y=d.pred.rf$y, fill=d.pred.rf$pred), alpha =0.15) +
  geom_point(aes(x=d$x, y=d$y, colour=d$type), alpha = 0.5) +
  scale_fill_gradientn(colours = c("darkorange", "grey95", "darkblue")) + 
  scale_color_manual(values = c("darkorange", "darkblue")) +
  theme(legend.position = "none") +
  xlab("x") + ylab("y")

Random forestは決定木をベースにした解析法なので,TJO氏も指摘しているようにカクカクとした決定境界になっていますね.しかも,縦と横の線が組み合わさった境界になっています.また,(x,y)=(0,-0.8) では複雑なギザギザな境界が形成されて,yが-2.8付近では横に伸びた青の境界が形成されており,明らかに過学習していてます.

解析に使用する場合には注意したいところです.

2018年03月14日

機械学習2 Support vector machine (SVM)

機械学習1 Random forestにつづいてSupport vector machine (SVM) による分類を試してみます.WikipediaのSVM記事によれば,SVMは境界に最も近いサンプルとの距離(これをmargin)を最大化するように境界を設定して,その計算の際にカーネルトリックと呼ばれるカーネル関数を用いて再生核ヒルベルト空間に写像することで非線形の分類問題にも対応できるようです.機械学習1 Random forestで生成したcsvファイルを作業ディレクトリに保存し,適当に読み込んでください.

d <- read.csv("double spiral.csv")
d <- d[,2:4]
d$type <- as.factor(d$type)

px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")

SVMによる予測モデルの作成と予測

{e1071}を使用してSVMモデルの作成とテストデータの予測を行ってみます.

gammaのデフォルト値は 1/パラメーター数 ですが,その場合だときれいな決定境界が引けなかったのでいくつか検討してみて10倍しています

library(e1071)
## Warning: package 'e1071' was built under R version 3.3.3
set.seed(8800)
gamma = 1/ncol(df)
fit.svm <- svm(type~., d, gamma = gamma*10)

pred.svm <- predict(fit.svm, newdata=grid)
summary(pred.svm)
##      0      1 
##      0 188871
d.pred.svm <- as.data.frame(cbind(grid, as.numeric(pred.svm)))
colnames(d.pred.svm) <- c("x", "y", "pred")
d.pred.svm$pred <- d.pred.svm$pred -1 # 0,1にする

テーマの設定はRandom forestの時と同じです


library(ggplot2)

ggplot(NULL) + # 異なるデータ数のデータフレームを重ねる時は初めにggplot(NULL) としておく
  geom_raster(aes(x=d.pred.svm$x, y=d.pred.svm$y, fill=d.pred.svm$pred), alpha =0.15) +
  geom_point(aes(x=d$x, y=d$y, colour=d$type), alpha = 0.5) +
  scale_fill_gradientn(colours = c("darkorange", "gray90","darkblue")) + 
  scale_color_manual(values = c("darkorange", "darkblue")) +
  theme(legend.position = "none") +
  xlab("x") + ylab("y")

これは非常にうまく分けることが出来ている印象です.しかも境界を「ヒトの手で書いたらこんなものかな~」というところに引けていますし,滑らかな曲線になっています.
我々の研究ではまだ試したことは無かったのですが,ここまでキレイに決定境界が計算できるのであれば研究にも応用したいですね.SV回帰 (SVR) だと回帰係数として変数の重要度を計算できるようなので,今後試してきたいですね.

2018年03月18日

機械学習3 XGBoost

Xgboost (eXtreme Gradient Boosting) は,Random Forest (決定木のアンサンブル学習) とGradient Boosting (重み付きアンサンブル学習) とを組み合わせたアルゴリズムだそうです.

@aaatsushi_bb によると,「世界的なデータ解析コンテストサイト"Kaggle"では、2015年に出された問題29問のうち、なんと17問がXGBoostを用いたモデルが1位となっています。」ということで非常に人気のある分類器です.たくさんのblogで紹介されています.

# データの読み込みとテストデータとしてグリッドを作成
d <- read.csv("double spiral.csv")
d <- d[,2:4]
d$type <- as.factor(d$type)

px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")
# XGBoostモデルの作成と予測
library(xgboost)
library(caret)
library(plyr); library(dplyr)
library(tidyverse)
# {caret}を使用してパレメーターのチューニングを行う
set.seed(88)
# trainControlを設定
ctrl <- trainControl(
  method = "cv", # cross validation
  number = 5, # 5-fold
  selectionFunction = "best"
  )

# 格子探索用のグリッド
tuneGrid <- expand.grid(
  nrounds = 500,
  eta = seq(0.01, 0.4, 0.01), # 0.01~0.40まで, 0.01刻み
  max_depth = 3:8,
  gamma = 0,
  colsample_bytree = 1,
  min_child_weight = 1,
  subsample = 1
  )

# xgboostのTuning
start <- Sys.time()

xgb.tune <- train(
  d[,2:3],
  d[,1],
  method = "xgbTree",
  metric = "Accuracy",
  trControl = ctrl,
  tuneGrid = tuneGrid
  )

end <- Sys.time()
time <- end-start
time
## Time difference of 29.23513 mins
xgb.tune$bestTune
##     nrounds max_depth eta gamma colsample_bytree min_child_weight
## 119     500         7 0.2     0                1                1
##     subsample
## 119         1
# Tune結果はデータフレームとして保存されている ggplot(xgb.tune$results, aes(x=eta, y=Accuracy, colour=factor(max_depth))) + geom_line(size=1) + geom_point(aes(shape=factor(max_depth), fill=factor(max_depth)), size=3, alpha=0.75)

max_depthが7,etaが0.2の時に最もAccucacyが高かったので,これらのパラメーターを使用してmodelを作成します.

ちなみに,max_depthとetaのデフォルト値はそれぞれ,6, 0.3です.max_depthが深くなるほど,etaが小さくなるほど過学習になるようです.チューニングしなくても似たような結果が得られるので,そこまで過学習しているというわけでは無いと思われます.

 

# Tuning結果を用いてXGboostを実行
label <- as.numeric(d[,1]) - 1
max_depth <- xgb.tune$bestTune$max_depth
eta <- xgb.tune$bestTune$eta

model.xgb <- xgboost(
  data=as.matrix(d[, 2:3]),
  label=label,
  num_class=2,
  params=list(max_depth = max_depth, eta = eta, objective="multi:softmax"), #  
  eval_metric="mlogloss",
  nrounds=1000
  )
## [1]  train-mlogloss:0.625895 
## [2]  train-mlogloss:0.565360 
## [3]  train-mlogloss:0.519274 
## [4]  train-mlogloss:0.404485 
## [5]  train-mlogloss:0.323599 
長いので中略

## [996]    train-mlogloss:0.002423 
## [997]    train-mlogloss:0.002422 
## [998]    train-mlogloss:0.002421 
## [999]    train-mlogloss:0.002421 
## [1000]   train-mlogloss:0.002420
print(model.xgb)
## ##### xgb.Booster
## raw: 992.4 Kb 
## call:
##   xgb.train(params = params, data = dtrain, nrounds = nrounds, 
##     watchlist = watchlist, verbose = verbose, print_every_n = print_every_n, 
##     early_stopping_rounds = early_stopping_rounds, maximize = maximize, 
##     save_period = save_period, save_name = save_name, xgb_model = xgb_model, 
##     callbacks = callbacks, num_class = 2, eval_metric = "mlogloss")
## params (as set within xgb.train):
##   max_depth = "7", eta = "0.2", objective = "multi:softmax", num_class = "2", eval_metric = "mlogloss", silent = "1"
## xgb.attributes:
##   niter
## callbacks:
##   cb.print.evaluation(period = print_every_n)
##   cb.evaluation.log()
##   cb.save.model(save_period = save_period, save_name = save_name)
## niter: 1000
## evaluation_log:
##     iter train_mlogloss
##        1       0.625895
##        2       0.565360
## ---                    
##      999       0.002421
##     1000       0.002420
summary(model.xgb)
##                Length  Class              Mode       
## handle               1 xgb.Booster.handle externalptr
## raw            1016155 -none-             raw        
## niter                1 -none-             numeric    
## evaluation_log       2 data.table         list       
## call                15 -none-             call       
## params               6 -none-             list       
## callbacks            3 -none-             list
pred.xgb <- predict(model.xgb, as.matrix(grid))

d.pred.xgb <- cbind(grid, as.numeric(pred.xgb))
d.pred.xgb <- as.data.frame(d.pred.xgb)
colnames(d.pred.xgb) <- c("x", "y", "pred")


# {ggplot}で予測結果をグラフに
ggplot(NULL) + # 異なるデータ数のデータフレームを重ねる時は初めにggplot(NULL) としておく
  geom_raster(aes(x=d.pred.xgb$x, y=d.pred.xgb$y, fill=d.pred.xgb$pred), alpha =0.15) +
  geom_point(aes(x=d$x, y=d$y, colour=d$type), alpha = 0.5) +
  scale_fill_gradientn(colours = c("darkorange", "grey95", "darkblue")) + 
  scale_color_manual(values = c("darkorange", "darkblue")) +
  theme(legend.position = "none") +
  xlab("x") + ylab("y")

基本的には縦・横の境界を組合わせた境界で,その印象はやはりRandom forestと似ていますね.やはり決定木をベースしたモデルはこのようなパターンになってしまうようです.しかし,Random forestとXGBoostの結果を細かく見ると違いがあるのが分かります.x=-2.5, y=2.0あたりを少し拡大して比較してみると,

XGBoostの方がRandom forestよりも,細かな境界の長さが長くなる傾向があります.XGBoostは上手く分類した場合に重みが付く一方で,Random forestは決定木の出力の多数決で評価することから,このような違いが生じると思われます.

Random forestは,一部ではあるものの,決定木をベースでありながら非線形っぽい境界が引けています.

Kaggleでの実績から考えると,XGBoostは非常に魅力的な分類器なんですけれども,非線形的な分類が必要な場合にはRandam forestよりも悪くなる可能性があるんじゃないかと思います.こうやって実際に動かして比較してみることは非常に重要ですね.

2018年04月01日

機械学習4 Deep Neural Network-1

機械学習シリーズもようやくDeep neural networkを紹介できるところまで来ました.Rで実行する場合,いくつかあるようですが,有名どころでは{mxnet}と{keras}があります.まずは,アマゾンが開発した{mxnet}で,古典的なDNN, 誤差逆伝播つき多層パーセプトロンを試していきたいと思います.

インストールに関しては

を参考にしました.また,基本的なRコードはTJO氏のブログを参考にしました.

隠れ層4, 各層のニューロンはそれぞれ5の多層パーセプトロンが,NN Playgroundで上手くいくことが多いかなという感じなので,そのように設定します.

d <- read.csv("double spiral.csv")
d <- d[,-1]
d$type <- as.factor(d$type)

sd(d$x)
## [1] 1.691215
sd(d$y)
## [1] 1.597088
px <- seq(-4.7, 4.7, 0.02)
py <- seq(-4, 4, 0.02)
grid <- expand.grid(px, py)
names(grid) <- c("x", "y")

z.grid <- grid
z.grid$x <- (grid$x - mean(d$x))/sd(d$x)
z.grid$y <- (grid$y - mean(d$y))/sd(d$y)
# Deep Neural Network
# 古典的なDNN, 誤差逆伝播つき多層パーセプトロン
library(mxnet)
## Warning: package 'mxnet' was built under R version 3.4.0
train.x <- data.matrix(scale(d[,-1]))
# train.y <- as.numeric(df[, 1])
train.y <- d[, 1]
test <- data.matrix(z.grid) # 数値でもファクターでも同じ

data <- mx.symbol.Variable("data")
fc1 <- mx.symbol.FullyConnected(data, name = "fc1", num_hidden = 5)
act1 <- mx.symbol.Activation(fc1, name = "relu1", act_type = "relu")

# 'relu': Rectified Linear Unit, :math:'y = max(x, 0)' 
# 'sigmoid': :math:'y = \frac11 + exp(-x)' 
# 'tanh': Hyperbolic tangent, :math:'y = \fracexp(x) - exp(-x)exp(x) + exp(-x)' 
# 'softrelu': Soft ReLU, or SoftPlus, :math:'y = log(1 + exp(x))'

fc2 <- mx.symbol.FullyConnected(act1, name = "fc2", num_hidden = 5)
act2 <- mx.symbol.Activation(fc2, name = "relu2", act_type = "relu")
fc3 <- mx.symbol.FullyConnected(act2, name = "fc3", num_hidden = 5)
act3 <- mx.symbol.Activation(fc3, name = "relu3", act_type = "relu")
fc4 <- mx.symbol.FullyConnected(act3, name = "fc4", num_hidden = 5)
#act4 <- mx.symbol.Activation(fc4, name = "relu4", act_type = "relu")
#fc5 <- mx.symbol.FullyConnected(act4, name = "fc5", num_hidden = 2)
softmax <- mx.symbol.SoftmaxOutput(fc4, name = "softmax")

devices <- mx.cpu()

mx.set.seed(123)

dnn1 <- mx.model.FeedForward.create(
  softmax,
  X = train.x, 
  y = train.y, 
  ctx = devices, 
  num.round = 1000, 
  array.batch.size = 25, 
  learning.rate = 0.03, 
  momentum = 0.8,  
  eval.metric = mx.metric.accuracy, 
  initializer = mx.init.uniform(1), # The scale of uniform distribution
  optimizer = "sgd", 
  array.layout = "rowmajor", 
  allow.extra.params = TRUE,
  epoch.end.callback = mx.callback.log.train.metric(100)
  )
## Start training with 1 devices
## [1] Train-accuracy=0.468717948717949
## [2] Train-accuracy=0.577
## [3] Train-accuracy=0.582
## [4] Train-accuracy=0.579
## [5] Train-accuracy=0.586
長いので中略
## [996] Train-accuracy=0.998
## [997] Train-accuracy=0.998
## [998] Train-accuracy=0.998
## [999] Train-accuracy=0.998
## [1000] Train-accuracy=0.998
pred.dnn1 <- t(predict(dnn1, test, array.layout = "rowmajor"))

d.pred.dnn1 <- cbind(grid, as.numeric(pred.dnn1[, 3]))
d.pred.dnn1 <- as.data.frame(d.pred.dnn1)
colnames(d.pred.dnn1) <- c("x", "y", "pred")

library(ggplot2)
ggplot(NULL) + # 異なるデータ数のデータフレームを重ねる時は初めにggplot(NULL) としておく
  geom_raster(aes(x=d.pred.dnn1$x, y=d.pred.dnn1$y, fill=d.pred.dnn1$pred), alpha =0.25) +
  geom_point(aes(x=d$x, y=d$y, colour=d$type), alpha = 0.5) +
  scale_fill_gradientn(colours = c("darkorange", "grey95", "darkblue")) + 
  scale_color_manual(values = c("darkorange", "darkblue")) +
  theme(legend.position = "none") +
  xlab("x") + ylab("y")

一応,それっぽい境界を作成することが出来ました.Trainingデータに対する予測精度は0.998でした.

こうみてみるとDNNの学習も基本的には,インプットであるx,yだけに依存するので線形分離していく感じですね.Random forestやXGBoostと異なるのは,境界が斜線にもなるところでしょうか.

最終的なアウトプットはまあ,なんとかそれっぽくなったのですが,ここまで来るのは大分苦労しました.初めのうちは,NN Playgroundと同じ設定にしても,全く学習が進まないこともあって,チューニングには試行錯誤しました.

ある程度学習するのが分かってからきちんとログを取ったので参考にして下さい.

DNN_comp <- read.csv("DNN/DNN comp.csv")
rownames(DNN_comp) <- DNN_comp$X
DNN_comp <- DNN_comp[,-1]

knitr::kable(DNN_comp)
Train.accuracy Hidden.layers neurons Activation num.round lerning.rate batch.size momentum initializer optimizer Appendix
0.6617 4 5,5,5,5 ReLU 3000 0.002 16 0.95 mx.init.uniform(1) sgd
0.6438 4 5,5,5,5 ReLU 3000 0.010 16 0.95 mx.init.uniform(1) sgd
0.5823 4 5,5,5,5 ReLU 1000 0.030 16 0.90 mx.init.uniform(1) sgd
0.8570 4 5,5,5,5 ReLU 1000 0.030 20 0.80 mx.init.uniform(1) sgd
0.8390 4 5,5,5,5 ReLU 3000 0.030 20 0.80 mx.init.uniform(1) sgd
0.8600 4 5,5,5,5 ReLU 1000 0.030 20 0.70 mx.init.uniform(1) sgd
0.8060 4 5,5,5,5 ReLU 1000 0.050 20 0.70 mx.init.uniform(1) sgd
0.8105 4 5,5,5,5 ReLU 1000 0.030 16 0.80 mx.init.uniform(1) sgd
0.9980 4 5,5,5,5 ReLU 1000 0.030 25 0.80 mx.init.uniform(1) sgd
0.7570 4 5,5,5,5 ReLU 1000 0.010 25 0.80 mx.init.uniform(1) sgd
0.6250 4 5,5,5,5 ReLU 1000 0.050 25 0.80 mx.init.uniform(1) sgd
0.9360 4 5,5,5,5 ReLU 3000 0.010 25 0.80 mx.init.uniform(1) sgd
0.9350 4 5,5,5,5 ReLU 5000 0.010 25 0.80 mx.init.uniform(1) sgd
0.9863 4 5,5,5,5 ReLU 1000 0.030 32 0.80 mx.init.uniform(1) sgd
0.8027 4 5,5,5,5 ReLU 1000 0.030 64 0.80 mx.init.uniform(1) sgd
0.9580 4 5,5,5,5 ReLU 1000 0.030 25 0.70 mx.init.uniform(1) sgd
0.8610 4 5,5,5,5 ReLU 1000 0.040 25 0.70 mx.init.uniform(1) sgd
0.9080 4 5,5,5,5 ReLU 1000 0.020 25 0.70 mx.init.uniform(1) sgd
0.9500 4 5,5,5,5 ReLU 3000 0.030 25 0.80 mx.init.uniform(1) sgd 2300ぐらいで下がる
0.8900 4 5,5,5,5 ReLU 3000 0.030 25 0.75 mx.init.uniform(1) sgd
0.8630 4 5,5,5,5 ReLU 3000 0.030 25 0.85 mx.init.uniform(1) sgd
0.8490 4 5,5,5,5 ReLU 3000 0.010 25 0.85 mx.init.uniform(1) sgd
0.9980 4 5,5,5,5 ReLU 5000 0.030 25 0.80 mx.init.uniform(1) sgd

普段は表をアウトプットすることはあまりないのですが(Excelで充分かなという場合が多い),再現性の高い解析結果にするためにも,解析に用いるデータセットをこのように一緒にknitしてアウトプットしておく方がよさそうですね.

本当は{DT}を使って,インタラクティブな表にしたかったのですが,埋め込むことが出来ませんでした.このようにアウトプットすることが可能です.-> リンク

追記1 (2018/04/10)

活性化関数ReLUについてとReLU一族【追記あり】| @hokekiyoo

上記のブログでは活性化関数の違いによる挙動が比較されています.直線を組合わせた決定境界になるのは活性化関数にReLUを使った際の特徴のようです.

2018年04月02日