# [R] Problem in while loop

R. Michael Weylandt michael.weylandt at gmail.com
Wed Dec 7 01:02:35 CET 2011

```You want a plot of s = sum(log(b^2 + (x-a)^2)) as a function of  a or
am I missing something? You do it just like any other R plot: pick
some values of a, evaluate s for each of them (a little tricky if you
use my vectorized version due to the implicit use of the recycling
rule but very easy with your looped version)and plot(a,
loglikelihood(a)). Of course, since log is monotonic, it's pretty easy
to get an analytic solution for this.

If you want to be more creative about the graph this would also work:

x.fix <- c(-4.2, -2.85, -2.3, -1.02, 0.7, 0.98, 2.72, 3.5)
b <- 0.1

curve(sum(log(b^2 + (x.fix - x)^2))) # Note that curve has to take a
function of "x" so you have to rewrite it a little.

Michael

On Tue, Dec 6, 2011 at 10:40 AM, Gyanendra Pokharel
<gyanendra.pokharel at gmail.com> wrote:
> Yes Michael, it works well and I got the result what I want but it totally
> depends on how reliable result do I want. When I take very high rho (near
> about 1) and very low psi, it takes very long time may be it gives us more
> accurate result. But for lower rho and higher psi, it gives immediately, and
> very poor result. Now as I asked before, how I can get the plot of "a"
> versus "s" which determines the cooling curve with the global minimum of s
> at a specific value of a.
>
> On Tue, Dec 6, 2011 at 6:52 AM, R. Michael Weylandt
> <michael.weylandt at gmail.com> wrote:
>>
>> So I just ran your code verbatim with this one change and it finished
>> in less than 10 seconds. However, even without the change it doesn't
>> take more than 15 seconds: what exactly lead you to believe this was
>> an infinite loop?
>>
>> Michael
>>
>> On Tue, Dec 6, 2011 at 12:03 AM, R. Michael Weylandt
>> <michael.weylandt at gmail.com> <michael.weylandt at gmail.com> wrote:
>> > 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]]
>> >>> >> >
>> >>> >> > ______________________________________________
>> >>> >> > R-help at 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.
>> >>> >
>> >>> >
>> >>
>> >>
>> >
>
>

```