The loglikelihood() looks ok and gives some value. But I am using this
function for the simulated annealing, generating the random samples from
uniform proposal density., The codes are as follows


epiannea <- function(T0 = 1, N = 500,beta = 0.1,x0 = 0.1, rho = 0.90, eps =
0.1, loglikelihood, data){

moving <- 1

count <- 0

Temp <- T0

alpha <- x0

while(moving > 0){

mmoving <- 0

for(i in 1:N){

r <- alpha + runif(1,-eps,eps)

if(r > 0){

a <- min(1,exp((loglikelihood(alpha,beta)-loglikelihood(r,beta))/Temp))

}

if(runif(1)< a){

moving <- moving +1

alpha <- r

}

}

Temp <- Temp*rho

count <- count + 1

}

#plot(a,loglikelihood(alpha,beta), type = "l")

fmin <- loglikelihood(alpha, beta)

return(c(fmin, alpha, count))

}

epiannea(loglikelihood=loglikelihood)

Some times it shows the warnings that logp[i] produces NaNs, that means
p[i] is negative, but it should not be that because p[i] is the probabiliy
and can't be negative. Some times the code runs but never stop.
On Wed, Nov 30, 2011 at 7:34 AM, R. Michael Weylandt <
michael.weylandt@gmail.com> <michael.weylandt@gmail.com> wrote:

> I'd suggest you do some leg-work and figure out why you are getting values
> >1. If your algorithm is motivated by some approximation then a min() or
> pmin() *might* be the right fix, but if there are no approximations you may
> need to start debugging properly to see why you are getting an out of
> bounds value.
>
> Since there's no random number generation involved, I'd hesitate to just
> throw out the result without knowing its source. Also keep in mind the
> limitations of floating point arithmetic if you expect alpha*d^beta to be
> small.
>
> Michael
>
> On Nov 29, 2011, at 6:58 PM, Sarah Goslee <sarah.goslee@gmail.com> wrote:
>
> > On Tue, Nov 29, 2011 at 6:55 PM, Gyanendra Pokharel
> > <gyanendra.pokharel@gmail.com> wrote:
> >> yes, log of negative number is undefined and R also do the same and
> produces
> >> NaNs. Here I want to reject the value of exp(-alpha*d^(-beta)) when
> greater
> >> than 1, and want to run the loop otherwise.
> >> Thanks
> >
> > Then just add another if() statement checking for that condition.
> >
> >> On Tue, Nov 29, 2011 at 6:48 PM, Sarah Goslee <sarah.goslee@gmail.com>
> >> wrote:
> >>>
> >>>> Here p[i] <- 1 - exp(-alpha*d^(-beta))> so,  log(p[i]) produces NaNs
> >>>> when exp(-alpha*d^(-beta)) is greater than 1.> How can I remove
> it.After
> >>>> generating the out put we can omit it, but the> problem is different.
> >>>
> >>> Wait... you're complaining that you can't take the natural log of a
> >>> negative
> >>> number in R?
> >>>
> >>> You can't do that anywhere. What do you expect to happen? The log of a
> >>> negative number IS NaN.
> >>>
> >>> Sarah
> >>> On Tue, Nov 29, 2011 at 6:28 PM, Gyanendra Pokharel
> >>> <gyanendra.pokharel@gmail.com> wrote:
> >>>> I have following code:
> >>>> loglikelihood <- function(alpha,beta= 0.1){
> >>>>    loglh<-0
> >>>>    d<-0
> >>>>    p<-0
> >>>>    k<-NULL
> >>>>    data<-read.table("epidemic.txt",header = TRUE)
> >>>>    attach(data, warn.conflicts = F)
> >>>>    k <-which(inftime==1)
> >>>>    d <- (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
> >>>>    p<-1 - exp(-alpha*d)
> >>>>    for(i in 1:100){
> >>>>        if(i!=k){
> >>>>            if(inftime[i]==0){
> >>>>                loglh<-loglh +log(1-p[i])
> >>>>            }
> >>>>            if(inftime[i]==2){
> >>>>                loglh<-loglh + log(p[i])
> >>>>            }
> >>>>        }
> >>>>    }
> >>>>    return(loglh)
> >>>> }
> >>>> Here p[i] <- 1 - exp(-alpha*d^(-beta))
> >>>> so,  log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater
> than
> >>>> 1.
> >>>> How can I remove it.After generating the out put we can omit it, but
> the
> >>>> problem is different.
> >>>>
> >>>> On Tue, Nov 29, 2011 at 5:22 PM, Gyanendra Pokharel <
> >>>> gyanendra.pokharel@gmail.com> wrote:
> >>>>
> >>>>> No, that,s not a problem Michael,
> >>>>> I have following code:
> >>>>> loglikelihood <- function(alpha,beta= 0.1){
> >>>>>     loglh<-0
> >>>>>     d<-0
> >>>>>     p<-0
> >>>>>     k<-NULL
> >>>>>     data<-read.table("epidemic.txt",header = TRUE)
> >>>>>     attach(data, warn.conflicts = F)
> >>>>>     k <-which(inftime==1)
> >>>>>     d <- (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
> >>>>>     p<-1 - exp(-alpha*d)
> >>>>>     for(i in 1:100){
> >>>>>         if(i!=k){
> >>>>>             if(inftime[i]==0){
> >>>>>                 loglh<-loglh +log(1-p[i])
> >>>>>             }
> >>>>>             if(inftime[i]==2){
> >>>>>                 loglh<-loglh + log(p[i])
> >>>>>             }
> >>>>>         }
> >>>>>     }
> >>>>>     return(loglh)
> >>>>> }
> >>>>> Here p[i] <- 1 - exp(-alpha*d^(-beta))
> >>>>> so,  log(p[i]) produces NaNs when exp(-alpha*d^(-beta)) is greater
> than
> >>>>> 1.
> >>>>> How can I remove it.After generating the out put we can omit it, but
> >>>>> the
> >>>>> problem is different.
> >>>>>
> >>>>>
> >>>
> >>> --
> >>> Sarah Goslee
> >>> http://www.functionaldiversity.org
> >>
> >>
> >
>  > ______________________________________________
> > 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]]

