2014-01-10 31 views
0

對於給定的觀察向量x和給定的參數向量(長度爲5),θ0,下面的代碼產生這些參數θ1的更好估計的向量。R中的迭代

part1=(1-theta0[5])*dnorm(x,theta0[1],theta0[3]) 
part2=theta0[5]*dnorm(x,theta0[2],theta0[4]) 
gam=part2/([part1+part2) 
denom1=sum(1-gam) 
denom2=sum(gam) 
mu1=sum((1-gam)*x)/denom1 
sig1=sqrt(sum((1-gam)*((x-mu1)^2))/denom1) 
mu2=sum(gam*x)/denom2 
sig2=sqrt(sum(gam*((x-mu2)^2))/denom2 
p=mean(gam) 
theta[1] <- c(mu1,mu2,sig1,sig2,p) 

我的問題是如何迭代過程而不必重新編寫代碼?我希望theta1取代原始估計的矢量,theta0,theta2來代替theta1等。收斂速度相當慢,所以我假設需要執行10000次迭代。有一個方便的捷徑嗎?

PS正如你可能會說,我是R的初學者,所以請儘量保持簡單。謝謝。

+0

您可以使用for循環進行固定次數的迭代。基本的語法是:for(i in 1:10000){do something}。 如果您可以測試收斂,則可以使用while循環來檢查每次迭代是否達到足夠的收斂。這最大限度地減少了計算時間(但要確保你的算法確實收斂!)。 –

+0

@AlexanderVosdeWael我將如何替換我的參數向量? – JohnK

+0

閱讀關於[申請](http://stackoverflow.com/questions/3505701/r-grouping-functions-sapply-vs-lapply-vs-apply-vs-tapply-vs-by-vs-aggrega/7141669# 7141669)功能。 – zx8754

回答

1

未經檢驗

theta0 <- c() 
x <- 

for (i in 1:10000) { 
part1 = (1 - theta0[5])*dnorm(x,theta0[1],theta0[3]) 
part2 = theta0[5] * dnorm(x,theta0[2],theta0[4]) 
gam = part2/(part1 + part2) 
denom1 = sum(1-gam) 
denom2 = sum(gam) 
mu1 = sum((1-gam)*x)/denom1 
sig1 = sqrt(sum((1-gam)*((x-mu1)^2))/denom1) 
mu2 = sum(gam*x)/denom2 
sig2 = sqrt(sum(gam * ((x - mu2)^2))/denom2) 
p = mean(gam) 

theta0 <- c(mu1,mu2,sig1,sig2,p) 
if (i %% 500 == 0) print(theta0) 
} 

真的沒有太大變化。將它包裝成一個循環,確保保存新的估計以替換舊的估計,以便循環執行某些操作。你有一些缺失的括號和拼寫錯誤。不要害怕使用空格,並確保關閉所有開放的([, {等。

+0

好的,謝謝。所以當我輸入theta0時,我會得到10000次迭代的估計,每次都替換參數向量? – JohnK

+1

是的。你將覆蓋0到10000次(或9999,我在數學上不好):)。如果你想在每次迭代中看到它,請添加'print(theta0)'(不推薦)。或者你可以說'if(i %% 500 == 0)print(theta0)'。那可能會更好 – rawr

+0

好了,它現在起作用了。非常感謝你! – JohnK