2018年5月29日火曜日

Rいろは:グラフ作成ggplot2編

###### ggplot2とは? #################
# 標準仕様でもそれらしい科学的グラフを作ってくれるRパッケージ
# (plot関数の場合、相当なアレンジをしないと使い物にならない)
# 色分けや複合グラフ、重ね合わせなどに強い、応用の利きやすいメリットがある
# さらにcowplotパッケージを適用して、より学術用に適したグラフ作成を紹介します

# 項や足し算形式でグラフ要素を追加する形態になっている。例えば、散布図、棒グラフのようなグラフ形式や、回帰直線など。
# 項の指定方法などが(plot関数と違って)英単語を元にしているものが多く連想しやすい

# 使い方の早見表も参考になります

install.packages("ggplot2", dep=T) # 未インストールの場合
install.packages("cowplot", dep=T) # 未インストールの場合
install.packages("reshape", dep=T) # 未インストールの場合
require(ggplot2) # 呼び出し

###### まず定番のdata"iris"を用いて、originalのggplot2で描いてみます #########

data(iris) # iris呼び出し
head(iris) # 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



##### 散布図 ####################
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)

# ggplot(data=iris, # データ元の指定
#  aes(x=Sepal.Length, y=Sepal.Width, # x, yに用いる変数の指定(aes()内で変数を指定)
#      fill=Species)) # Speciesによって異なる色を使用
# + geom_point( # + geom_xxx: xxxの部分でグラフ種類を指定
#  pch=21,  # pch=21: 塗りと枠に異なる色を適用可能(21~25) # 全種に適用する場合、aes外に記述
#  size=2) # プロットサイズ(少し大きくしました)

# 以下でも同様のグラフになります(fill=Speciesを後半に記述)
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2)

# このままでは背景がうるさく、学術用に不向き
require(cowplot) # cowplotパッケージを呼び出し、同じようにグラフを作成します
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)
# 背景が白地のスッキリしたものに変わりました。文字のサイズもちょうどいいです。



##### 箱ひげ図 ####################
ggplot(data=iris, aes(x=Species, y=Sepal.Width, fill=Species)) + geom_boxplot()



##### 棒グラフ+SD ####################
# サンプルデータの準備
MEAN1 <- with(iris, tapply(Sepal.Width, Species, mean)) # 各種の平均("対象データ", "グループ分け", "関数名")
SD1 <- with(iris, tapply(Sepal.Width, Species, sd)) # 各グループの標準偏差
iris2 <- data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1)
# グラフの作成
ggplot(data=iris2, aes(x=Species, y=me, fill=Species)) + geom_bar(stat="identity") + geom_linerange(aes(ymax=me + sd, ymin=me - sd))



##### 2要因の棒グラフ + SD ####################
# サンプルデータの準備
iris3 <- rbind( # 1.2倍したサンプルデータを追加
data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1, treat="C"),
data.frame(Species=names(MEAN1), me=1.2*MEAN1, sd=1.2*SD1, treat="T") )
# グラフの作成
ggplot(data=iris3, aes(x=Species, y=me, fill=treat)) + geom_bar(stat="identity", position="dodge") + geom_linerange(aes(ymax=me + sd, ymin=me - sd), position=position_dodge(width=0.9))
# position, dodge: 接したグラフを作成するためのkeyword



##### ヒストグラム ####################
ggplot(iris, aes(x=Sepal.Width, fill=Species)) + geom_histogram()



##### 時系列データのグラフ ####################
data(airquality) # data"airquality"を使用
head(airquality) # データの中身をチェック
  #Ozone Solar.R Wind Temp Month Day
#1    41     190  7.4   67     5   1
#2    36     118  8.0   72     5   2
#3    12     149 12.6   74     5   3
#4    18     313 11.5   62     5   4
#5    NA      NA 14.3   56     5   5
#6    28      NA 14.9   66     5   6

# 時間形式に変換
airquality$Time <- strptime(with(airquality, paste(2018, Month, Day, sep="-")), "%Y-%m-%d")
ggplot(airquality, aes(x=Time, y=Ozone)) + geom_point(pch=21, cex=2, fill="purple") + geom_line()



##### 時系列データのグラフ(複数系列) ####################
require(reshape) # 要インストール
airquality2 <- melt( # データ構成を変換
airquality[,c("Ozone","Temp","Time")], id.var="Time")
head(airquality2)
        #Time variable value
#1 2018-05-01    Ozone    41
#2 2018-05-02    Ozone    36
#3 2018-05-03    Ozone    12
#4 2018-05-04    Ozone    18
#5 2018-05-05    Ozone    NA
#6 2018-05-06    Ozone    28

# 2パネル構成
base <- ggplot(airquality2, aes(x=Time, y=value)) # ごちゃごちゃしてきたので分割します
poi <- geom_point(pch=21, cex=2, aes(fill=variable)) 
lin <- geom_line(aes(colour=variable)) 
base + poi + lin + facet_wrap(~ variable, scale="free_y", ncol=1) # variableでパネルを分割




###### 散布図+回帰直線 ##################
summary(iris) # 値の幅をチェック
  #Sepal.Length    Sepal.Width     Petal.Length    Petal.Width          Species
 #Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100   setosa    :50
 #1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300   versicolor:50
 #Median :5.800   Median :3.000   Median :4.350   Median :1.300   virginica :50
 #Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199                
 #3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800                
 #Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500                

md1 <- lm(Sepal.Length ~ Sepal.Width, iris[iris$Species=="setosa",]) # 1種に限定
coef(md1)
#(Intercept) Sepal.Width
  #6.5262226  -0.2233611

# 散布図 + 回帰直線
x.interval <- seq(2,5,0.1)
est.Y <- predict(md1, newdata=data.frame(Sepal.Width=x.interval))
fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y)
# fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=coef(md1)[1] + coef(md1)[2]*x.interval) # 上と同じ
base2 <- ggplot(iris[iris$Species=="setosa",], aes(x=Sepal.Width, y=Sepal.Length)) 
poi2 <- geom_point(fill="turquoise", pch=21, size=2) 
lin2 <- geom_line(data=fitted, colour="turquoise", lwd=1) # ここは"data="を明記する必要あり
base2 + poi2 + lin2



##### 散布図 + 回帰直線 + 信頼区間 ##################
est.Y2 <- predict(md1, newdata=data.frame(Sepal.Width=x.interval), interval="confidence")
head(est.Y2)
       #fit      lwr      upr
#1 4.019981 3.753101 4.286860
#2 4.089030 3.839589 4.338470
#3 4.158079 3.925981 4.390177
#4 4.227128 4.012251 4.442004
#5 4.296177 4.098369 4.493984
#6 4.365226 4.184291 4.546160

fitted2 <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y2[,1], upperCI=est.Y2[,3], lowerCI=est.Y2[,2])
band <- geom_ribbon(data=fitted2, aes(ymin=lowerCI, ymax=upperCI), alpha=0.5, fill="turquoise", linetype="blank") # alpha: 半透明(0~1)
base2 + poi2 + band + lin2 # 重ねたい順で



##### 散布図 + 回帰直線(複数グループ) ##################
x.interval3 <- seq(1,5,0.1)
new.dat <- data.frame( # 種数分、繰り返す
Species=rep(unique(iris$Species),each=length(x.interval3)), 
Sepal.Width=rep(x.interval3, length(unique(iris$Species))))
md3 <- lm(Sepal.Length ~ Sepal.Width*Species, iris) # 全種を対象、Speciesによる違い(ANCOVA)
est.Y3 <- predict(md3, newdata=new.dat)
fitted3 <- data.frame(Species=new.dat$Species, Sepal.Width=new.dat$Sepal.Width, Sepal.Length=est.Y3) # Species列を追加
base3 <- ggplot(iris, aes(x=Sepal.Width, y=Sepal.Length)) 
poi3 <- geom_point(aes(fill=Species), pch=21, size=2) 
lin3 <- geom_line(data=fitted3, aes(colour=Species), lwd=1) # ここは"data="を明記する必要あり
base3 + poi3 + lin3




##### 複数タイプのグラフの組み合わせ(cowplot利用) ##################
g1 <- ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2) 
g2 <- ggplot(data=iris, aes(x=Species, y=Petal.Width)) + geom_boxplot(aes(fill=Species)) # 箱ひげ図
+ theme(axis.text.x = element_text(angle=90, vjust=0.5)) # label方向

g12 <- plot_grid(g1, g2, labels="auto", align="h") 
# labels="auto":小文字でラベリング(大文字は"AUTO")
# 直接入力の例:labels=c("a","b")
# align="h":図の端を揃える("h", "v", "hv")

#### グラフの調整(例:凡例を統一、サイズ調整)
g1n <- g1 + theme(legend.position="none") # 一方の凡例消去
bd <- panel_border(colour=1, size=1) # 外枠追加
g1n2 <- plot_grid(g1n + bd, g2+ bd, labels="auto", 
   rel_widths=c(1, 1.4), # 横幅を1:1.4に変更
   scale=0.9, # パネル間の間隔を広げる
   vjust=0.2) # ラベル位置(a, b)の調整
 
 
 
#### ファイルの保存 ##############################
ggsave("directry/filename.pdf", g1n2, height=4, width=10, unit="cm") # 保存


#### 水平線・垂線 ##############################
abl <- geom_hline(yintercept=0, linetype="dotted") # 水平 ("dotted": 点線)
abl2 <- geom_vline(xintercept=0, linetype="dotted") # 垂直


#### 軸とラベルの調整 ############################
x.ax <- scale_x_continuous(breaks=c(1:5), limits=c(1, 5), expand=c(0.001, 0.001), name="Sepal width (mm)")
# ex.
base3 + poi3 + lin3 + x.ax


#### 色の指定 ##############################
# グレースケール
set.fillG <- scale_fill_grey(start=0.5, end=1) # グレー 〜 白の場合(塗りつぶし色)
# "fill"の部分を"colour"にすると線の色の指定("shape": 形状、

set.fillC <- scale_fill_manual(values=c("turquoise", "salmon", "blue")) # 塗りつぶし色
set.colC <- scale_colour_manual(values=c("turquoise", "salmon", "blue")) # 線の色
# ex.
base3 + poi3 + lin3 + set.fillC + set.colC

#### cf. 色コードの選択 ##########################
"見やすい"色の組み合わせ例の参照サイト(ColorBrewer)
fill="#2ca25f" のようにしてカラーコードで色の指定が可能



#### 90° 回転(グラフを横向きにする)#######################
 + coord_flip()


#### 凡例の調整 ##############################
+ theme(legend.position="none") # 凡例消し
+ theme(legend.position="bottom") # ラベルの右揃え(cf. vjust)

# 長すぎるラベルを改行する( \n が使える)
scale_x_continuous(name = "Length\n (cm)")

# 上付き、ギリシャ文字など
expression(m^2) 
expression(Delta) # Δ



#### cowplotの調整 ##############################
+ panel_border(colour=1, size=1) # defaultの線・文字色がなぜかグレーなので黒に変更
+ theme_cowplot(font_size = 12) # 文字が大きく感じる場合はサイズ指定

# 日本語の利用(ただしアウトライン化されてしまう模様 )
ggsave(..., family="Japan1GothicBBB",...) 

2017年3月30日木曜日

Rで機種依存文字の混じったデータファイルを読み込む方法

Rへ読み込む際に非常に厄介な機種依存文字(=環境依存文字)、データを作成する際には絶対に避けてほしい代物ですが、丸に数字など、日本のあらゆる文書で愛されており、根絶は非常に困難です。

見渡せる程度の小さなファイルならば手作業で特定し、検索・置換で除去するのですが、巨大なファイルの場合、いったいどんな文字が悪さをしているのか見つけることさえできない場合がありますね。

今回、到底手作業で対処できないサイズのファイルに取り組んでいた際にいい方法を思いついたので挙げておきます(小ネタです)。

1) まず、機種依存文字のファイルは.xlsxファイルとして保存しておきます(ここでは"データファイル名.xlsx"とします)。

2) 通常、read.csv関数などでデータファイルを読み込むところを、openxlsxパッケージのread.xlsx関数を使用します(要インストール、install.packages("openxlsx") などで)。

3) 以下で読み込むことができます、fileEncodingの指定も特に不要の様子。
data <- read.xlsx("データファイル名.xlsx")

cf. もちろんファイル出力にも対応しています(outputというデータフレームを出力する場合の例)write.xlsx(output, "出力ファイル名.xlsx")

Excelファイルを直接読むのはまだ抵抗はあるのですが、気にしなければ汎用性のある方法なので便利です。あとはこのパッケージがちゃんと存続してくれることを祈るばかりですが。なお、他にもExcelファイルを読めるパッケージは複数存在しているようですが、他のは基本的にJavaに依存しているようなので却下しました。

cf. 機種依存文字が混じっているファイルをread.csvで読もうとすると、"In scan(file, what, nmax, sep, dec, quote, skip, nlines, na.strings,  : 入力コネクション 'ファイル名.csv' に不正な入力がありました "のようなエラーメッセージが出る。
Rじゃなくとも、メールの送受信でも文字化けの温床となります。データ解析に限らず、トラブルの原因にしかなりません。

2016年12月21日水曜日

macOS10.12 (Sierra) でWinBUGSを動かす (Wine, R2WinBUGS使用)

メインマシンをクリーンインストールする羽目に会いました。今度こそWinBUGSは卒業したくはあるのですが、依然StanやJAGSでは通らないコードも残っておりまだ手放せずにいます。早くStanが離散パラメータを直接扱えるようになってほしい…

以下は、いったんMountain lion (OSX10.8)をクリーンインストール → Sierraにアップグレードしたマシンにインストールした場合の報告例です。El Capitan(OSX10.11)の時とほぼ同様です。

まず、XQuartzとHomebrew(パッケージ化されていないアプリを容易にインストールするための補助ツール?)を入れ、Homebrewを用いてWine(非Windows OS上でWin専用アプリケーションを実行する環境)をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。

cf. 以下の「ターミナル」の使い方:
「アプリケーション」→「ユーティリティ」にある「ターミナル」を立ち上げる。

コンピュータ名:~ ユーザ名$

このドルマーク $ の後にコマンドを打っていく。
なお、インストールに関わるところでパスワードを求められるが、その都度、自分のアカウントのパスワードを入れる。
(以降、パスワードを入れる作業は説明を省略)




<以下、作業手順(もしかするとSierraでは1〜4の行程は不要かもしれない>
0)実行環境
・Mountain lion (OSX10.8)をクリーンインストール → SierraにアップグレードしたiMac


1) 下準備の開始、/usr/local/フォルダを作る、ロックをいったん外して操作をするという動作をするのですが、その前にシステムの基本的なセキュリティを一旦外します。推奨されていない動作だということをお忘れなく。

リカバリモード(⌘+R を押しながら起動)で起動し、ターミナルを立ち上げる


2) ターミナルの $ マークの直後に、下記のコードを打ち込む(セキュリティを外す作業)これはコピペでOK、以降も同様。

csrutil disable


3) 通常の再起動をする


4) ターミナルに下記を打ち込む(改行されて見えているだろうが、改行無しで打ち込む)

sudo mkdir /usr/local && sudo chflags norestricted /usr/local && sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local

ただし、今回の実行環境では、このディレクトリは既に存在すると言われた。なので、この危なっかしい工程はもはや省略できるかもしれない。


5) 再度、リカバリモードで起動


6) 下記のコードを打ち込む(セキュリティを元に戻す)

csrutil enable


7) 今一度、通常の再起動をする


8) Xcodeのインストール、ターミナルに以下を打ち込む

xcode-select --install


9) homebrewをインストールする

ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

(せいぜい数分でインストールは終了するはず)


10) XQuartzのインストール。もはやアプリとして付属しなくなった。コードでインストール可能であるようなのでインストール。

brew cask install xquartz


11) Wineのインストール、下記のコードをターミナルに打ち込む。

brew install wine


12) 回線速度にもよるが、良好な環境では1時間程度でインストールは終わるだろう。これが終わったら、次のコマンドを打つ。

winecfg


13) XQuartzのウインドウが表示され、"Gecko package"が必要だから入れてもいいか?と聞かれた。表示されているInstallボタンで承諾するとインストールされる。

しばらく処理がされた後、XQuartzからWineの環境設定のようなウインドウが表示される。単に一番下のOKをクリックすればよい。


14) WinBUGSのインストール
パッチなどを当ててある展開済みのWinBUGSフォルダを用意する。Windows7、Windows8へのインストールも現在はこのやり方で行くしかないことを考えれば、Macでも同様にすればいいだろう。適用済みのWinBUGSも公開されている。

下記のコードで不可視フォルダにあるProgram Filesフォルダを開く

open ~/.wine/drive_c/Program\ Files/

ここへWinBUGSフォルダを入れればインストール完了

15)RからWinBUGSを実行する。下記のような単純なサンプルコードで試してみる。Wineを経由するので、bugs()内にそのためのコードがたくさん必要。


# R2WinBUGSのインストールをお忘れなく
require(R2WinBUGS)
# 真の値は、a=3, b=2, sd=1
X <- c(1:100)
Y <- rnorm(100, mean=(3 + 2*X), sd=1)
data <- list(X=X, Y=Y)
inits <- function() list(a=0, b=0, tau=1) 
parameters <- c("a", "b", "sigma")

model <- function() {
a ~ dnorm(0, 1.0E-6)
b ~ dnorm(0, 1.0E-6)
tau ~ dgamma(1.0E-2, 1.0E-2)
for (i in 1:100) {
Y[i] ~ dnorm(mean[i], tau)
mean[i] <- a + b*X[i] }
sigma <- 1/sqrt(tau)
}
modelpath <- file.path(tempdir(), "model.bug")
write.model(model, modelpath)

mcmc <- bugs(
data=data, inits=inits, parameters=parameters, model.file=modelpath, 
n.chains=3, n.iter=5000, debug=T,
working.directory=NULL, clearWD=T, useWINE=T, newWINE=T,
WINE="/usr/local/bin/wine", WINEPATH="/usr/local/bin/winepath")

print(mcmc) # ちゃんと真の値(a=3, b=2, sigma=1)が推定できたかチェックしよう

# 今回、opt/localではなくusr/localにパスを通すよう変更する必要が出た。以前のWINE、WINEPATHは/opt/local/bin/になっていたが、ここは/usr/local/bin/に変更していることに注意。


17)まだR上で下記のエラーコードが出るが、これはこちら(http://ggorjan.blogspot.jp/2008/10/runnning-r2winbugs-on-mac.html)によると害のないエラーコードらしい。要は推定計算さえ無事に行われていればよいだろう。
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject no class object {0003000a-0000-0000-c000-000000000046} could be created for context 0x3
err:ole:CoReleaseMarshalData IMarshal::ReleaseMarshalData failed with error 0x8001011d


2016年11月20日日曜日

潜水から無事に帰還する装備について考える

海中での潜水調査・作業をよくやる人の中には、ヒヤッとした経験のある人も多いだろう。私自身もかつてBC(浮力調節)ジャケットの弁のトラブルで強制浮上になってしまったことがあり、それ以降はとくに気になるようになってきた。

先週、知っている人の中で重大な潜水事故が起こってしまった。まったくもって今更ではあるけれど、ここでは海中または浮上後の海上でトラブルが起こった際に無事に帰還するための装備について考えてみたい。ご意見も歓迎します。

まず、流されてしまったか何かで、岸や船から遠い地点に浮上してしまった場合、BCに付属している笛を吹く程度では、よほど凪でもなければ気づいてもらえないだろう。他に入手しやすいものとしては、海面上に長く伸びるフロート(レギュレータからのエアを入れて膨らませられる)があるだろう。私も念のため持っているけれど、しかし実際これで十分に目立つのだろうか。

他の手段…例えばスマホを完全防水で携帯して、浮上したら電話するというのはどうだろうか。場合によってはアカウントを利用して、PCからスマホを探すことも可能な気がするのだけれど。

スマホに近いものとして、ココセコムというのがあるそうだ。子供やお年寄りに持たせておいて携帯電話網でマップ上で追跡が可能、料金的にもスマホ利用より安そうです。

無線を使うという手段は電波法上かなり難しいようだが、合法的にダイビングで使用できる製品もあるようだ(トランスポンダ SEAKER_1
船からのダイビングという形態を前提としているようだけれど、通常使用(自分の母船に連絡)と、緊急使用(母船とはぐれている場合に無線を受けられるすべての船に緊急発信が可能)の2通りの発信が可能というのが優れている。距離的には20 kmまで届くようなので頼もしい。

(2016.12.16 追記)意識があって海上で操作することが前提ですが、その条件ではかなり有効な製品が出たようです(個人用救難信号発信機 PLB)。精度100mで全世界対応。機器そのものは今年3月に日本で認可、防水ケースが先月末に発売された模様。無線局として申請・許可を得る必要があるなど取り扱いは注意ではあります。

では、意識がないけれど海面には浮かんでいる場合どうするか。上記の携帯を利用した方法が可能であったら使えそうではあるけれど、岸から遠くへ流されてしまっていたとしたら(携帯電話エリア外)アウトだろう。
登山用にはもう少し手軽な製品が出ているようだ(ヒトココ)。
こちらは最大距離 1 kmと短めだけれど、親機を安全なところに確保しておいて、わずか20gの子機を携帯していくという形で使用できる模様。完全生活防水とのことだけれど、ダイビングで使うなら何かしらの防水ケースが別途必要だろう。距離の制限はあるけれど、水面に浮いてさえしてくれればいいところは汎用性は高そうだ(意識の有無にかかわらず…)。

その他、ビーコンという案もいただいた。

またバイオロギングで海獣にGPS発信機を付けるという話もあるけれど、あれはダイビングの安全管理には使えないのだろうか?

しかし、もし沈んだままで浮いていなかったら、これらの手段はどうにもならないですね。。海底で作業している間に意識を失ったとしたら、中性浮力でなくエアは抜いた状態のはず。

2016年9月2日金曜日

(cowplotパッケージ)研究用にスッキリ簡潔にggplotを描画 & 複数パネル化

最近、Rでのグラフ作成の標準になりつつある気がするggplot2パッケージですが、デフォルトのテーマは研究用としては装飾過剰なので、自分用にアレンジしたテーマを使っている人も多いと思います。でもそのテーマを図示の度に毎回引っ張ってくるのはとても面倒。

それから、複数の異なる種類のグラフを組み合わせて描画するときに、gridExtraパッケージを使うというのがありますが、図示するまではいいけれど、保存する時に行・列数の指定ができないなど、こちらもいろいろ厄介でした。

この両方を一気に解決してくれたのがcowplotパッケージ、研究目的のスッキリしたグラフを作成するのに特化したようなパッケージです。ggplot2を基本にして拡張したようなものなのですが、使い方はかなり簡単です。以下、irisデータを用いて図示してみます。

data(iris) # iris呼び出し

# いずれも要インストール
require(ggplot2)
require(cowplot)

初めにggplot2とcowplotの両方を呼び出しておけば、あとは普段通りにggplotで図示するだけで論文ライクなスッキリしたグラフが描かれます。

# x=Sepal.Length, y=Sepal.Widthの散布図をSpecies毎に色分けして描く
g1 <- ggplot(iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2) 
g1 # 図示



# 保存も普通にggplot2と同様にやるだけ
ggsave("cowsample1.pdf", g1, height=10, width=12, unit="cm")





次は、複数種類のグラフを並べる場合

# 先ほどのg1と新たにg2を用意し、横並びにします。
g2 <- ggplot(iris, aes(x=Petal.Length, y=Petal.Width)) + geom_point(aes(fill=Species), pch=21, size=2)

# 複数グラフを合成するにはplot_grid関数を使用します
g12 <- plot_grid(g1, g2, labels="auto", align="h") 

# labels="auto":小文字で各パネルにラベリング("AUTO"だと大文字)
# labels=c("a","b") のように直接入力することも可能
# align="h":グラフの上下が揃うようにサイズ調整してくれる



# このグラフのように凡例が共通の場合、一方は消してしまいましょうか。
g1n <- g1 + theme(legend.position="none") # ggplotのtheme関数と合わせ技が可能、1パネル目の凡例を消去
g1n2 <- plot_grid(g1n, g2, labels="auto", rel_widths=c(1, 1.4), scale=0.9, vjust=0.2
# rel_widthsで横幅を1:1.4に変更
# それからデフォルトのままだと、パネル間がくっつきすぎてしまうので、scaleパラメータで調整します(1→0.9)。vjustはラベル位置(a, b)の調整






# もう一つ増やして、二段重ねにしてみます(nrow=2で2行になる:gridExtraパッケージのarrageGrob関数と違って、行・列数の指定が可能。grid.arrange関数があるじゃないかと思うかもですが、ファイル保存できないので実用不可)
g3 <- ggplot(iris, aes(x=Species, y=Petal.Width)) + geom_boxplot(aes(fill=Species)) + theme(axis.text.x = element_text(angle=90, vjust=0.5))
# 追加のg3はx軸のラベルが重なってしまうので縦にする↑(themeを使用)
g123 <- plot_grid(g1, g2, g3, labels="auto", nrow=2, scale=0.9, vjust=0.2

# ちゃんとこの状態で保存も可能です
ggsave("cowsample123.pdf", g124, height=15, width=19, unit="cm")





さらに、ggplotらしく要因で分けたグラフと普通のグラフを組み合わせてみます。縦長と普通サイズの組み合わせをしようとしているのですが、この場合、plot_gridを重ねることで実現できます。ラベルをautoのままにすると重なってしまうので注意が必要です。

# g1, g3で縦に連結(横に持ってくるので、ラベルがb, cになるようにする)
# 幅も揃えておきます(align="v")
g13 <- plot_grid(g1, g3, labels=c("b","c"), ncol=1, vjust=0.5, scale=0.9, align="v") 

# 次にfacet_wrapで種毎のグラフ(ncol=1で1列になる)
g4n <- ggplot(iris, aes(x=Sepal.Length, y=Petal.Length)) + geom_point(aes(fill=Species), pch=21, size=2) + facet_wrap(~ Species, ncol=1) + theme(legend.position="none")

# 両者をさらにplot_gridの中に入れます(これによりplot_gridが入れ子になる)。
# 新たな左のグラフのみにラベルをつけるので、c("a","")のように指定しました。
g134 <- plot_grid(g4n, g13, labels=c("a",""), nrow=1, rel_widths=c(1, 1.4), scale=0.9, vjust=0.1)




単体のggplot2と比べて、とても便利。日常的に使いたくなるパッケージです。

2017.05.19追記:
theme関数で外枠を追加しようとすると上手く行かず。代わりにcowplotパッケージのpanel_border関数を使えば良い模様。
+ panel_border(colour=1, size=1) # defaultの線色がなぜかグレーなので黒に変更、sizeで線太さを調節

2016年8月31日水曜日

RのforループをC++で高速化する(Rcppパッケージ)

"時間が掛かるからRでforループを使うな"、"applyファミリーを使え(ほかにsapply, lapplyなど)"とよく言われる。かといって無理に使うと難解なコードになる場合もあるし、せっかく実現しても計測してみたら、むしろforループの方が早かったなんてこともあった(今回の動機)。並列化(複数のCPUコアを用いる)するという手も考えたが、大きなー繋がりのタスクを1回やるだけならよいが、入れ子になっていると呼び出す時間の方が律速になったりもする。

applyファミリーやforeachなどで置換えにくい計算例として、値を毎回更新するようなプログラムがあるだろう。例えば以下のような計算例("0.9* "さえ無ければapply(x, 2, cumsum)で一発なのだが、sapplyとapplyを組み合わせて書いたらforより遅くなった…)。

NN <- 10^6
x <- matrix(0, NN, 10) 
# NN回分の結果を格納するための入れ物(1回あたりは長さ10のベクトル)
system.time( for(r in 2:(nrow(x))) x[r,] <-  0.9*x[r-1,] + rnorm(10) )
# ひとつ前の結果に0.9を掛け、正規乱数を足していくマルコフ連鎖
#  ユーザ   システム       経過
#    5.457      0.363      5.778

手元の環境では6秒弱掛かった。

いろいろ調べて行き着いたのが Rcpp というパッケージ(要インストール)。RからC++("しーぷらぷら"と呼ぶみたい、なのでcppなのだろう)で作成した関数を呼び出せるというもの。Cで作られているというパッケージは最近よく見かけるようになったが、自分でカスタマイズできるのはすごい。もっとも、C++を書けるならばですが。それでも部分的にだけでも記述できればだいぶ速くすることができるはず。目的のforループ作成を達成するのに丸一日掛かりましたが、最終目的の計算に掛かる予想時間が数ヶ月だったことを考えれば大幅な時間短縮になりました。

手順として、まず開発環境のセットアップが必要です。Macの場合はXcode(AppStoreからフリーでダウンロード)とX11(ユーティリティ内にあるX11のアイコンをダブルクリックして立ち上げてライセンスにOKすれば大丈夫なはず)。Windows版はすみませんよく分かりません…コマンドプロンプトでC++をコンパイルする必要があると思います。
次に肝心の、C++で記述したコードが必要です。ただし随所に純粋のC++と異なる書き方を要する箇所があります(私はここでつまづきました)。C++で記述したファイルをRの作業ディレクトリに置いて読み込むと、Rの関数として使用できるようになります。





まずは最小限のコード例:(以下、C++のコードは紫にしておきます)
require(Rcpp) # Rcppパッケージを呼び出す(この行はRに直打ち)

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

double testF(double x) {
    return(x*x);
}

紫で書いた6行分をコードエディタなどで書いて test.cpp というファイル名で作業ディレクトリに保存します( .cpp ってC++ファイルの拡張子なのだと思います)。
C++用コードの初めの3行は決まり文句のようなものです。3行目は用途によって変える場合があるようですが、とりあえず忘れてよさそう。
まずRとの大きな違いはすべての変数型を逐一定義すること。この性質は嫌いじゃないです。むしろ勝手に変数型が変わるのはRが嫌われる点でもありますね。

double は実数の意味、整数なら int 、のように、使う変数の前に半角スペースを置いて書いて定義する必要があります。

testF は、いま定義しようとしている自作関数の名前です。括弧の中で x という変数が与えられた時に { } の処理を行う変数を作成することを宣言しています。

return(x*x); これは x*x の結果を返すことを示しています。もうひとつRと大きな違いとして各行の最後に ;(セミコロン)が入ることです(中括弧 { の後ろ以外)。


次にR上で以下のようにして、この自作関数を呼び出します(初回の呼び出しには数秒掛かるかもしれない)。
sourceCpp("testF.cpp")

すると、testF 関数が使用できるようになります。例えば…
testF(1.41421356) # 一夜一夜に人見頃
結果はちゃんと 2 が返ってくるはず。


ちなみにRの上でC++コードを走らせることも可能です。""で挟んで、code=code.testで指定してやるだけです。
sourceCpp(code="
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

double testF(double x) {
    return(x*x);
} ")

ただし、Rcpp.h以外のファイルを読み込む必要がある場合のやり方が今ひとつ分からず。inlineパッケージのcxxfunctionを使ってできそうなのですが。






次の例では段階を上げて行きます。
rnorm 関数をC++で書いてみましょう。なお、乱数の発生アルゴリズムとしてRの場合はMersenne Twister法が使用されていますが、C++でやろうとすると一工夫必要でした。これについては、詳しい解説とコードも配布されているサイトがありましたのでリンクを貼っておきます。MT.hファイルと内包されるinit_genrand関数、genrand_real3関数はこちらから取得しました。(追記:他にも、匿名さんから頂いたコメントで乱数発生アルゴリズムのまとめサイトがあるとのことです。

r_normal という名前でC++版rnormを作成します。
下準備としてMt.hと名付けたファイルをRの作業ディレクトリに置きます。下記がRcppから読み込む用のC++コードになります(もはや純粋なC++コードではないので、こう呼びます)。

#include <Rcpp.h>
#include "MT.h"

using namespace Rcpp;
// [[Rcpp::export]]

NumericVector r_normal( int NN, double mu, double sigma ) {  
  NumericVector x(NN);
init_genrand((unsigned)time(NULL));

for(int i=0; i<NN; i++) {
x[i] =  mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );
}
return(x);
}

新たに加わった include "MT.h" で上記のファイルを詠みこみます

#include <stdio.h>
#include "MT.h"

次に、この行ですが、これがクセモノです。
NumericVector r_normal( int NN, double mu, double sigma ) {

NumericVector は返り値がベクトルになるように定義しますという意味です(Rcpp特有の表現:cf. Rcpp Quick reference)。そして r_normal という名前の自作関数を作ろうとします。

int NN, double mu, double sigma は、ベクトルの要素数 NN(整数)、平均が mu(実数)、標準偏差が sigma(実数)の値を用いることを宣言しています。

NumericVector x(NN);
次にこれですが、xという、要素数が NN のベクトル(初期値は全部0になっている)を用意することを意味しています(Rcpp特有の定義のようだ)。

init_genrand((unsigned)time(NULL));
これは乱数のタネの指定、現在時刻を秒精度で使用します。

for(int i=0; i<NN; i++) {
C++版のforループ。整数 i を 0からNN未満まで適用しますという意味。0から始めるのが基本である模様(誤解があるかもしれないけれど、とりあえずこのコードではOK)。

x[i] =  mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );

xのi番目の要素に正規乱数を格納するコード。右辺はごちゃごちゃしていますが、genrand_real3が一様乱数を返す関数です。詳しい説明が参照元にあります。
return(x); で x を返り値で返すようにすれば、r_normal関数が使えるようになります。

rnormと速度比較をしてみます。

NN <- 10^7
system.time(rnorm(n, 0, 1))
#   ユーザ   システム       経過
#    0.677      0.006      0.679

system.time(r_normal(n, 0, 1))
#   ユーザ   システム       経過
#    0.363      0.005      0.367

速度は2倍弱といったところ。もう少し期待したかったですが、それでもベクトル化されているRコードよりもさらに高速です。


一方、Rの関数をC++側に渡すことも可能です。Rの関数をC++で計算させることができます。やり方としては、cppFunction関数を用いてrnormのC++版、rnormCppを定義します。
cppFunction( "
NumericVector rnormCpp(int N, double me, double sd) { 
return(rnorm(N, me, sd)); 
}" )

定義の仕方などはだいたい同じですが、rnormが使われていることに注意です。これはC++の関数ではなく、たしかにRのrnorm関数です。
これを通すと、rnormCppという関数が使えるようになります。気になる速度ですが…

system.time(rnormCpp(10^7, 0, 1))
#   ユーザ   システム       経過
#     0.577      0.019      0.594

Rのrnormより少し速くはなりましたが、C++で1から作るより遅いです。それでも手軽に計算速度を向上させることができるので、何かと試したくなる技です。





ではいよいよ、本題のforループを r_markov 関数として記述してみます(同様にファイル保存します)。

#include <Rcpp.h>
#include "MT.h"

using namespace Rcpp;
// [[Rcpp::export]]

NumericMatrix r_markov( int repl, int length, double mu, double sigma ) {

NumericMatrix x(repl, length);
int i,j;

init_genrand((unsigned)time(NULL));

for(i=1; i<repl; i++) {
for(j=0; j<length; j++) {
x(i,j) =  0.9*x(i-1,j) + mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );
}
}
return(x);
}

繰り返しを各行に割り当てるため、int repl、各回のベクトル要素数は int length です。
ひとつ前のr_normal関数との違いは、まず毎回の結果を行列に格納している点です。
NumericMatrixが行列型を示しています。
一行目でもxをrepl行、length列の行列として定義します。NumericMatrix x(repl, length);
毎回の繰り返しでは、ひとつ前の結果に0.9を掛けて正規乱数を足していくので、iのforループは1から始まるようにしています(1行目は0のままで放置、Rコードで除去します)。
もう一つ、分かりにくい違いが添字の書き方です。さっきはx[i] = ...と肩のある括弧を使いましたが、行列の場合はx(i,j)と普通の括弧を使用するようです。C++の解説を見るとx[i][j]とありますが、こうするとエラーになりました。x[i,j]もダメです。



それでは冒頭のR版forループと速度比較してみます。
NN <- 10^6
system.time(r_markov(NN+1, 10, 0, 1)[-1,])
 #  ユーザ   システム       経過
 #   0.439      0.024      0.460

R版では5.778秒だったので10倍以上早いです。C++の速さを思い知らされました。

Rcpp、少しずつ部分的にC++化できるので、練習にも持って来いだと思いました。たぶん純粋C++に移行することはなくハイブリッドでやっていくことになりそうな気がします。

2015年12月11日金曜日

MacOSX10.11 (El Capitan) でWinBUGSを動かす (Wine, R2WinBUGS使用)

OSX10.8以降、長らくアップデートしていなかったので、更新することにしました。もういい加減WinBUGSは卒業したくはあるのですが、依然StanやJAGSでは通らないコードも残っておりまだ手放せずにいます。早くStanが離散パラメータを直接扱えるようになってほしい…

クリーンインストールのEl Capitanに入れることを想定して手順を挙げておきます。まず、X11とHomebrew(パッケージ化されていないアプリを容易にインストールするための補助ツール?)を入れ、Homebrewを用いてWine(非Windows OS上でWin専用アプリケーションを実行する環境)をインストール、Wineのディレクトリ内にWinBUGSをインストールするという流れです。以前はMacPortsを使用していましたが、Homebrewの方がはるかに簡単のようです。以前利用できていたユーザの上の階層がEl Capitanでは使いにくくなったというのも大きな理由です。

なお、インストール作業はターミナルからUNIXコマンドを打ちながらのもの。sudoなどのコマンドは注意深く扱う必要があるようなので、チャレンジする際には慎重に。またXcodeもソフトウェア開発に使うような類のツールなので取り扱い注意です。参考にする際は、この辺りを理解の上、自己責任でお願いします…。

(下記、Rコードは緑、ターミナルのコードは紫にしてみます)

cf. Windows10でWinBUGSを使用するには一手間必要。Program Filesが変更不能になったため、User以下のフォルダにWinBUGSを入れるしか無い。
bugs(..., bugs.directory="C:/Users/ゆーざぁ名/Documents/WinBUGS14/")
のようにして、WinBUGS14.exeファイルを置いている階層を指定してやる必要がある(例は、ドキュメントフォルダ内にWinBUGS14フォルダを置いて、その直下のWinBUGS14.exeを呼び出す場合。"ゆーざぁ名"は御自身の使用しているものに置き換えてください)。

cf. usr/local階層に変更を加えるための認証は、OSのマイナーアップグレードの度に必要になるかもしれない。その場合、ターミナルで以下のコードを打ち込む
sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local



**********************
0)実行環境
・OSX10.11 (El Capitan) 搭載のMac

1)App StoreのApple IDを設定しておく


2)「アプリケーション」→「ユーティリティ」にある「ターミナル」を立ち上げる。
すると、冒頭にこのように出ている。

コンピュータ名:~ ユーザ名$

このドルマーク $ の後にコマンドを打っていく。
なお、インストールに関わるところでパスワードを求められるが、その都度、自分のアカウントのパスワードを入れる。
(以降、パスワードを入れる作業は説明を省略)


3) 下準備、/usr/local/フォルダを作る、ロックをいったん外して操作をするという動作のようなので、推奨されていない動作だということをお忘れなく。

リカバリモード(⌘+R を押しながら起動)で起動し、ターミナルを立ち上げる


4) ターミナルの $ マークの直後に、下記のコードを打ち込む(コピペでOK、以降も同様)

csrutil disable


5) 通常の再起動をする


6) ターミナルに下記を打ち込む(改行されて見えているだろうが、改行無しで打ち込む)

sudo mkdir /usr/local && sudo chflags norestricted /usr/local && sudo chown $(whoami):admin /usr/local && sudo chown -R $(whoami):admin /usr/local


7) 再度、リカバリモードで起動


8) 下記のコードを打ち込む

csrutil enable


9) 今一度、通常の再起動をする


10) homebrewをインストールする

ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

(せいぜい数分でインストールは終了するはず)


11) Wineのインストールを試みる、下記のコードをターミナルに打ち込む

brew install wine

数分ほどでエラーメッセージが出る、メッセージの最後に下記のように書かれている
To continue, you must install Xcode from the App Store, or the CLT by running:
xcode-select --install

先にXcodeをインストールしておけば回避できるのだろうが、このやり方でインストールするほうがむしろ手間が省けるだろう。


12) 上記のメッセージ通り、下記のコードを打ち込む
xcode-select --install

すると、AppStoreからXcodeをインストールしてもよいかと聞かれるのでOKをする。


13) インストールが終わったら、再びWineのインストールを試みる

brew install wine


14) 回線速度にもよるが、良好な環境では30分程度でインストールは終わるだろう。これが終わったら、次のコマンドを打つ。

winecfg

しばらく処理音が聞こえた後、X11からWineの環境設定のようなウインドウが表示される。単に一番下のOKをクリックすればよい。
ターミナルにはエラーメッセージがいくつか出ているが気にしなくてよさそうだ。
Wineのインストールはたったこれだけで終わり…MacPortsの時の面倒を思えば、Wine自体のインストールはずっと容易になった。


15) WinBUGSのインストール
パッチなどを当ててある展開済みのWinBUGSフォルダを用意する。Windows7、Windows8へのインストールも現在はこのやり方で行くしかないことを考えれば、Macでも同様にすればいいだろう。適用済みのWinBUGSも公開されている。

下記のコードで不可視フォルダにあるProgram Filesフォルダを開く

open ~/.wine/drive_c/Program\ Files/

ここへWinBUGSフォルダを入れればインストール完了


16)RからWinBUGSを実行する。下記のような単純なサンプルコードで試してみる。Wineを経由するので、bugs()内にそのためのコードがたくさん必要。

# R2WinBUGSのインストールをお忘れなく
require(R2WinBUGS)
# 真の値は、a=3, b=2, sd=1
X <- c(1:100)
Y <- rnorm(100, mean=(3 + 2*X), sd=1)
data <- list(X=X, Y=Y)
inits <- function() list(a=0, b=0, tau=1) 
parameters <- c("a", "b", "sigma")

model <- function() {
a ~ dnorm(0, 1.0E-6)
b ~ dnorm(0, 1.0E-6)
tau ~ dgamma(1.0E-2, 1.0E-2)
for (i in 1:100) {
Y[i] ~ dnorm(mean[i], tau)
mean[i] <- a + b*X[i] }
sigma <- 1/sqrt(tau)
}
modelpath <- file.path(tempdir(), "model.bug")
write.model(model, modelpath)

mcmc <- bugs(
data=data, inits=inits, parameters=parameters, model.file=modelpath, 
n.chains=3, n.iter=5000, debug=T,
working.directory=NULL, clearWD=T, useWINE=T, newWINE=T,
WINE="/usr/local/bin/wine", WINEPATH="/usr/local/bin/winepath")

print(mcmc) # ちゃんと真の値(a=3, b=2, sigma=1)が推定できたかチェックしよう

# 今回、opt/localではなくusr/localにパスを通すよう変更する必要が出た。以前のWINE、WINEPATHは/opt/local/bin/になっていたが、ここは/usr/local/bin/に変更していることに注意。


17)まだR上で下記のエラーコードが出るが、これはこちら(http://ggorjan.blogspot.jp/2008/10/runnning-r2winbugs-on-mac.html)によると害のないエラーコードらしい。要は推定計算さえ無事に行われていればよいだろう。
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject class {0003000a-0000-0000-c000-000000000046} not registered
err:ole:CoGetClassObject no class object {0003000a-0000-0000-c000-000000000046} could be created for context 0x3
err:ole:CoReleaseMarshalData IMarshal::ReleaseMarshalData failed with error 0x8001011d