2012-03-07 95 views
6

僅供參考:自從我的第一版以來,我已經對其進行了重大編輯。這種模擬從14小時減少到14分鐘。R中的代碼更快

我是編程新手,但我做了一個模擬試圖跟隨生物體中的無性複製,並量化父母和子女生物之間染色體數量的差異。模擬運行非常緩慢。大約需要6個小時才能完成。我想知道什麼是使模擬運行更快的最佳方法。

這些數字生物具有x個染色體。與大多數生物體不同,染色體彼此獨立,因此它們有可能被轉移到任何一個女性生物體中。

在這種情況下,染色體在子細胞中的分佈遵循二項分佈,概率爲0.5。

函數sim_repo取一個具有已知染色體數目的數字生物矩陣,並將它們通過12代複製。它複製這些染色體,然後使用rbinom函數隨機生成一個數字。然後將該號碼分配給女兒單元。由於在無性生殖期間沒有染色體丟失,所以其他子細胞接收剩餘的染色體。然後這個重複的代數爲G。然後從矩陣的每一行中採樣一個值。

​​

我我的研究很感興趣,在最初的祖先生物,在這個模擬的最後時間點的染色體變異中的變化。以下功能代表將細胞轉移到新的生活環境中。它採用函數sim_rep的輸出並使用它來生成更多代。然後它找出第一個和最後一個矩陣列中的行之間的差異,並找出它們之間的差異。

# The following function is mostly the same as I talked about in the description. 
    # The only difference is I changed some aspects to take into account I am using 
    # matrices and not lists. 
    # The function outputs the difference between the intial variance component between 
    # 'cell lines' with the final variance after t number of transfers 

sim_exp = function(x1, G=12, k=1, t=25, h=1000) { 

    xn <- matrix(NA, nrow(x1), t) 
    x <- x1 
    xn[,1] <- x1 
    for (l in 2:t) { 
     x <- sim_repo(x, G, k, t, h) 
     xn[, l] <- x 
    } 

    colvar <- matrix(apply(xn,2,var),ncol=ncol(xn)) 
    ivar <- colvar[,1] 
    fvar <- colvar[,ncol(xn)] 
    deltavar <- fvar - ivar 
    deltavar 
} 

我需要複製這個模擬^h的次數。因此,我做了以下功能,將調用功能sim_exph次數。

sim_1000 = function(x1, G=12, k=1, t=25, h=1000) { 
    xn <- vector(length=h) 
    for (l in 2:h) { 
     x <- sim_exp(x1, G, k, t, h) 
     xn[l] <- x 
    } 
     xn 
} 

當我打電話與像6倍的值需要花費大約52秒來完成的值sim_exp功能。

x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) 
system.time(sim_1000(x1,h=1)) 
    user system elapsed 
    1.280 0.105 1.369 

如果我能得到更快然後我可以完成更多的這些模擬和仿真應用選擇模型。

我的輸入將類似於以下x1,在其自己的行每祖生物基質:

x1 <- matrix(data=c(100,100,100,100,100,100),ncol=1) # a matrix of 6 organisms 

當我運行:

a <- sim_repo(x1, G=12, k=1) 

我的預期產出將是:

a 
    [,1] 
[1,] 137 
[2,] 82 
[3,] 89 
[4,] 135 
[5,] 89 
[6,] 109 

system.time(sim_repo(x1)) 
    user system elapsed 
    1.969 0.059 2.010 

當我調用sim_exp函數時,

b < - sim_exp(X1,G = 12,k = 1時,T = 25)

它調用sim_repo函數G倍和輸出:

b 
[1] 18805.47 

當我調用sim_1000功能,我通常會將h設置爲1000,但這裏我將它設置爲2.所以這裏sim_1000將調用sim_exp並複製它兩次。

c <- sim_1000(x1, G=12, k=1, t=25, h=2) 
c 
[1] 18805.47 18805.47 
+0

乍一看,我敢打賭,最大的原因你的代碼是緩慢的是,你不預先分配你的對象:特別是'sim_1000()'裏的'sim_exp()'和'c()'裏面的'cbind()'必須非常昂貴。 – flodel 2012-03-07 02:20:32

+0

@ flodel,謝謝你的提示。你有一個例子如何預先分配在我的代碼?例如,在'sim_exp()'中,我會在最終輸出中創建一個與我預期的相同數量的列和行的矩陣,但是用NULL來填充值? – Kevin 2012-03-07 02:50:09

+0

The R Inferno的一章專門爲此而設:http://www.burns-stat.com/pages/Tutor/R_inferno.pdf – 2012-03-07 03:01:03

回答

8

正如在評論中提到別人,如果我們在函數sim_repo只能看,並更換行:

dup <- x1 * 2 

dup <- apply(x1, c(1,2),"*",2) 

z <- apply(dup,c(1,2),rbinom,n=1,prob=0.5) 

z <- matrix(rbinom(n=rep(1,length(dup)),size = as.vector(dup),prob = 0.5),nrow = nrow(dup)) 

和內部的循環與

x1 <- matrix(apply(x1,1,sample,size = 1), ncol=1) 

,我收到了,好了,大的速度增加:

system.time(sim_exp(x1)) 
    user system elapsed 
    0.655 0.017 0.686 
> system.time(sim_expOld(x1)) 
    user system elapsed 
21.445 0.128 21.530 

我驗證了它在做同樣的事情:

set.seed(123) 
out1 <- sim_exp(x1) 

set.seed(123) 
out2 <- sim_expOld(x1) 

all.equal(out1,out2) 
> TRUE 

而這甚至沒有削減進入預先分配,根據您構建代碼的方式,如果不完全重新設計,可能實際上很難實現。

而這也甚至沒有開始看你是否真的需要所有三種功能...

+0

我需要使用您的電腦。我仍然得到: 'system.time(sim_exp(x1,G = 12,k = 1,t = 25,h = 1))' '用戶系統過去了' '23.598 0.767 24.390' – Kevin 2012-03-07 05:05:03

+0

@Kev My電腦不快。這是一臺一年前的macbook air。隨着兩個處理器選項中較慢的。更有可能你沒有完全正確地修改代碼。 – joran 2012-03-07 05:15:20

+0

你是對的,忘了適用於那個for循環。 – Kevin 2012-03-07 06:09:43