R 的 vectorization, again
Submitted by kzeng on Sun, 2011-03-27 21:26思考的方式，大概会受一个人正在说的语言影响，写 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++ 的影响吧。看来习惯也是可以改的：）
Free Tags:
Blog分类: