I ma sorry, I miss typed the function, it should be "loglikelihood" instead.
On Mon, Dec 5, 2011 at 10:37 PM, Gyanendra Pokharel <
gyanendra.pokha...@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.weyla...@gmail.com> <michael.weyla...@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.pokha...@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.weyla...@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.pokha...@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.weyla...@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.pokha...@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.pokha...@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@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<http://www.r-project.org/posting-guide.html>
>>> >> > and provide commented, minimal, self-contained, reproducible code.
>>> >
>>> >
>>>
>>
>>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
