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〜４の行程は不要かもしれない＞

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

data(iris) # iris呼び出し

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

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

→（2020.08.29追記）cowplotの仕様が変更になっており、下記のようなcowplotスタイルのグラフにするには、theme_set(theme_cowplot())、を通しておく必要があります。またはggplot()... + theme_cowplot()のように追加してもOKです。

# 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関数と合わせ技が可能、１パネル目の凡例を消去

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で２行になる：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で１列になる）

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で線太さを調節

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

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

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("test.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++で１から作るより遅いです。それでも手軽に計算速度を向上させることができるので、何かと試したくなる技です。

ではいよいよ、本題の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++に移行することはなくハイブリッドでやっていくことになりそうな気がします。