Thanks Michael Lets figure out the problem by using the following function. I found the same problem in this code too.
loglikehood <- function(a, b = 0.1, x = c(-4.2, -2.85, -2.3, -1.02, 0.7, 0.98, 2.72, 3.5)) { s <- 0 for(i in 1:length(x)){ s <- s + log(b^2 + (x[i] - a)^2) } s } loglikelihood(0.1) simann <- function(T0 = 1, N = 500, rho = 0.9, x0 = 0, eps = 0.1, f){ moving <- 1 count <- 0 Temp <- T0 x <- x0 while(moving > 0){ moving <- 0 for(i in 1:N){ y <- x + runif(1,-eps,eps) alpha <- min(1,exp((f(x) -f(y))/Temp)) if(runif(1)<alpha){ moving <- moving +1 x <- y } } Temp <- Temp*rho count <- count + 1 } fmin <- f(x) return(c(x,fmin,count)) } simann(f = loglikelihood) I hope we can analyze every problems from this function On Mon, Dec 5, 2011 at 10:27 PM, R. Michael Weylandt < michael.weyla...@gmail.com> <michael.weyla...@gmail.com> wrote: > It's not necessarily equivalent to your "loglikelihood" function but > since that function wasn't provided I couldn't test it. > > My broader point is this: you said the problem was that the loop ran > endlessly: I showed it does not run endlessly for at least one input so at > least part of the problem lies in loglikelihood, which, being unprovided, I > have trouble replicating. > > My original guess still stands: it's either 1) a case of you getting stuck > at probaccept = 1 or 2) so slow it just feels endless. Either way, my > prescription is the same: print() > > Michael > > > On Dec 5, 2011, at 9:30 PM, Gyanendra Pokharel < > gyanendra.pokha...@gmail.com> wrote: > > Yes, your function out<- epiann(f = function(a,b) > log(dnorm(a)*dnorm(b))), N = 10) works well. > > But why you are changing the loglikelihood function to f = function(a,b) > log(dnorm(a)*dnorm(b))? how it is equivalent to loglikelihood? is there any > mathematical relation? I also want to see the plot of aout and bout versus > loglikelihood, and see the cooling rate in different temperature. how is > this possible? > > On Mon, Dec 5, 2011 at 6:07 PM, R. Michael Weylandt < > michael.weyla...@gmail.com> wrote: > >> If you run >> >> out<- epiann(f = function(a,b) log(dnorm(a)*dnorm(b))), N = 10) >> >> It takes less than 0.5 seconds so there's no problem I can see: >> perhaps you want to look elsewhere to get better speed (like Rcpp or >> general vectorization), or maybe your loglikihood is not what's >> desired, but there's no problem with the loop. >> >> Michael >> >> On Mon, Dec 5, 2011 at 5:29 PM, Gyanendra Pokharel >> <gyanendra.pokha...@gmail.com> wrote: >> > Yes, I checked the acceptprob, it is very high but in my view, the while >> > loop is not stopping, so there is some thing wrong in the use of while >> loop. >> > When I removed the while loop, it returned some thing but not the result >> > what I want. When i run the while loop separately, it never stops. >> > >> > >> > >> > On Mon, Dec 5, 2011 at 5:18 PM, R. Michael Weylandt >> > <michael.weyla...@gmail.com> wrote: >> >> >> >> Your code is not reproducible nor minimal, but why don't you put a >> >> command print(acceptprob) in and see if you are getting reasonable >> >> values. If these values are extremely low it shouldn't surprise you >> >> that your loop takes a long time to run. >> >> >> >> More generally, read up on the use of print() and browser() as >> debugging >> >> tools. >> >> >> >> Michael >> >> >> >> On Mon, Dec 5, 2011 at 3:47 PM, Gyanendra Pokharel >> >> <gyanendra.pokha...@gmail.com> wrote: >> >> > I forgot to upload the R-code in last email, so heare is one >> >> > >> >> > epiann <- function(T0 = 1, N=1000, ainit=1, binit=1,rho = 0.99, >> amean = >> >> > 3, >> >> > bmean=1.6, avar =.1, bvar=.1, f){ >> >> > >> >> > moving <- 1 >> >> > count <- 0 >> >> > Temp <- T0 >> >> > aout <- ainit >> >> > bout <- binit >> >> > while(moving > 0){ >> >> > moving <- 0 >> >> > for (i in 1:N) { >> >> > aprop <- rnorm (1,amean, avar) >> >> > bprop <- rnorm (1,bmean, bvar) >> >> > if (aprop > 0 & bprop > 0){ >> >> > acceptprob <- min(1,exp((f(aout, bout) - >> >> > f(aprop,bprop))/Temp)) >> >> > u <- runif(1) >> >> > if(u<acceptprob){ >> >> > moving <- moving +1 >> >> > aout <- aprop >> >> > bout <- bprop >> >> > } >> >> > else{aprob <- aout >> >> > bprob <- bout} >> >> > } >> >> > } >> >> > Temp <- Temp*rho >> >> > count <- count +1 >> >> > >> >> > } >> >> > fmin <- f(aout,bout) >> >> > return(c(aout, bout,fmin, count) ) >> >> > >> >> > } >> >> > out<- epiann(f = loglikelihood) >> >> > >> >> > On Mon, Dec 5, 2011 at 3:46 PM, Gyanendra Pokharel < >> >> > gyanendra.pokha...@gmail.com> wrote: >> >> > >> >> >> Hi all, >> >> >> I have the following code, >> >> >> When I run the code, it never terminate this is because of the while >> >> >> loop >> >> >> i am using. In general, if you need a loop for which you don't know >> in >> >> >> advance how many iterations there will be, you can use the `while' >> >> >> statement so here too i don't know the number how many iterations >> are >> >> >> there. So Can some one suggest me whats going on? >> >> >> I am using the Metropolis simulated annealing algorithm >> >> >> Best >> >> >> >> >> > >> >> > [[alternative HTML version deleted]] >> >> > >> >> > ______________________________________________ >> >> > R-help@r-project.org mailing list >> >> > https://stat.ethz.ch/mailman/listinfo/r-help >> >> > PLEASE do read the posting guide >> >> > http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html> >> >> > and provide commented, minimal, self-contained, reproducible code. >> > >> > >> > > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.