[R] Problem in while loop
R. Michael Weylandt
michael.weylandt at gmail.com
Wed Dec 7 05:14:04 CET 2011
I'm not here to walk you through learning how to plot: but once again
does this work?
curve(loglikelihood(x, bout = 0.1))
If you want to plot a function of two independent variables, you'll
have to look into the various ways of doing so. Please do read section
12 of "An Introduction to R" which you can access by typing
help.start() at the prompt for the basics of graphics in R.
Michael
On Tue, Dec 6, 2011 at 9:47 PM, Gyanendra Pokharel
<gyanendra.pokharel at gmail.com> wrote:
>
> Thanks Michael,
> I am able to find the very nice plot of the function we discussed but still
> have the problem of ploting the function loglikelihood(aout,bout) versus
> aout posted in the initial massage.
> Best
> 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.
>>> >>> >
>>> >>> >
>>> >>
>>> >>
>>> >
>>
>>
>
More information about the R-help
mailing list