4

我不知道为什么这样的事情应该很慢：

steps=500
samples=100000
s_0=2.1
r=.02
sigma=.2
k=1.9

at<-matrix(nrow=(steps+1),ncol=samples)
at[1,]=s_0

for(j in 1:samples)
{
  for(i in 2:(steps+1))
  {
    at[i,j]=at[(i-1),j] + sigma*sqrt(.0008)*rnorm(1)
  }
}

我尝试使用 sapply 重写它，但从性能的角度来看它仍然很糟糕。

我在这里错过了什么吗？这将是 c++ 甚至是臃肿的 c# 中的秒数。

4

3 回答 3

4

R 可以向量化某些操作。在您的情况下，您可以通过进行以下更改来摆脱外循环。

for(i in 2:(steps + 1))
{
    at[i,] = at[(i - 1),] + sigma * sqrt(.0008) * rnorm(samples)
}

按system.time原版计算samples = 1000耗时6.83s，修改后耗时0.09s。

于 2013-02-13T21:49:29.257 回答
4

怎么样：

at <- s_0 + t(apply(matrix(rnorm(samples*(steps+1),sd=sigma*sqrt(8e-4)),
                   ncol=samples),
                    2,
                    cumsum))

（还没有仔细测试过，但我认为它应该是正确的，而且速度更快。）

于 2013-02-13T21:56:52.753 回答
1

要编写快速的 R 代码，您确实需要重新考虑如何编写函数。您想要对整个向量进行操作，而不仅仅是一次观察单个观察值。

如果您真的不喜欢编写 C 风格的循环，您也可以尝试 Rcpp。如果您非常习惯于 C++ 并且更喜欢以这种方式编写函数，那么可能会很方便。

library(Rcpp)
do_stuff <- cppFunction('NumericMatrix do_stuff(
  int steps,
  int samples,
  double s_0,
  double r,
  double sigma,
  double k ) {

  // Ensure RNG scope set
  RNGScope scope;

  // allocate the output matrix
  NumericMatrix at( steps+1, samples );

  // fill the first row
  for( int i=0; i < at.ncol(); i++ ) {
    at(0, i) = s_0;
  }

  // loop over the matrix and do stuff
  for( int j=0; j < samples; j++ ) {
    for( int i=1; i < steps+1; i++ ) {
      at(i, j) = at(i-1, j) + sigma * sqrt(0.0008) * R::rnorm(0, 1);
    }
  }

  return at;

}')

system.time( out <- do_stuff(500, 100000, 2.1, 0.02, 0.2, 1.9) )

给我

   user  system elapsed 
  3.205   0.092   3.297

因此，如果您已经具备一些 C++ 背景，请考虑学习如何使用 Rcpp 将数据映射到 R 或从 R 映射数据。

于 2013-02-14T00:51:10.163 回答