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 <
[email protected]> <[email protected]> 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 <
> [email protected]> 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 <
> [email protected]> 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
>> <[email protected]> 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
>> > <[email protected]> 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
>> >> <[email protected]> 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 <
>> >> > [email protected]> 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]]
>> >> >
>> >> > ______________________________________________
>> >> > [email protected] mailing list
>> >> > https://stat.ethz.ch/mailman/listinfo/r-help
>> >> > PLEASE do read the posting guide
>> >> > http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>
>> >> > and provide commented, minimal, self-contained, reproducible code.
>> >
>> >
>>
>
>
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
