0

For a given vector of observations x and a given vector of parameters(of length 5) ,theta0, the code below produces a vector of better estimates of these parameters, theta1.

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)

My question is how can I iterate the procedure without having to write the code anew? I would like theta1 to replace my vector of original estimates, theta0, theta2 to replace theta1 etc etc. The rate of convergence is rather slow, so I suppose 10000 iterations will have to be performed. Is there a convenient shortcut?

PS As you can probably tell, I am a beginner in R so please keep it as simple as possible. Thank you.

Matthew Lundberg
  • 42,009
  • 6
  • 90
  • 112
JohnK
  • 1,019
  • 3
  • 14
  • 30
  • You can use a for-loop for a fixed number of iterations. The basic syntax is: for(i in 1:10000) { do something }. If you can test for convergence, you could use a while loop that checks if sufficient convergence is reached on each iteration. This minimizes computation time (but make sure your algorithm does converge at all!). – Alexander Vos de Wael Jan 10 '14 at 19:30
  • @AlexanderVosdeWael How would I replace my parameter vector that way though? – JohnK Jan 10 '14 at 19:32
  • Read on about [apply](http://stackoverflow.com/questions/3505701/r-grouping-functions-sapply-vs-lapply-vs-apply-vs-tapply-vs-by-vs-aggrega/7141669#7141669) functions. – zx8754 Jan 10 '14 at 22:41

1 Answers1

1

untested

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)
}

Really didn't change much. Wrap it in a loop, make sure you save the new estimates to replace the old ones so your loop does something. You had some missing parentheses and typos. Don't be afraid to use spaces, and make sure you close all opening (, [, {, etc.

rawr
  • 20,481
  • 4
  • 44
  • 78
  • Okay thank you. So when I input theta0, I will get the estimates of 10000 iterations, having replaced the vector of parameters each time? – JohnK Jan 10 '14 at 19:46
  • 1
    Yes. You are overwriting theta0 10000 times (or 9999, I'm bad at math :)). If you want to see it at each iteration, add `print(theta0)` (not recommended). Or you could say `if (i %% 500 == 0) print(theta0)`. That might be better – rawr Jan 10 '14 at 19:49