ナンクル力学系

学んだ事を書き連ねていこう。

Posts Tagged ‘R

エラーバーと平均を元データと一緒にプロットするRスクリプト

with 2 comments

こんな感じの出力になるやつ:

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) [/sourcecode]

Written by tkf

July 7, 2009 at 1:45 am

Posted in PC

Tagged with

Rで離散時間の力学系を計算するためにC拡張を書いた

leave a comment »

x_0 を与えたときに、 x_{t+1} = f(x_t) を使って 時系列 (x_0, x_1, x_2, \dots ) を計算するって単純なことが したいだけだったのに、上手くRでやる方法が思いつかなかった。

ので、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)&#91;time*i&#93; = REAL(x0)&#91;i&#93;;
  }
  for(t = 1; t < time; t++) {
    PROTECT( x1 = allocVector(REALSXP, dim) );
    for (i = 0; i < dim; ++i){
      REAL(x1)&#91;i&#93; = REAL(xt)&#91;t-1 + time*i&#93;;
    }
    defineVar(install("x"), x1, rho);
    UNPROTECT(1);

    PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
    for(i = 0; i < dim; i++){
      REAL(xt)&#91;t + time*i&#93; = REAL(x2)&#91;i&#93;;
    }
    UNPROTECT(1);
  }
  UNPROTECT(1);
  return(xt);
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>Rの関数(ラッパ?):</h5>

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) )
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>実行例:</h5>

> 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&#91;t-1,&#93;) })
   ユーザ   システム       経過
     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)&#91;time*i&#93; = REAL(x0)&#91;i&#93;;
  }
  for(t = 1; t < time; t++) {
    for (i = 0; i < dim; ++i){
      REAL(x1)&#91;i&#93; = REAL(xt)&#91;t-1 + time*i&#93;;
    }

    PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
    for(i = 0; i < dim; i++){
      REAL(xt)&#91;t + time*i&#93; = REAL(x2)&#91;i&#93;;
    }
    UNPROTECT(1);
  }
  UNPROTECT(2);
  return(xt);
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>実行例:</h5>

> 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&#91;t-1,&#93;) })
   ユーザ   システム       経過
     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

Written by tkf

March 11, 2009 at 7:17 pm

Posted in PC

Tagged with , ,

Rで複数の3D曲線を描く関数を少し直した

leave a comment »

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))) [/sourcecode]

「異様に遅い」は嘘でした><

id:syou6162さんに system.time を教えてもらって 関数を実行する部分だけ時間はかってみたら、全然遅くなかった。

疑問

  • 「system.time(lines3d <- function(略))」 とやってみた
    • 実時間でかなり時間かかる
    • でも、 system.time が表示する実行時間は0秒
    • なぜ? Rは関数の評価に時間かるんだろうか?
  • 「source(‘ファイルへのパス’)」で実行してみると瞬殺
    • RじゃなくてESSの問題とか?

他やること

  • まじめにやるなら、
    • lineごとに色や線種を変える
    • 凡例を出す
    • 外で定義したpmatを使えるようにする
    • 。。。とかやるべきだけど、[あとでやる]
      そんなに複雑なグラフは今は書かないはず。

Written by tkf

March 11, 2009 at 2:15 am

Posted in PC

Tagged with

Rで複数の3D曲線を描く関数を作ってみたけど異様に遅い

leave a comment »

球の三次元プロット (r-help 記事より, 2003.12.27) << グラフィックス参考実例集:自作グラフィックス投稿欄 – RjpWiki を参考にして、作ってみた。

関数と実行例:

# # デフォルトで定義されているっぽい?
# trans3d <- function(x,y,z, pmat) { # tr <- cbind(x,y,z,1) %*% pmat # list(x = tr[,1]/tr[,4], y= tr[,2]/tr[,4]) # } 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 <- list(c( min(linelist[[1]][[1]]), max(linelist[[1]][[1]]) ) , c( min(linelist[[1]][[2]]), max(linelist[[1]][[2]]) ) , c( min(linelist[[1]][[3]]), max(linelist[[1]][[3]]) ) ) for (lin in linelist) { for (i in seq(length(limit))) { limit[[i]] <- c(min(min(lin[[i]]), limit[[i]][1]) , max(max(lin[[i]]), limit[[i]][2]) ) } } # # 上書き(この辺に,xlim=NULLあたりの処理をかくつもり) # 上限と下限の差が0より大きいかチェック for (i in seq(length(limit))) { 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(list(x0,y0,z0), list(x1,y1,z1))) [/sourcecode]

なんで遅い?

  • たぶん、
    • forを使いまくっている
      • Rっぽくないよな
    • listを使いまくっているデータ構造が、そもそもRっぽい処理に向かない?
    • linelistの一つ目の要素を二回しらべている
      • [あとでなおす] 手を抜きました><

Written by tkf

March 10, 2009 at 11:52 pm

Posted in PC

Tagged with

Tsukuba.R#4のLTに使ったRのソースとかやっとUPした。あと+何か。

leave a comment »

発表スライドは Tsukuba.R#4で力学系の分岐現象についてLTしてきた « ナンクル力学系 にある。

なぜかgithubにまとめておくとにした

+何か

今更なんだよ、と思うけれどw

  • Tsukuba.RのWikiにトラバ飛ばすのかなと思ったら俺だけだったので悲しかった。
    • id:syou6162さんがまとめたりしてるけど、もれがないように自主的に送るのが 良いのではと思ったりする
    • なんで俺は同じ記事でふたつも飛ばしてるんだ
  • 聞きたい人がいるような内容ではなかった(気がする)ので、次やるとしたらニーズが ある内容をやりたい
  • ふらふらしている時に良い発表はできません
    • 元気がなかったのでスライド飛ばした(ぉ
  • Tsukuba.Rの運営に関してもらった反応に応えていくよ!! – Seeking for my unique color.
    とそのリンク先を読んで、

    • 「Rというツールを通して他の分野のこと知る」みたいな感じになれば面白いなあ
    • まったく違う分野でも数学的な構造が同じ、ってことは良くある。この「数学的な構造」を 「Rの使い方」に置き換えても成り立つはず
    • これは、Rの使い方(テクニック)を披露しあう場になれば当然成り立つはず
    • だから「背景説明 + テクニックの説明」ってな感じが良いのかな
      • って、良く考えたらだいたいの人は一応そんな感じだったなw
      • 俺のLTはテクニックの説明がなかったので反省
    • 。。。と思った。id:syou6162の考えとかなりズレていることが分かりますねハイ。id:flashingwindさんの意見と若干似てるのかな。

Written by tkf

March 10, 2009 at 9:34 pm

Posted in 研究日誌

Tagged with , ,

Rの練習のために力学系を描いてみた

with 2 comments

Rを使う練習。使い方間違ってる気がしてならない。

odesolveとscatterplot3dが要ります。

Read the rest of this entry »

Written by tkf

February 27, 2009 at 10:27 pm

Posted in PC, 数学

Tagged with ,