Problem in while loop
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.