2013-02-21 89 views
0

我沒有太多的經驗,在R.我想寫一個Gibbs採樣在那裏我有一個循環是這樣的:如何向量化這個循環中的R

for (iNum in 1:totNum) { 
    rateNum <- Y3[iNum] 
    if(Y3[iNum] > 0) { 
     yStar3[iNum] <- rtnorm(1, mean = Mean3[iNum], sd = sqrt(Var3), 
       lower = gz[rateNum], upper = gz[rateNum + 1]) 
    } else if(Y3[iNum] == 0) { 
    yStar3[iNum] <- rtnorm(1, mean = Mean3[iNum], sd = sqrt(Var3), 
        lower = -Inf, upper = Inf); 
    } 
} 

這是服用過多的時間。我試圖使用lapply,但這還不夠快。有沒有一種方法來實現這個循環的矢量化?

謝謝你,最好的問候!

+1

我會說是的,因爲'rtnorm'被矢量,但你似乎對每個值使用不同的截斷值(上,下),而且我認爲這個參數沒有被矢量化。 – joran 2013-02-21 19:54:26

回答

0

這是一個沒有你的變量值的問題,但你想要做的事情是非常簡單的。你想堅持所有的矢量化語句,並儘量不要佔用過多的內存。這是基本策略:

第1步:找出如何計算所有數字。

# The number of values you need from 'rtnorm' 
sum(Y3 > 0) 
sum(Y3 == 0) 

# The means you need from the 'Mean3' array 
Mean3[Y3 > 0] 
Mean3[Y3 == 0] 

# Lower and upper limits for Y3 > 0 
gz[Y3[Y3 > 0]] 
gz[Y3[Y3 > 0] + 1] 

步驟2:上yStar3的向量濾波器使用這些值。我不能絕對肯定我所有的語法是完美的沒有一些樣本數據和變量值,但它應該是這個樣子:

yStar3[Y3 > 0] <- rtnorm(
    sum(Y3 > 0), 
    mean = Mean3[Y3 > 0], 
    sd = sqrt(Var3), 
    lower = gz[Y3[Y3 > 0]], 
    upper = gz[Y3[Y3 > 0] + 1]) 

yStar3[Y3 == 0] <- rtnorm(
    sum(Y3 == 0), 
    mean = Mean3[Y3 == 0], 
    sd = sqrt(Var3), 
    lower = -Inf, 
    upper = Inf) 
+0

感謝大家的好評!我已經能夠得到它的工作,並顯着提高速度。這一步的一次迭代現在需要0.05秒而不是1秒。5秒,實際上爲整個算法節省了數小時的計算時間! – user2096864 2013-02-22 10:18:52

+0

很高興我們可以提供幫助。正如你可以說的那樣,優化問題對我們很多人來說實際上是非常有趣的。下一次,如果您可以提供一些示例數據並要求優化而不是向量化,那麼您可能會對您的問題產生更多興趣和活動。 – Dinre 2013-02-22 18:19:23

+0

此外,如果您已達到解決方案,請繼續選擇您最滿意的答案,或者通過單擊其中一個答案旁邊的複選標記回答您的原始問題最有幫助。這將有助於未來找到此頁面的其他人並將結束該問題。 – Dinre 2013-02-25 12:36:36

2

因此,它不會出現,你必須迭代間的依賴性,這使得它非常簡單的向量化

lhs = rtnorm(length(Y3), mean = Mean3, sd = sqrt(Var3), lower = gz[Y3], 
       upper = gz[Y3 + 1]) 
    rhs = rtnorm(length(Y3), mean=Mean3, sd = sqrt(Var3), lower=-Inf, upper=Inf) 

    ifelse(Y3 > 0, lhs, rhs) 

這裏的問題是,rtnorm已經超過它的輸入參數進行矢量化,意味着,下,和上。情況可能並非如此,在這種情況下,你將不得不做更多的工作。

+0

打敗我吧。我相信'rtnorm'是msm包中的一個截斷正態分佈,所以它應該已經被矢量化了。 – user295691 2013-02-21 20:04:10

2

最簡單的方法是生成兩個條件的一半,並選擇你想要的。該mean參數將採取載體手段,使你得到的東西是這樣的:

yStar3 <- ifelse(
    Y3 > 0, 
    rtnorm(totNum, mean=Mean3, sd=sqrt(Var3), lower=gz[ratenum], upper=gz[rateNum+1]), 
    rtnorm(totNum, mean=Mean3, sd=sqrt(Var3), lower=-Inf, upper=Inf)) 

你必須與其中Y3小於零的情況下附加條件細化ifelse,潛在的,但是這是總體思路。

更新:@hadley表明移動ifelse的rtnorm內:

yStar3 <- rtnorm(totNum, mean=Mean3, sd=sqrt(Var3), 
    lower=ifelse(Y3>0,gz[rateNum], -Inf), 
    upper=ifelse(Y3>0,gz[rateNum+1], Inf)) 

現在有基本上是零不必要的計算回事。

更新:當然,1是錯誤的,正如評論者指出的;它應該是totNum。

+1

'ifelse'必須完全計算兩個語句;你最好在'rtnorm'裏面使用它' – hadley 2013-02-21 22:39:47

+0

你只是想到了我。 – user295691 2013-02-21 22:43:40

+1

您是否需要爲'rtnorm'功能指定除'1'以外的數字?第一個參數告訴它返回多少個值。 – Dinre 2013-02-22 01:19:23