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.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
>> >> > 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.
