僅供參考:自從我的第一版以來,我已經對其進行了重大編輯。這種模擬從14小時減少到14分鐘。R中的代碼更快
我是編程新手,但我做了一個模擬試圖跟隨生物體中的無性複製,並量化父母和子女生物之間染色體數量的差異。模擬運行非常緩慢。大約需要6個小時才能完成。我想知道什麼是使模擬運行更快的最佳方法。
這些數字生物具有x個染色體。與大多數生物體不同,染色體彼此獨立,因此它們有可能被轉移到任何一個女性生物體中。
在這種情況下,染色體在子細胞中的分佈遵循二項分佈,概率爲0.5。
函數sim_repo
取一個具有已知染色體數目的數字生物矩陣,並將它們通過12代複製。它複製這些染色體,然後使用rbinom
函數隨機生成一個數字。然後將該號碼分配給女兒單元。由於在無性生殖期間沒有染色體丟失,所以其他子細胞接收剩餘的染色體。然後這個重複的代數爲G。然後從矩陣的每一行中採樣一個值。
我我的研究很感興趣,在最初的祖先生物,在這個模擬的最後時間點的染色體變異中的變化。以下功能代表將細胞轉移到新的生活環境中。它採用函數sim_re
p的輸出並使用它來生成更多代。然後它找出第一個和最後一個矩陣列中的行之間的差異,並找出它們之間的差異。
# 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_exp
h次數。
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
乍一看,我敢打賭,最大的原因你的代碼是緩慢的是,你不預先分配你的對象:特別是'sim_1000()'裏的'sim_exp()'和'c()'裏面的'cbind()'必須非常昂貴。 – flodel 2012-03-07 02:20:32
@ flodel,謝謝你的提示。你有一個例子如何預先分配在我的代碼?例如,在'sim_exp()'中,我會在最終輸出中創建一個與我預期的相同數量的列和行的矩陣,但是用NULL來填充值? – Kevin 2012-03-07 02:50:09
The R Inferno的一章專門爲此而設:http://www.burns-stat.com/pages/Tutor/R_inferno.pdf – 2012-03-07 03:01:03