Off the bat I'd suggest you vectorize loglikelihood as a simple one liner:

sum(log(b^2 + (x-a)^2))

That alone will speed up your function many times over: I'll look at the big function in more detail tomorrow. 

Michael

On Dec 5, 2011, at 10:37 PM, Gyanendra Pokharel <gyanendra.pokharel@gmail.com> wrote:

> 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.weylandt@gmail.com> <michael.weylandt@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.pokharel@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.weylandt@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.pokharel@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.weylandt@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.pokharel@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.pokharel@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
>> >> > and provide commented, minimal, self-contained, reproducible code.
>> >
>> >
>> 
> 

	[[alternative HTML version deleted]]

