遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

RとC++:乱数の生成２

R Rcpp inline C++ 乱数
  • こちらには、"3.3 Can I use code from the Rmath header and library with Rcpp ? Can I call functions defined in the Rmath header file and the standalone math library for R–as for example the random number generators?" という質問がある
  • Rで使っている(もともとはCのソースを使っている)関数の呼び出しに関すること
  • RNGscope()っていうのはRandom Number Generator Scopeのはずだが…何かしら、シードをコントロールする仕組みのようだ
  • たとえば、次のように２通りを作る
fx <- cxxfunction(signature(),'RNGScope();return rnorm(5, 0, 100);',plugin="Rcpp")
fx2 <- cxxfunction(signature(),'return rnorm(5, 0, 100);',plugin="Rcpp")
  • RNGScope z;というのもあるらしい…(こちらを参照)
  • いずれにしろ、コンパイルするときにシードをセットしてしまうもののようなのだが…
fx3 <- cxxfunction(signature(),'RNGScope z;return rnorm(5, 0, 100);',plugin="Rcpp")
  • 結果を見てみよう。RNGScope()がなんなのかの確認が必要だが…
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> tmp<-runif(10)
> fx()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556
> fx()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556

> fx2()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx2()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556
> fx2()
[1] -14.58635 127.54590 157.75196  68.11208 254.32503
> fx2()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx2()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> fx()
[1]   90.12386  -89.43372  -81.24207  -42.35224 -104.18918
> tmp<-runif(10)
> fx()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556
> fx()
[1]  188.14939 -111.80350  139.42556  -28.01398   19.08556
> fx2()
[1] -14.58635 127.54590 157.75196  68.11208 254.32503
> fx2()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> tmp<-runif(10)
> fx2()
[1] -14.58635 127.54590 157.75196  68.11208 254.32503
> fx2()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx2()
[1]  -30.88188 -100.02309  160.37386  -12.11866   84.19376
> fx2()
[1] -129.85674  -87.10575  -96.05648   52.80472  -66.53405
> tmp
 [1] 0.9700476 0.4158934 0.1317760 0.3208837 0.9183798 0.5305914
 [7] 0.3896851 0.5477884 0.5756806 0.1227024
> runif(10)
 [1] 0.4420146 0.6912912 0.8989266 0.3812118 0.9426620 0.7244207
 [7] 0.7521025 0.2034892 0.9945087 0.7337561
> fx2()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx()
[1] 155.97619 259.53612  61.05704 -52.29941  51.76066
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx2()
[1]  -30.88188 -100.02309  160.37386  -12.11866   84.19376
> fx2()
[1] -129.85674  -87.10575  -96.05648   52.80472  -66.53405
> runif(10)
 [1] 0.9405919 0.1702227 0.9952754 0.8375223 0.7292580 0.6240177
 [7] 0.3004892 0.9802543 0.6976336 0.7320461
> fx2()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx3 <- cxxfunction(signature(),'RNGScope z;return rnorm(5, 0, 100);',plugin="Rcpp")
> fx3()
[1]  -96.65321   20.55194  103.07353   91.79665 -207.66010
> fx3()
[1]  -30.88188 -100.02309  160.37386  -12.11866   84.19376
> fx3()
[1] -129.85674  -87.10575  -96.05648   52.80472  -66.53405
  • 指数乱数も
fx.rexp <- cxxfunction(signature(N="integer",M="numeric"),'int n=as<int>(N);double m=as<double>(M);return rexp(n, m);',plugin="Rcpp")
fx.rexp.2 <- cxxfunction(signature(N="integer",M="numeric"),'RNGScope z;int n=as<int>(N);double m=as<double>(M);return rexp(n, m);',plugin="Rcpp")
x<-fx.rexp(10000,1)
mean(x)
var(x)

> x<-fx.rexp(10000,1)
> mean(x)
[1] 0.9938479
> var(x)
[1] 0.9802089
  • ポアッソン乱数は
fx.rpois <- cxxfunction(signature(N="integer",M="numeric"),'int n=as<int>(N);double m=as<double>(M);return rpois(n, m);',plugin="Rcpp")
fx.rpois.2<-cxxfunction(signature(N="integer",M="numeric"),'RNGScope z;int n=as<int>(N);double m=as<double>(M);return rpois(n, m);',plugin="Rcpp")
x.pois<-fx.rpois(10000,1)
table(x.pois)

> table(x.pois)
x.pois
   0    1    2    3    4    5    6 
3617 3728 1838  642  138   36    1