Posts Tagged ‘R’
エラーバーと平均を元データと一緒にプロットするRスクリプト
こんな感じの出力になるやつ:
tapply と factor を使うとかなり楽にかけた．Rすごい!
ソースコード
plot.raw.mean.errorbar <- function(x, y, new=FALSE, col=par("fg"), axes=TRUE,
col.raw=FALSE, col.mean=FALSE, col.eb=FALSE,
xlim=FALSE, ylim=FALSE,
arrows.length=0.1 ){
# calc mean and error bar
fc <- factor(x) ## <<-- ここ！
ux <- as.numeric(levels(fc))
my <- tapply(y, fc, mean) ## <<-- ここ！
sy <- tapply(y, fc, sd) ## <<-- ここ！
ucl <- my + sy
lcl <- my - sy
# plot
if ( col.raw == FALSE ) col.raw <- col
if ( col.mean == FALSE ) col.mean <- col
if ( col.eb == FALSE ) col.eb <- col
if ( xlim == FALSE ) xlim <- range(x)
if ( ylim == FALSE ) ylim <- range(y)
par(new=new)
plot( x, y, xlim=xlim, ylim=ylim, col=col.raw, axes=F, xlab="", ylab="" )
arrows( ux, ucl, ux, lcl,
length=arrows.length, angle=90, code=3, lwd=2, col=col.eb )
par(new=T)
plot( ux, my, xlim=xlim, ylim=ylim, col=col.mean,
axes=axes, type = "l", xlab="", ylab="" )
}
tanh.g <- function(x){
y <- tanh(x)
return( y + (1-y^2) * rnorm(length(x)) )
}
n <- 20
x <- rep(-5:5,n)
y <- tanh.g(x)
plot.raw.mean.errorbar(x,y)
Rで離散時間の力学系を計算するためにC拡張を書いた
ので、Cで拡張作った。
とりあえず動いたver
Cのコード
#include <R.h>
#include <Rinternals.h>
SEXP dtds(SEXP _time, SEXP x0, SEXP fn, SEXP rho)
{
int i, dim = length(x0);
int t, time = INTEGER(_time)[0];
SEXP xt, x1, x2;
PROTECT(xt = allocMatrix(REALSXP, time, dim));
for(i = 0; i < dim; i++){
REAL(xt)[time*i] = REAL(x0)[i];
}
for(t = 1; t < time; t++) {
PROTECT( x1 = allocVector(REALSXP, dim) );
for (i = 0; i < dim; ++i){
REAL(x1)[i] = REAL(xt)[t-1 + time*i];
}
defineVar(install("x"), x1, rho);
UNPROTECT(1);
PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
for(i = 0; i < dim; i++){
REAL(xt)[t + time*i] = REAL(x2)[i];
}
UNPROTECT(1);
}
UNPROTECT(1);
return(xt);
}
Rの関数（ラッパ？）:
dyn.load("dtds.so")
dtds <- function(N, x0, fmap){
.Call("dtds", as.integer(N), as.double(x0), quote(fmap(x)), new.env())
}
fmap <- function(x){
x * (x - 1)
}
rms <- function(x,y){
d <- c(x)-c(y)
sqrt( d %*% d/length(d) )
}
実行例:
> N <- 100000
> x0 <- c(-0.1,0.1)
>
> xt <- matrix(rep(x0, N),N,byrow=T)
> system.time(for (t in 2:N){ xt[t,] <- fmap(xt[t-1,]) })
ユーザ システム 経過
1.936 0.040 1.977
> print(xt[N,])
[1] 0.002237959 -0.002232851
>
> system.time(xs <- dtds(N, x0, fmap))
ユーザ システム 経過
0.42 0.00 0.42
# for使って計算したやつの4倍くらい?
> print(xs[N,])
[1] 0.002237959 -0.002232851
>
> print(rms(xt,xs))
[,1]
[1,] 0
# for使って計算したやつと同じ結果！
あれ？
4倍しかはやくならない？
くやしかったので少し改良したver
- 毎ステップ allocVector しなくて良いはず
- defineVar もボトルネックか？
Cのコード
#include <R.h>
#include <Rinternals.h>
SEXP dtds(SEXP _time, SEXP x0, SEXP fn, SEXP rho)
{
int i, dim = length(x0);
int t, time = INTEGER(_time)[0];
SEXP xt, x1, x2;
PROTECT(xt = allocMatrix(REALSXP, time, dim));
PROTECT( x1 = allocVector(REALSXP, dim) );
defineVar(install("x"), x1, rho);
for(i = 0; i < dim; i++){
REAL(xt)[time*i] = REAL(x0)[i];
}
for(t = 1; t < time; t++) {
for (i = 0; i < dim; ++i){
REAL(x1)[i] = REAL(xt)[t-1 + time*i];
}
PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
for(i = 0; i < dim; i++){
REAL(xt)[t + time*i] = REAL(x2)[i];
}
UNPROTECT(1);
}
UNPROTECT(2);
return(xt);
}
実行例:
> N <- 100000
> x0 <- c(-0.1,0.1)
> #x0 <- c(0.1)
>
> xt <- matrix(rep(x0, N),N,byrow=T)
> system.time(for (t in 2:N){ xt[t,] <- fmap(xt[t-1,]) })
ユーザ システム 経過
1.944 0.040 1.981
> #print(length(xt[1,]))
> #print(length(xt[,1]))
> print(xt[N,])
[1] 0.002237959 -0.002232851
>
> system.time(xs <- dtds(N, x0, fmap))
ユーザ システム 経過
0.388 0.004 0.392
> #print(length(xs[1,]))
> #print(length(xs[,1]))
> print(xs[N,])
[1] 0.002237959 -0.002232851
>
> print(rms(xt,xs))
[,1]
[1,] 0
うーん
- ほとんど変化なし
- やはり、evalがだめな感じがする
感想
PythonのC拡張比べてみて。（昔書いたPythonのC拡張はこんな感じ。 この例と同じ離散時間の系の計算だけど、fmapにあたる写像は決め打ち。）
- PythonでC拡張書くより、全然書く量が少ないのは良い
- ドキュメントがなさすぎるw
- Writing R Extensions だけ
- リファレンスが無いので、関数の役目が謎なままつかうことに
- 関数とかマクロとかフリーダム
- PythonのCコードは必ず「Py」から始まっていたり、マクロならそれ以下が大文字 になっていたりと安心する
- まー、Cだし、コンパイラが怒るので問題なし？
- 何回でも dyn.load(“dtds.so”) と dyn.unload(“dtds.so”) できるのはいいね！
- Pythonだと、一度読み込んだCのモジュールは reload しても変化なしなので，Python自体を再起動する必要があった
- デバッグ楽
- Cコードレベルのミス（？）に気づくRいいね
- 何度、 「 エラー： VECTOR_ELT() can only be applied to a ‘list’, not a ‘double’」 と怒られたことか。。。
- 慣れたらサクサクかけそう！
- 慣れる ＝ Rのソース読めるくらいか？w
Rで複数の3D曲線を描く関数を少し直した
Rで複数の3D曲線を描く関数を作ってみたけど異様に遅い « ナンクル力学系 で作った関数を少し変更。
関数と実行例:
lines3d <-
function(linelist,
theta=25, phi=30, expand=1,
xlab="X", ylab="Y", zlab="Z",
ticktype="detailed",
xlim=NULL, ylim=NULL, zlim=NULL ){
# x,y,z の上限と下限を求める
limit <- matrix(c(range(linelist[[1]][,1]) ,
range(linelist[[1]][,2]) ,
range(linelist[[1]][,3]) ),
c(3,2), byrow=T)
for (lin in linelist) {
for (i in seq(3)) {
limit[i,] <- range( lin[,i], limit[i,] )
}
}
# limitが指定されていれば上書き
if ( !is.null(xlim) ){ limit[1,] <- xlim }
if ( !is.null(ylim) ){ limit[2,] <- ylim }
if ( !is.null(zlim) ){ limit[3,] <- zlim }
# 上限と下限の差が非零かチェック
for (i in seq(3)) {
if ( limit[i,1] == limit[i,2] ){
limit[i,1] <- limit[i,1] - 1
limit[i,2] <- limit[i,2] + 1
}
}
pmat <-
persp(0:1, 0:1, matrix(,2,2),
xlim=limit[1,], ylim=limit[2,], zlim=limit[3,],
theta=theta, phi=phi, expand=expand,
xlab=xlab, ylab=ylab, zlab=zlab,
ticktype=ticktype )
for (lin in linelist) {
lines(trans3d(lin[,1], lin[,2], lin[,3], pmat))
}
return(pmat)
}
x0 <- seq(-1,1,length=1000)
y0 <- cos(x0*6*pi)
z0 <- rep(0,length(x0))
x1 <- seq(-1,1,length=200)
y1 <- rep(0,length(x1))
z1 <- sin(x1*2*pi)
lines3d(list(cbind(x0,y0,z0),
cbind(x1,y1,z1),
cbind(x1,y1+1,z1*0.2),
cbind(x1,y1-1,z1*0.5)))
「異様に遅い」は嘘でした＞＜
id:syou6162さんに system.time を教えてもらって 関数を実行する部分だけ時間はかってみたら、全然遅くなかった。
疑問
- 「system.time(lines3d <- function（略）)」 とやってみた
- 実時間でかなり時間かかる
- でも、 system.time が表示する実行時間は0秒
- なぜ？ Rは関数の評価に時間かるんだろうか？
- 「source(‘ファイルへのパス’)」で実行してみると瞬殺
- RじゃなくてESSの問題とか？
他やること
- まじめにやるなら、
- lineごとに色や線種を変える
- 凡例を出す
- 外で定義したpmatを使えるようにする
- 。。。とかやるべきだけど、[あとでやる]
そんなに複雑なグラフは今は書かないはず。
