Problem in while loop

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 at 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 at gmail.com> <michael.weylandt at 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 at 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 at 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 at 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 at 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 at 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 at 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]]