Message-ID: <CAAmySGPzdXk3iE5t+4K+2OQ=iRt6s2udYovMH3B=fNtDc-Cd6A@mail.gmail.com>
Date: 2011-11-16T18:13:00Z
From: R. Michael Weylandt
Subject: Error in random walk Metroplis-hasting
In-Reply-To: <CAK=huh7UCFBbzYNtym8-u=gjMX9SzvoEXtXi5wbTTT4CLV-w_w@mail.gmail.com>
Three things:
1) This isn't reproducible without your data file. Work out a minimal
reproducible example -- I bet you'll find your answer along the way --
and supply it.
2) What do the warnings say: they are usually pretty good at directing
you to trouble.
3) Don't use attach. Seriously -- just load the data in once and pass
it around where needed. (It's going to slow you down to re-load it
each time (~400 times) you call likelihood.)
Michael
PS -- As a style point, might I suggest you put more spaces around the
assignment operator "<-" it's hard (for both a human and the R
interpreter) to tell is x<-3 means to set x equal to 3 or to test if x
is less than "-3" and sometimes this leads to very tricky bugs.
On Wed, Nov 16, 2011 at 10:36 AM, Gyanendra Pokharel
<gyanendra.pokharel at gmail.com> wrote:
> Hi R community,
> I have some data set and construct the likelihood as follows
> likelihood <- function(alpha,beta){
> ? ?lh<-1
> ? ?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){
> ? ? ? ? ? ? ? ?lh<-lh*(1-p[i])
> ? ? ? ? ? ?}
> ? ? ? ? ? ?if(inftime[i]==2){
> ? ? ? ? ? ? ? ?lh<-lh*p[i]
> ? ? ? ? ? ?}
> ? ? ? ?}
> ? ?}
> ? ?return(lh)
> }
> ?Then, I want to simulate this by using random walk Metropolis- Hasting
> algorithm with the single parameter update. i have the following R function
> mh.epidemic <-function(m,alpha, beta){
> ? ? ?z<- array(0,c(m+1, 2))
> ? ?z[1,1] <- alpha
> ? ?z[1,2] <- beta
> ? ?for(t in 2:m){
> ? ? ? ?u <- runif(1)
> ? ? ? ?y1 <- rexp(1,z[t-1,1])
> ? ? ? ?y2 <-rexp(1,z[t-1,2])
> ? ? ? ?z[t,1] <- likelihood(y1,z[t-1,2])/likelihood(z[t-1,1],z[t-1,2])
> ? ? ? ?a1 <-min(1,z)
> ? ? ? ?if(u<=a1) z[t,] <- y1 else
> ? ? ? ?z[t,2] <-z[t-1,1]
>
> ? ? ? ?z[t,2]<-likelihood(z[t,1], y2)/likelihood(z[t,1],z[t-1,2])
> ? ? ? ?accept <-min(1,z)
> ? ? ? ?if(u<=accept) z[t,] <- y2 else
> ? ? ? ?z[t,2] <-z[t-1,1]
> ? ?}
> ? ?z
> }
> mh.epidemic(100, 1,2)
> when I run it shows the following error:
> Error in if (u <= accept) z[t, ] <- y2 else z[t, 2] <- z[t - 1, 1] :
> ?missing value where TRUE/FALSE needed
> I know it is due to the NaN produce some where. So I tried by running the
> code separately, as follows
> m<- 100
>> ? ? alpha <-1
>> ? ? beta <- 2
>> ? ? z<- array(0,c(m+1, 2))
>> ? ? z[1,1] <- alpha
>> ? ? z[1,2] <- beta
>> ? ? for(t in 2:m){
> + ? ? ? ? u <- runif(1)
> + ? ? ? ? y1 <- rexp(1,z[t-1,1])
> + ? ? ? ? y2 <-rexp(1,z[t-1,2])
> + ? ? ? ? z[t,1] <- likelihood(y1,z[t-1,2])/likelihood(z[t-1,1],z[t-1,2])
> + ? ? ? ? a1 <-min(1,z)
> + }
> There were 50 or more warnings (use warnings() to see the first 50)
>> y1
> [1] NaN
>> y2
> [1] NaN
> why, this simulation from exponential function is produce NaN?
> If some one help me it will be great.
>
> ? ? ? ?[[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.
>