海中での潜水調査・作業をよくやる人の中には、ヒヤッとした経験のある人も多いだろう。私自身もかつてBC（浮力調節）ジャケットの弁のトラブルで強制浮上になってしまったことがあり、それ以降はとくに気になるようになってきた。
先週、知っている人の中で重大な潜水事故が起こってしまった。まったくもって今更ではあるけれど、ここでは海中または浮上後の海上でトラブルが起こった際に無事に帰還するための装備について考えてみたい。ご意見も歓迎します。
まず、流されてしまったか何かで、岸や船から遠い地点に浮上してしまった場合、BCに付属している笛を吹く程度では、よほど凪でもなければ気づいてもらえないだろう。他に入手しやすいものとしては、海面上に長く伸びるフロート（レギュレータからのエアを入れて膨らませられる）があるだろう。私も念のため持っているけれど、しかし実際これで十分に目立つのだろうか。
他の手段…例えばスマホを完全防水で携帯して、浮上したら電話するというのはどうだろうか。場合によってはアカウントを利用して、PCからスマホを探すことも可能な気がするのだけれど。
スマホに近いものとして、ココセコムというのがあるそうだ。子供やお年寄りに持たせておいて携帯電話網でマップ上で追跡が可能、料金的にもスマホ利用より安そうです。
無線を使うという手段は電波法上かなり難しいようだが、合法的にダイビングで使用できる製品もあるようだ（トランスポンダ SEAKER_1）
船からのダイビングという形態を前提としているようだけれど、通常使用（自分の母船に連絡）と、緊急使用（母船とはぐれている場合に無線を受けられるすべての船に緊急発信が可能）の２通りの発信が可能というのが優れている。距離的には20 kmまで届くようなので頼もしい。
（2016.12.16 追記）意識があって海上で操作することが前提ですが、その条件ではかなり有効な製品が出たようです（個人用救難信号発信機 PLB）。精度100mで全世界対応。機器そのものは今年３月に日本で認可、防水ケースが先月末に発売された模様。無線局として申請・許可を得る必要があるなど取り扱いは注意ではあります。
では、意識がないけれど海面には浮かんでいる場合どうするか。上記の携帯を利用した方法が可能であったら使えそうではあるけれど、岸から遠くへ流されてしまっていたとしたら（携帯電話エリア外）アウトだろう。
登山用にはもう少し手軽な製品が出ているようだ（ヒトココ）。
こちらは最大距離 1 kmと短めだけれど、親機を安全なところに確保しておいて、わずか20gの子機を携帯していくという形で使用できる模様。完全生活防水とのことだけれど、ダイビングで使うなら何かしらの防水ケースが別途必要だろう。距離の制限はあるけれど、水面に浮いてさえしてくれればいいところは汎用性は高そうだ（意識の有無にかかわらず…）。
その他、ビーコンという案もいただいた。
またバイオロギングで海獣にGPS発信機を付けるという話もあるけれど、あれはダイビングの安全管理には使えないのだろうか？
しかし、もし沈んだままで浮いていなかったら、これらの手段はどうにもならないですね。。海底で作業している間に意識を失ったとしたら、中性浮力でなくエアは抜いた状態のはず。
2016年11月20日日曜日
2016年9月2日金曜日
（cowplotパッケージ）研究用にスッキリ簡潔にggplotを描画 & 複数パネル化
最近、Rでのグラフ作成の標準になりつつある気がするggplot2パッケージですが、デフォルトのテーマは研究用としては装飾過剰なので、自分用にアレンジしたテーマを使っている人も多いと思います。でもそのテーマを図示の度に毎回引っ張ってくるのはとても面倒。
それから、複数の異なる種類のグラフを組み合わせて描画するときに、gridExtraパッケージを使うというのがありますが、図示するまではいいけれど、保存する時に行・列数の指定ができないなど、こちらもいろいろ厄介でした。
この両方を一気に解決してくれたのがcowplotパッケージ、研究目的のスッキリしたグラフを作成するのに特化したようなパッケージです。ggplot2を基本にして拡張したようなものなのですが、使い方はかなり簡単です。以下、irisデータを用いて図示してみます。
data(iris) # iris呼び出し
# いずれも要インストール
require(ggplot2)
require(cowplot)
初めにggplot2とcowplotの両方を呼び出しておけば、あとは普段通りにggplotで図示するだけで論文ライクなスッキリしたグラフが描かれます。
それから、複数の異なる種類のグラフを組み合わせて描画するときに、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です。
# 複数グラフを合成するにはplot_grid関数を使用します
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)
# 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)
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を使用）
# 追加の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")
# 幅も揃えておきます（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で１列になる）
# 両者をさらに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で線太さを調節
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コアを用いる）するという手も考えたが、大きなー繋がりのタスクを１回やるだけならよいが、入れ子になっていると呼び出す時間の方が律速になったりもする。
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++) {
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++に移行することはなくハイブリッドでやっていくことになりそうな気がします。
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]]
return(x*x);
}
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 <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() );
return(x); で x を返り値で返すようにすれば、r_normal関数が使えるようになります。
rnormと速度比較をしてみます。
NN <- 10^7
system.time(rnorm(n, 0, 1))
# ユーザ システム 経過
# 0.677 0.006 0.679
# ユーザ システム 経過
# 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++に移行することはなくハイブリッドでやっていくことになりそうな気がします。
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、以降も同様）
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をクリックすればよい。
ターミナルにはエラーメッセージがいくつか出ているが気にしなくてよさそうだ。
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
クリーンインストールの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自体のインストールはずっと容易になった。
パッチなどを当ててある展開済みの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
2015年1月29日木曜日
RでGIS その 1：シェープファイル操作、図示
RでのGIS操作、いずれまとめようと思いつつ放ったらかしてました。だれにでも有用そうなものから少しずつアップしていく予定です。
基本の関数の備忘録、とくにshapefileの読み込みと書き出しの関数が長くて忘れてしまいがちです。。
基本の関数の備忘録、とくにshapefileの読み込みと書き出しの関数が長くて忘れてしまいがちです。。
require(maptools) # shapefileの読み込みなどに用いるパッケージ（要インストール）
shape <- readShapeSpatial(file.choose()) # 読み込み & .shpファイルをメニューで選択
plot(shape) # 図示もできます
zoom(shape) # 自分で選んだ範囲を拡大する場合。
# このコマンドを打った後に、図のウィンドウ上で対角線の端と端をクリックで選択するという、Rらしからぬ操作法でズームします
plot(shape2, add=T) # 他のファイルshape2を重ねて図示したい場合
str(shape, 5) # shapeの中身を眺める場合（"5"くらいに制限しておかないとコンソールが溢れて大変なことになる）
shape@data # shapeファイルのデータを取り出す場合（@dataの中身はデータフレーム）
# 手持ちのデータフレームDataからshapefileを作る場合（LongitudeとLatitudeの列を含むデータとします）
# GPSデータはWGS84（133.33333のような表記）にするのが原則です
require(sp) # maptoolsを使用していれば、新たに呼び出さなくてよいはず
coordinates(Data) <- c("Longitude", "Latitude") # このようにGPS列を指定すると空間データ化する
# 変な感じがするかもしれないですが、x, yの順番なのでLongitudeを先に書きます
bbox(Data) # cf. これをやるとデータの四隅（最少・最大）が分かります
#もしデータがグリッド状に揃っている場合はグリッドに変換することができる
gridded(Data) <- TRUE # そうでない場合はエラーになるはず
# データの書き出し
writeSpatialShape(Data, "ファイル名.shp") # readの場合とはShapeとSpatialの順番が逆！
2014年6月10日火曜日
累積ロジットとGLM二項分布の比較・再＆続
（うっかり、同様の検証記事を消してしまったので、ついでにアップデートします）
前回に引き続き、段階的なカテゴリーデータのモデリングに用いられる累積ロジット（cumulative logit）、
例えば、悪い、ふつう、よい、のようなデータを関連しそうな要因で解析するような場合に用いる。
しかし、悪い＝0、ふつう＝1、よい＝2、のように数値化してしまえば、二項分布型のGLMではダメなのだろうか？たぶん、間隔のいびつなデータでは累積ロジットの方が適しているのだろうが。
テストデータを用いて、両者の推定と推定値の求め方を比較する。
# まず解析用データの設定
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 3 # 最大3、つまり 0, 1, 2, 3 の値を取りうる
X <- rep(c(1:10), each=10)
Y1 <- rbinom(100, N, prob=logistic(-5 + 0.8*X)) # logisticの中身を基にして、二項乱数を発生
Y2 <- factor(Y1, ordered=T) # 累積ロジット用に、ランク化した応答変数を用意
D <- data.frame(X, Y1, Y2)
require(VGAM)
M1 <- glm(cbind(Y1, N-Y1) ~ X, family=binomial, data=D) # 二項分布のGLM
M2 <- clm(Y2 ~ X, data=D) # 累積ロジット
M3 <- vglm(Y2 ~ X, family=cumulative(parallel=T), data=D)
# 比較用にvglm版、parallel=Tで切片のみ複数になる、Fにすると回帰係数もランクごとに推定（切片のみランク毎の場合、比例オッズモデルともいうようだ）
# clmでは、
M2@fitted.values で同様にランク毎の確率が得られる。
#### ついでにベイズ版でも計算
前回に引き続き、段階的なカテゴリーデータのモデリングに用いられる累積ロジット（cumulative logit）、
例えば、悪い、ふつう、よい、のようなデータを関連しそうな要因で解析するような場合に用いる。
しかし、悪い＝0、ふつう＝1、よい＝2、のように数値化してしまえば、二項分布型のGLMではダメなのだろうか？たぶん、間隔のいびつなデータでは累積ロジットの方が適しているのだろうが。
テストデータを用いて、両者の推定と推定値の求め方を比較する。
# まず解析用データの設定
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 3 # 最大3、つまり 0, 1, 2, 3 の値を取りうる
X <- rep(c(1:10), each=10)
Y1 <- rbinom(100, N, prob=logistic(-5 + 0.8*X)) # logisticの中身を基にして、二項乱数を発生
Y2 <- factor(Y1, ordered=T) # 累積ロジット用に、ランク化した応答変数を用意
D <- data.frame(X, Y1, Y2)
# 解析の実行
require(ordinal)require(VGAM)
M1 <- glm(cbind(Y1, N-Y1) ~ X, family=binomial, data=D) # 二項分布のGLM
M2 <- clm(Y2 ~ X, data=D) # 累積ロジット
M3 <- vglm(Y2 ~ X, family=cumulative(parallel=T), data=D)
# 比較用にvglm版、parallel=Tで切片のみ複数になる、Fにすると回帰係数もランクごとに推定（切片のみランク毎の場合、比例オッズモデルともいうようだ）
### GLM
summary(M1)
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -5.9324 0.6644 -8.929 <2e-16 ***
# X 0.9241 0.1011 9.141 <2e-16 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# (Dispersion parameter for binomial family taken to be 1)
#
# Null deviance: 287.723 on 99 degrees of freedom
# Residual deviance: 84.861 on 98 degrees of freedom
# AIC: 139.14
### 累積ロジット（clm）
summary(M2)
# link threshold nobs logLik AIC niter max.grad cond.H
# logit flexible 100 -67.06 142.12 5(0) 2.20e-08 2.6e+03
#
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# X 1.2295 0.1659 7.41 1.26e-13 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Threshold coefficients:
# Estimate Std. Error z value
# 0|1 5.955 0.913 6.522
# 1|2 8.082 1.130 7.153
# 2|3 9.716 1.283 7.575
### 累積ロジット（vglm）# 回帰係数の正負が逆になる
summary(M3)
# Coefficients:
# Estimate Std. Error z value
# (Intercept):1 5.9548 0.92423 6.4430
# (Intercept):2 8.0821 1.14814 7.0393
# (Intercept):3 9.7158 1.31086 7.4117
# X -1.2295 0.16832 -7.3046
#
# Number of linear predictors: 3
#
# Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3])
#
# Dispersion Parameter for cumulative family: 1
#
# Residual deviance: 134.1208 on 296 degrees of freedom # 過小分散気味、clmには無い情報！
#
# Log-likelihood: -67.06038 on 296 degrees of freedom
### 推定値を得るには少し手間がかかる
pre1 <- round(3*fitted(M1)) # GLM：試行回数*確率、を整数値に丸める
pre2 <- as.numeric(predict(M2,type="class")$fit) - 1
# 累積ロジット（clm）：predict関数で推定ランクを求め、
# ラベルに合うように 1 を引く（ランクは1,2,3,4、元の値は0,1,2,3なので）
pre3 <- apply(fitted(M3), 1, which.max) - 1
# 累積ロジット（vglm）：何番目のランクの確率が最大かを求め、
# ラベルに合うように 1 を引く
cbind(X, Y1, pre1, pre2, pre3)
# 推定は完全ではないが、GLMと累積ロジットの推定結果は近いものになった。
X Y1 pre1 pre2 pre3
1 1 0 0 0 0
2 1 0 0 0 0
3 1 0 0 0 0
4 1 0 0 0 0
5 1 0 0 0 0
（中略）
51 6 1 1 1 1
52 6 1 1 1 1
53 6 2 1 1 1
54 6 2 1 1 1
55 6 2 1 1 1
56 6 1 1 1 1
57 6 2 1 1 1
58 6 1 1 1 1
59 6 1 1 1 1
60 6 0 1 1 1
61 7 3 2 2 2
62 7 1 2 2 2
63 7 1 2 2 2
64 7 3 2 2 2
65 7 3 2 2 2
66 7 1 2 2 2
67 7 3 2 2 2
68 7 2 2 2 2
69 7 1 2 2 2
70 7 3 2 2 2
71 8 3 2 3 3
72 8 2 2 3 3
73 8 3 2 3 3
74 8 3 2 3 3
75 8 2 2 3 3
76 8 1 2 3 3
77 8 1 2 3 3
78 8 3 2 3 3
79 8 3 2 3 3
80 8 3 2 3 3
81 9 2 3 3 3
82 9 1 3 3 3
83 9 3 3 3 3
（後略）
### では、推定確率の曲線を描くにはどうしたらよいか？
# これがけっこう面倒くさい。まずはvglmを用いてチェックする。
# vglmをfitted()すると、ランク毎の確率が出てくる
head(fitted(M3)) # head()で頭出しする
0 1 2 3
1 0.991209875 0.007734525 0.0008493615 0.0002062383
2 0.991209875 0.007734525 0.0008493615 0.0002062383
3 0.991209875 0.007734525 0.0008493615 0.0002062383
4 0.991209875 0.007734525 0.0008493615 0.0002062383
5 0.991209875 0.007734525 0.0008493615 0.0002062383
6 0.991209875 0.007734525 0.0008493615 0.0002062383
# clmでは、
M2@fitted.values で同様にランク毎の確率が得られる。
# まず、この出来合いの推定値を図示してみる
plot(X, fitted(M3)[,1], xlim=c(0,10), ylim=c(0,1), col="lightblue", ylab="probability") # 0
points(X, fitted(M3)[,2], col="turquoise") # 1
points(X, fitted(M3)[,3], col="royalblue") # 2
points(X, fitted(M3)[,4], col="darkblue") # 3
# モデルの構造を考えると推定確率曲線は次の計算でいいはず
#（coef(M3)[4]は回帰係数だが、正負が逆なので - を付ける）
# 累積ロジットの名前通り、累積で表されている
curve(1 - logistic(-coef(M3)[4]*x - coef(M3)[1]),
add=T, lwd=2, col="lightblue") # ランク0: 1 - ランク1の累積確率
add=T, lwd=2, col="lightblue") # ランク0: 1 - ランク1の累積確率
curve(logistic(-coef(M3)[4]*x - coef(M3)[1]) - logistic(-coef(M3)[4]*x - coef(M3)[2]),
add=T, lwd=2, col="turquoise") # ランク1: ランク1の累積確率 - ランク2の累積確率
curve(logistic(-coef(M3)[4]*x - coef(M3)[2]) - logistic(-coef(M3)[4]*x - coef(M3)[3]),
add=T, lwd=2, col="royalblue") # ランク2: ランク2の累積確率 - ランク3の確率
curve(logistic(-coef(M3)[4]*x - coef(M3)[3]),
add=T, lwd=2, col="darkblue") # ランク3: ランク3の確率（もはや累積でない）
curve(logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="tomato", lwd=2) # 比較用にGLMも
add=T, lwd=2, col="darkblue") # ランク3: ランク3の確率（もはや累積でない）
curve(logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="tomato", lwd=2) # 比較用にGLMも
# 色が薄い〜濃いにかけて、それぞれ、0, 1, 2, 3 になる確率、赤はGLMの確率
# ちゃんと関数による推定プロットと、自前で作った推定曲線が一致した。計算の仕方はこれでよさそうだ。
# ちなみにclmの場合はこう計算する（回帰係数の前の - が不要）
curve(1 - logistic(coef(M2)[4]*x - coef(M2)[1]), add=T, lwd=2, col="lightblue") # 0
curve(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2]),
add=T, lwd=2, col="turquoise") # 1
curve(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3]),
add=T, lwd=2, col="royalblue") # 2
curve(logistic(coef(M2)[4]*x - coef(M2)[3]), add=T, lwd=2, col="darkblue") # 3
### 次に、推定値の曲線を図示してみる（上のは推定確率でした）
# 比較対照用にGLMの曲線と比較する
# 元データ
plot(Y1 ~ X, data=D, xlim=c(0,10), ylim=c(0,3))
# GLM
curve(3*logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="red", lwd=2)
# 累積ロジット（全部足し合わせるような累積構造になる）
curve(
1*(1 - logistic(coef(M2)[4]*x - coef(M2)[1])) # ランク0
+ 4*(logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク3
1*(1 - logistic(coef(M2)[4]*x - coef(M2)[1])) # ランク0
+ 2*(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2])) # ランク1
+ 3*(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク2
- 1, # 0, 1, 2, 3なので、1,2,3,4から1を引く
add=T, col="blue", lwd=2)
# 赤：GLM、青：累積ロジット。この単純なケースではほとんど推定値は変わらない。
#### ついでにベイズ版でも計算
# 下準備、データを累積に変換。1以上、2以上、3以上にまとめる
rank <- matrix(0, nrow=100, ncol=3) # 3はランクの数-1、100はいわゆるN数
for (r in 1:3) rank[Y1 >= r, r] <- 1
# GLM用の用意していた数値データY1を利用、ランク1~3以上の場合に各列に1を入れる
model <- function() {
for (j in 1:3) { # ランクの数-1
for (i in 1:100) { # いわゆるN数
rank[i,j] ~ dbern(p[i,j]) # それぞれベルヌーイ分布
logit(p[i,j]) <- alpha[j] + beta*X[i] # 切片だけランク毎になってる
}
alpha[j] ~ dnorm(0, 1E-6) # 切片をランク毎に推定
}
beta ~ dnorm(0, 1E-6) # 回帰係数は共通
}
data <- list(X=X, rank=rank)
parm <- list(beta=0, alpha=rnorm(3))
source("WBUGS.R") # 自作ラッパー関数
out <- wbugs(data, parm, model)
# 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10
# n.sims = 3000 iterations saved
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# beta 1.328 0.166 1.022 1.217 1.318 1.434 1.680 1.001 2200
# alpha[1] -6.377 0.889 -8.217 -6.960 -6.327 -5.770 -4.703 1.002 1100
# alpha[2] -8.645 1.142 -11.040 -9.370 -8.563 -7.851 -6.528 1.002 1800
# alpha[3] -10.518 1.336 -13.331 -11.373 -10.420 -9.585 -8.139 1.001 3000
# deviance 154.787 2.934 151.200 152.600 154.050 156.200 162.000 1.002 1500
# DIC info (using the rule, pD = Dbar-Dhat)
# pD = 4.0 and DIC = 158.8
rank <- matrix(0, nrow=100, ncol=3) # 3はランクの数-1、100はいわゆるN数
for (r in 1:3) rank[Y1 >= r, r] <- 1
# GLM用の用意していた数値データY1を利用、ランク1~3以上の場合に各列に1を入れる
model <- function() {
for (j in 1:3) { # ランクの数-1
for (i in 1:100) { # いわゆるN数
rank[i,j] ~ dbern(p[i,j]) # それぞれベルヌーイ分布
logit(p[i,j]) <- alpha[j] + beta*X[i] # 切片だけランク毎になってる
}
alpha[j] ~ dnorm(0, 1E-6) # 切片をランク毎に推定
}
beta ~ dnorm(0, 1E-6) # 回帰係数は共通
}
data <- list(X=X, rank=rank)
parm <- list(beta=0, alpha=rnorm(3))
source("WBUGS.R") # 自作ラッパー関数
out <- wbugs(data, parm, model)
# 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10
# n.sims = 3000 iterations saved
# mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
# beta 1.328 0.166 1.022 1.217 1.318 1.434 1.680 1.001 2200
# alpha[1] -6.377 0.889 -8.217 -6.960 -6.327 -5.770 -4.703 1.002 1100
# alpha[2] -8.645 1.142 -11.040 -9.370 -8.563 -7.851 -6.528 1.002 1800
# alpha[3] -10.518 1.336 -13.331 -11.373 -10.420 -9.585 -8.139 1.001 3000
# deviance 154.787 2.934 151.200 152.600 154.050 156.200 162.000 1.002 1500
# DIC info (using the rule, pD = Dbar-Dhat)
# pD = 4.0 and DIC = 158.8
# BUGS版では、切片が正負が逆になった
# きちんとコードを書いた結果がこれなので、符号が逆なのはvglmやclmなのだが…じつにややこしい。
2014年6月9日月曜日
累積ロジットの汎用Rパッケージ {ordinal}
累積ロジット（cumulative logit model）を使う際に、今ひとつ使い勝手のいいRパッケージがないのが気になっていた。
（cf. 累積ロジット：よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ）
例えば、vglmではランダム効果が入れられない。
mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない（手計算すれば良い話ではあるが…）。
最近、ordinalというパッケージを見つけた。１つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。
ちなみに、GLM版にstepAICは使えたが、dredgeはダメだった。
GLMM版では、モデル選択関数はdrop1が使用できた。
require(ordinal)
example(ordinal) # パッケージで用意されている例を出力
# まずはGLM版
ordinl> ## A simple cumulative link model:
ordinl> fm1 <- clm(rating ~ contact + temp, data=wine)
ordinl> summary(fm1)
formula: rating ~ contact + temp
data: wine
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|) # 回帰係数部分（この例ではカテゴリカルだが）
contactyes 1.5278 0.4766 3.205 0.00135 **
tempwarm 2.5031 0.5287 4.735 2.19e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Threshold coefficients: # 各ランクの切片
Estimate Std. Error z value
1|2 -1.3444 0.5171 -2.600
2|3 1.2508 0.4379 2.857
3|4 3.4669 0.5978 5.800
4|5 5.0064 0.7309 6.850
# 比較対照用にvglmでの結果
require(VGAM)
vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)
Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.3444 0.50850 -2.6438
(Intercept):2 1.2508 0.43908 2.8487
(Intercept):3 3.4669 0.59711 5.8061
(Intercept):4 5.0064 0.72906 6.8669
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意
tempwarm -2.5031 0.53199 -4.7052
# ordinalに戻って、こちらはGLMM版
ordinl> ## A simple cumulative link mixed model:
ordinl> fmm1 <- clmm(rating ~ contact + temp + (1|judge), data=wine) # glmer()同様の構造
ordinl> summary(fmm1)
Cumulative Link Mixed Model fitted with the Laplace approximation
formula: rating ~ contact + temp + (1 | judge)
data: wine
link threshold nobs logLik AIC niter max.grad cond.H # ちゃんとAICも算出する
logit flexible 72 -81.57 177.13 331(996) 1.05e-05 2.8e+01
Random effects: # ランダム効果
Groups Name Variance Std.Dev.
judge (Intercept) 1.279 1.131
Number of groups: judge 9
Coefficients:
Estimate Std. Error z value Pr(>|z|)
contactyes 1.8349 0.5125 3.580 0.000344 ***
tempwarm 3.0630 0.5954 5.145 2.68e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Threshold coefficients:
Estimate Std. Error z value
1|2 -1.6237 0.6824 -2.379
2|3 1.5134 0.6038 2.507
3|4 4.2285 0.8090 5.227
4|5 6.0888 0.9725 6.261
（次の記事に続く）
（cf. 累積ロジット：よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ）
例えば、vglmではランダム効果が入れられない。
mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない（手計算すれば良い話ではあるが…）。
最近、ordinalというパッケージを見つけた。１つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。
ちなみに、GLM版にstepAICは使えたが、dredgeはダメだった。
GLMM版では、モデル選択関数はdrop1が使用できた。
require(ordinal)
example(ordinal) # パッケージで用意されている例を出力
# まずはGLM版
ordinl> ## A simple cumulative link model:
ordinl> fm1 <- clm(rating ~ contact + temp, data=wine)
ordinl> summary(fm1)
formula: rating ~ contact + temp
data: wine
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|) # 回帰係数部分（この例ではカテゴリカルだが）
contactyes 1.5278 0.4766 3.205 0.00135 **
tempwarm 2.5031 0.5287 4.735 2.19e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Threshold coefficients: # 各ランクの切片
Estimate Std. Error z value
1|2 -1.3444 0.5171 -2.600
2|3 1.2508 0.4379 2.857
3|4 3.4669 0.5978 5.800
4|5 5.0064 0.7309 6.850
# 比較対照用にvglmでの結果
require(VGAM)
vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)
Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.3444 0.50850 -2.6438
(Intercept):2 1.2508 0.43908 2.8487
(Intercept):3 3.4669 0.59711 5.8061
(Intercept):4 5.0064 0.72906 6.8669
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意
tempwarm -2.5031 0.53199 -4.7052
# ordinalに戻って、こちらはGLMM版
ordinl> ## A simple cumulative link mixed model:
ordinl> fmm1 <- clmm(rating ~ contact + temp + (1|judge), data=wine) # glmer()同様の構造
ordinl> summary(fmm1)
Cumulative Link Mixed Model fitted with the Laplace approximation
formula: rating ~ contact + temp + (1 | judge)
data: wine
link threshold nobs logLik AIC niter max.grad cond.H # ちゃんとAICも算出する
logit flexible 72 -81.57 177.13 331(996) 1.05e-05 2.8e+01
Random effects: # ランダム効果
Groups Name Variance Std.Dev.
judge (Intercept) 1.279 1.131
Number of groups: judge 9
Coefficients:
Estimate Std. Error z value Pr(>|z|)
contactyes 1.8349 0.5125 3.580 0.000344 ***
tempwarm 3.0630 0.5954 5.145 2.68e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Threshold coefficients:
Estimate Std. Error z value
1|2 -1.6237 0.6824 -2.379
2|3 1.5134 0.6038 2.507
3|4 4.2285 0.8090 5.227
4|5 6.0888 0.9725 6.261
（次の記事に続く）
登録:
投稿 (Atom)