2014-04-23 31 views
1

我想重複該過程,除非在某個時間存儲結果時滿足條件。循環在r中,直到值收斂並存儲所有輸出

下面是一個簡單的例子,我知道循環次數在循環執行:

# just example data 
smpls <- rnorm(100,50,50) 

ncycles <- 1000 
outm <- matrix(nrow=ncycles, ncol = 1) 

#repeate過程n個循環

for(i in 1:ncycles){  
outm[i] <- mean(sample(smpls, 50)) 
} 
# get average of outm 
outm <- mean(sample(smpls, 50)) 

但是我的情況是在這個意義不同我不認識ncyles。我想繼續採樣,除非樣本會變得非常低方差或收斂(我猜是「while」循環。例如,除非VSD在以下情況下小於1。

vsd <- NULL 
outm <- mean(sample(smpls, 50)) 
while (vsd > 1){ 
    outm[i] <- mean(sample(smpls, 50)) 
    vsd <- sd(outm) 
    } 

我不知道該值我在這裏要設置幫助表示讚賞

編輯:。

smpls <- rnorm(100,50,50) 
iter <- 0 
# maximum iteration 
itermax <- 1000 
outm <- rep(NA, itermax) 
vsd <- 2 
while((vsd > 1) && (iter < itermax)) { 
    outm[iter] <- mean(sample(smpls, 50)) 
    vsd <- sd(outm) 
    iter <- iter+1 
    } 
Error in while ((vsd > 1) && (iter < itermax)) { : 
    missing value where TRUE/FALSE needed 

當它達到收斂是爲了節省時間停止主要思想雖然ABO作爲一個簡單的平均函數的例子很快,我原來的函數需要大量的時間來做迭代,我想在收斂時停下來。在你的代碼

+2

你只要循環爲'標準偏差outm'大於1。然而,你的第一次迭代'outm'將有長度後1和標準差0,導致您的循環終止。 – josliber

+2

開始outm'作爲一個長矢量,說10K。如果迭代索引超出範圍,請在循環中檢查另一個10K的循環。收斂後修剪NA。不用說,你需要初始化並增加'i'。這將最大限度地減少'outm'從零開始重寫的次數。 – ilir

+0

看看是否有幫助:https://stat.ethz.ch/pipermail/r-help/2010-June/241386.html – SHRram

回答

1

這裏是一個解決方案:

數據

set.seed(123) # so that you can replicate what I did 
smpls <- rnorm(100,50,50) 

我想你需要一些初始化週期(最小迭代),使您會收到假收斂,因爲你有少量的樣品。所以運行幾個樣本 - 比如說小米。你還需要一個最大的迭代,以便你的循環不會變得瘋狂 - 比如說maxiter。

meanconverge <- function (data, miniter, maxiter, tolerance){ 
     outm <- rep(NA, maxiter) 
    for(i in 1:miniter){  
    outm[i] <- mean(sample(smpls, 50)) 
    } 
    # sd of initial cycles 
    vsd <- sd(outm, na.rm = TRUE) 
    if(vsd > tolerance) { 
        iter <- miniter+1 
        sdout <- rep(NA, maxiter) 
        while((vsd > tolerance) && (iter < maxiter)) { 
        iter <- iter + 1 
        outm[iter] <- mean(sample(smpls, 50)) 
        vsd <- sd(outm, na.rm = TRUE) 
        sdout[iter] <- vsd    
     } 
     out <- list(outm, sdout) 
     return(out) 
     } else { 
     return(outm) 
     } 
     } 

out <- meanconverge (data = smpls, miniter = 50, maxiter = 100000, tolerance = 3) 
plot(unlist(out[2]), pch = ".", col = "red") 

enter image description here

plot(unlist(out[1]), pch = ".", col = "red") 

enter image description here

2

兩個問題:

1)你需要SD(...,na.rm = TRUE)

2)你必須確保至少有兩個數字在outd爲sd(outm,na.rm = TRUE)!=不適用

順便說一句,給定sd你指定rnorm,我不認爲你需要超過幾十次迭代

sim <- function() { 
    smpls <- rnorm(100,50,5) 
    itermax <- 1000 
    outm <- rep(NA, itermax) 
    outm[1] <- mean(sample(smpls, 50)) 
    iter <- 1 
    vsd <- 2 
    while((vsd > 1) && (iter < itermax)) { 
     iter <- iter+1 
     outm[iter] <- mean(sample(smpls, 50)) 
     vsd <- sd(outm, na.rm = TRUE) 
     } 

    iter 
    } 

set.seed(666) 
iters <- replicate(100000, sim()) 
range(iters) # c(2, 11) 

乾杯。

+0

感謝您指出代碼中的兩個重要限制 - 但是我只是獲得了所有100000次迭代的值 - 雖然這個工作示例很小且很快 - 我需要做的迭代需要時間並且我想停止進一步的計算當價值達到收斂時 – jon

+0

@jon - 上面的第2-12行是否適合你?如果他們這樣做,你可以採取解決方案,並保持獎金(我沒有一個strackoverflow帳戶...)。如果它不起作用,讓我知道你看到了什麼,我會再去一次。 – Dale

+0

+1用於糾正問題(2-12行),但最初我不明白該函數的意義。 – jon

0

檢查收斂是一件棘手的事情。開始行動的一個好方法是觀察計算過程中值的變化。融合是關於任意接近邊界的;在編程上,你必須做出選擇「任意」的意思。您還需要決定如何衡量收斂。

爲了說明這一點,假設我想知道我的估計是否符合我的條件,是否真的很接近。我可能是這樣的:

# inside my function or method that performs this convergence feat 
while (while_condition && i < itermax)) { 

    outcome[i] <- some_complicated_foo(bar) 

    if (abs(outcomes[i-1] - outcomes[i]) <= tolerance) { 
     while_condition <- FALSE # i.e. STOP LOOPING 
     return outcomes 
    } 

    else {continue} 

} 

哪裏tolerance是你的任意親近的定義。現在,這對你的指甲看起來像是錘子吧?那麼,如果你收斂到錯誤的答案會發生什麼?你怎麼知道?這件事是否會收斂?這些問題的訣竅是讓切肉刀猜測你的功能或你正在分析的數據生成過程。但是,只要合理,具有最大迭代邊界將肯定會減少一些計算時間。知道自己是否正確的真正方法是使用測試(如統計測試或單元測試)來查看是否存在任何「垃圾進入垃圾」或獲得與您期望的不同之處一個衆所周知的答案。

查看優化算法並查看它們是如何實現的。請參閱?optim或其他優化包,瞭解專業人員如何執行此操作。