- こちらからの続き
- 乱数を発生するシミュレーションをしてみる
- 基本は一様整数乱数を発生させるrand()と計算機の時刻からシード整数を取るsrand()
- rand()は規定の大きな整数を最大値とする一様整数乱数を作るので、0-1の実数一様乱数を作るには
r = (double)rand()/((double)RAND_MAX+1);
srand(time(NULL)); // 現在時刻を乱数の種の設定
- これをRcpp,inlineで使うには、こちらにあるように、
#include <ctime> // for time()
#include <cstdlib> // for srand(), rand()
- と、関数をインクルードする必要があり、それをinlineで使うには、cxxfunction()の引数includeに渡す
- このようにすることで、こちらでは乱数をRで発生させて、C++に渡していたものを、C++の中で乱数発生させることになった
- さて、Rの強みは、いろいろな確率分布からの乱数発生が簡単なことだが、C++はその点が、「がっかり」な感じがする
- そもそも標準装備ではない
- アドオンとしてとってくるにしても、Rでの簡単作業『CRANで検索、パッケージ名を見つけたら、install.packages()』に慣れている身にはつらい
- そもそも、C、C++がどこにどう入ったのかもわからずに使っていることが問題なのだが
- Rのrnorm()関数は、以下のように.Internal呼び出しなので、どこかに、C/Fortranなどで書いたソースが埋まっているはず。それを呼び出せればよいように思うのだが…
- Rのこの関数をStandard Template Library式に呼び出して使っても、遜色ない結果が出るのなら、それでもよいだろう…
> rnorm
function (n, mean = 0, sd = 1)
.Internal(rnorm(n, mean, sd))
<bytecode: 0x02ab3c90>
<environment: namespace:stats>
my.diff2 <- cxxfunction(sig=signature(A="integer",B="integer"), include=c("#include <ctime>","#include <cstdlib>"),body=my.body2, plugin="Rcpp")
library(Rcpp)
library(inline)
my.body2 <- '
int a=as<int>(A);
int b=as<int>(B);
double r;
double r2;
double p;
Rcpp::NumericMatrix y(a,b);
srand(time(NULL)); // 現在時刻を乱数の種の設定
for(int j=0; j<b; j++){
y(0,j) = 0;
}
for (int i = 1; i < a; i++){
for(int j = 0; j<b; j++){
//r = rand();
r = (double)rand()/((double)RAND_MAX+1);
r2 = (double)rand()/((double)RAND_MAX+1);
p=1;
if(r2 < 0.5){
p =-1;
}
y(i,j) = y(i-1,j)+r*p;
}
}
return(y);
'
my.diff2 <- cxxfunction(sig=signature(A="integer",B="integer"), include=c("#include <ctime>","#include <cstdlib>"),body=my.body2, plugin="Rcpp")
A<-100000
A.less<-100
B<-4
my.out2<-my.diff2(A,B)
library(rgl)
plot3d(my.out2[1:A.less,1:3],cex=0.1)
plot(as.data.frame(my.out2),cex=0.01)
