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.pokha...@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@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.
>
______________________________________________
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.
