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.
