R 的 vectorization, again

  • Posted on: 27 March 2011
  • By: kzeng

思考的方式,大概会受一个人正在说的语言影响,写 codes 似乎也是一样。如果是写 C++,思考的习惯大概会和 Matlab 不一样——不是说 syntax 或者是 concept,而是解决问题的抽象方式。

一个简单的题目,模拟一个 Markov Pure Jump Process (MPJP), transition matrix 是一个 8x8 的矩阵,模拟一百万次,先在 R 中写,习惯的先写一个1:8的loop 模拟一条路径,然后重复一百万次:

   1:  exp_sim=function(n){          
   2:    sim=array(0,c(n,64)) 
   3:    for (k in 1:n){ 
   4:      for(i in 1:8){ 
   5:        j=i 
   6:        t=rexp(1,rate=lam[j])           
   7:        while (t<5){ 
   8:          j=sample(1:8,1,prob=P[j,])    
   9:          t=t+rexp(1,rate=lam[j])       
  10:          }                             
  11:        sim[k,((i-1)*8+j)]=1 
  12:        } 
  13:      } 
  14:    result=array(0,c(8,8)) 
  15:    for (i in 1:8){ 
  16:      for (j in 1:8){ 
  17:      result[i,j]=mean(sim[,((i-1)*8+j)]) 
  18:      } 
  19:      } 
  20:      result 
  21:  }

结果消耗的时间可想而知,耗费了大约 340 秒才完成,实在太慢了,以为是 R 的问题,又在 matlab 里重新写过,运行大约耗费了 23 秒,正感叹 Matlab 速度快,优化作的好。。。可马上也意识到在 matlab 里写 codes 的思路竟然完全不同——根本没有去作那个 1:10^6 的 loop,而是按照习惯,一开始就生成了一个一百万行的 vector,避开了这个耗时的 loop。。。所以又参照 matlab 的 codes,重新写 R 的 codes:

   1:  exp_simv2=function(n){ 
   2:    mysample=function(vec){ 
   3:      m=length(vec) 
   4:      temp=numeric(m) 
   5:      for (i in 1:8){ 
   6:     temp=temp+(vec==i)*sample(1:8,m,prob=P[i,], replace=TRUE) 
   7:      } 
   8:      temp 
   9:      } 
  10:      
  11:      result=array(0,c(8,8)) 
  12:      for (i in 1:8) { 
  13:          t=array(0,n) 
  14:          state=i*array(1,n) 
  15:          t=t+rexp(n,rate=lam[state]) 
  16:          while (min(t)<5){ 
  17:              newstate=mysample(state) 
  18:              state=newstate*(t<=5)+state*(t>5) 
  19:              t=t+rexp(n,rate=lam[state]) 
  20:              } 
  21:              for (j in 1:8){ 
  22:                  result[i,j]=sum(state==j) 
  23:              } 
  24:          } 
  25:          result=result/n 
  26:          result 
  27:  }

像 matlab 一样思考,避开了那个一百万次的 loop,虽然写的时候有些别扭,但是运行以后发现,速度果然快了不少,效率大约提高了10倍,最后只耗时30秒左右,基本上和 matlab 的速度相仿。

习惯的力量真可怕 –_-#

P.S. 不过也发现了一个新的让人欣慰的习惯,不管是写 R 还是 matlab,已经不自觉习惯于 encapsulation,这肯定是受 C++ 的影响吧。看来习惯也是可以改的:)

Blog分类: