R help - Dec 2007 - Metropolis-Hastings within Gibbs coding error

Dear list,

After running for a while, it crashes and gives the following error message: can
anybody suggest how to deal with this?

Error in if (ratio0[i] < log(runif(1))) { : 
  missing value where TRUE/FALSE needed



################### original program ########
p2 <- function (Nsim=1000){

x<- c(0.301,0,-0.301,-0.602,-0.903,-1.208, -1.309,-1.807,-2.108,-2.71) #
logdose
n<-c(19,20,19,21,19,20,16,19,40,81) # total subject in dose-response
experiment
y<-c(19,18,19,14,15,4,0,0,0,2) # success in each trials
dta<-cbind(x,n,y)
dta<-as.data.frame(dta) # creating data frame

proposal.b0 = current.b0 = ratio0 = double(Nsim) # blank vector
proposal.b1 = current.b1 = ratio1 = double(Nsim) # blank vector
index <- 1:Nsim # creating index
a0 <- 10 # initial value (assumed) for tau
b0 <- (.01) # initial value (assumed) for tau
fit <- glm((y/n)~x,family=binomial, weight = n, data=dta) # initial value for
beta

parameters <- c("Beta0", "Beta1", "Tau")
parameter.matrix <- array(NA,c(Nsim,3)) # blank array
parameter.matrix <- as.data.frame(parameter.matrix) # creating data frame
parameter.matrix[1,] <- c(fit$coef[1],fit$coef[2],rgamma(1, shape=a0, scale =
b0)) # putting initial values

for (i in 2:Nsim){
# generating Gibbs sampler

parameter.matrix[i,]<- c(rnorm(1, 0, (1/parameter.matrix[i-1,3])), rnorm(1,
0, (1/parameter.matrix[i-1,3])), rgamma(1, shape=(a0+1),
rate=(1/b0+(parameter.matrix[i-1,1]^2+parameter.matrix[i-1,2]^2)/2)))
# as the Gamma with specified parameters is the conditional for tau given beta,
data

# implementing Metropolis-Hastings within Gibbs to get estimates of beta0 and
beta1

proposal.b0[i]<-sum(log(
((exp(parameter.matrix[i,1])^y)/((1+exp(parameter.matrix[i,1])^n))*exp(-parameter.matrix[i-1,3]*(parameter.matrix[i,1]^2)/2))))
proposal.b1[i]<-sum(log( ((exp(parameter.matrix[i,2]*x)^y) /
((1+exp(parameter.matrix[i,2]*x))^n) *
exp(-parameter.matrix[i-1,3]*(parameter.matrix[i,2]^2)/2) )))
# as the above is the conditional for beta's given tau, data

# took logarithm to take care of the 0 problem in product space, but does not
help much

current.b0[i]<-sum(log( ((
exp(parameter.matrix[i-1,1])^y)/((1+exp(parameter.matrix[i-1,1])^n))*exp(-parameter.matrix[i-1,3]*(parameter.matrix[i-1,1]^2)/2))))
current.b1[i]<-sum(log( (( exp(parameter.matrix[i-1,2]*x)^y) /
((1+exp(parameter.matrix[i-1,2]*x))^n) *
exp(-parameter.matrix[i-1,3]*(parameter.matrix[i-1,2]^2)/2))))

# ratio0 id for beta0
if(current.b0[i]==0) {ratio0[i]=1} else {ratio0[i] <-
proposal.b0[i]-current.b0[i]}
if (ratio0[i] < log(runif(1))) {parameter.matrix[i,1] <-
parameter.matrix[i-1,1]}
# for beta0

# ratio1 id for beta1
if(current.b1[i]==0) {ratio1[i]=1} else {ratio1[i]=proposal.b1[i]-current.b1[i]}
if (ratio1[i] < log(runif(1))) {parameter.matrix[i,2] <-
parameter.matrix[i-1,2]}
# for beta1
cat("At Iteration ", i, "ratio0 and ratio1 are", ratio0[i],
ratio1[i], "\n" )
}

x11()
plot(parameter.matrix[,1], parameter.matrix[,2], type="b",
xlab="beta.0", ylab="beta.1")
write.table(parameter.matrix, file="z:\\paramaters.txt", quote = F,
sep = " ")
x11()

par(mfrow=c(3,1))
plot(index, parameter.matrix[index,1], type="l",
xlab="Index", ylab="beta0")
plot(index, parameter.matrix[index,2], type="l",
xlab="Index", ylab="beta1")
plot(index, parameter.matrix[index,3], type="l",
xlab="Index", ylab="tau")
}

p2(Nsim=1000)





     
____________________________________________________________________________________


-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: mhwg.txt
Url:
https://stat.ethz.ch/pipermail/r-help/attachments/20071204/78e50b58/attachment.txt

On 12/4/2007 11:06 AM, Mohammad Ehsanul Karim wrote:
> Dear list,
> 
> After running for a while, it crashes and gives the following error
message: can anybody suggest how to deal with this?
> 
> Error in if (ratio0[i] < log(runif(1))) { : 
>   missing value where TRUE/FALSE needed
Use options(error=recover) to drop to the debugger when the error 
occurs, and you could examine the value of ratio0.  runif(1) should not 
give a 0, so log(runif(1)) should succeed.

Duncan Murdoch
> 
> 
> 
> ################### original program ########
> p2 <- function (Nsim=1000){
> 
> x<- c(0.301,0,-0.301,-0.602,-0.903,-1.208, -1.309,-1.807,-2.108,-2.71) #
logdose
> n<-c(19,20,19,21,19,20,16,19,40,81) # total subject in dose-response
experiment
> y<-c(19,18,19,14,15,4,0,0,0,2) # success in each trials
> dta<-cbind(x,n,y)
> dta<-as.data.frame(dta) # creating data frame
> 
> proposal.b0 = current.b0 = ratio0 = double(Nsim) # blank vector
> proposal.b1 = current.b1 = ratio1 = double(Nsim) # blank vector
> index <- 1:Nsim # creating index
> a0 <- 10 # initial value (assumed) for tau
> b0 <- (.01) # initial value (assumed) for tau
> fit <- glm((y/n)~x,family=binomial, weight = n, data=dta) # initial
value for beta
> 
> parameters <- c("Beta0", "Beta1", "Tau")
> parameter.matrix <- array(NA,c(Nsim,3)) # blank array
> parameter.matrix <- as.data.frame(parameter.matrix) # creating data
frame
> parameter.matrix[1,] <- c(fit$coef[1],fit$coef[2],rgamma(1, shape=a0,
scale = b0)) # putting initial values
> 
> for (i in 2:Nsim){
> # generating Gibbs sampler
> 
> parameter.matrix[i,]<- c(rnorm(1, 0, (1/parameter.matrix[i-1,3])),
rnorm(1, 0, (1/parameter.matrix[i-1,3])), rgamma(1, shape=(a0+1),
rate=(1/b0+(parameter.matrix[i-1,1]^2+parameter.matrix[i-1,2]^2)/2)))
> # as the Gamma with specified parameters is the conditional for tau given
beta, data
> 
> # implementing Metropolis-Hastings within Gibbs to get estimates of beta0
and beta1
> 
> proposal.b0[i]<-sum(log(
((exp(parameter.matrix[i,1])^y)/((1+exp(parameter.matrix[i,1])^n))*exp(-parameter.matrix[i-1,3]*(parameter.matrix[i,1]^2)/2))))
> proposal.b1[i]<-sum(log( ((exp(parameter.matrix[i,2]*x)^y) /
((1+exp(parameter.matrix[i,2]*x))^n) *
exp(-parameter.matrix[i-1,3]*(parameter.matrix[i,2]^2)/2) )))
> # as the above is the conditional for beta's given tau, data
> 
> # took logarithm to take care of the 0 problem in product space, but does
not help much
> 
> current.b0[i]<-sum(log( ((
exp(parameter.matrix[i-1,1])^y)/((1+exp(parameter.matrix[i-1,1])^n))*exp(-parameter.matrix[i-1,3]*(parameter.matrix[i-1,1]^2)/2))))
> current.b1[i]<-sum(log( (( exp(parameter.matrix[i-1,2]*x)^y) /
((1+exp(parameter.matrix[i-1,2]*x))^n) *
exp(-parameter.matrix[i-1,3]*(parameter.matrix[i-1,2]^2)/2))))
> 
> # ratio0 id for beta0
> if(current.b0[i]==0) {ratio0[i]=1} else {ratio0[i] <-
proposal.b0[i]-current.b0[i]}
> if (ratio0[i] < log(runif(1))) {parameter.matrix[i,1] <-
parameter.matrix[i-1,1]}
> # for beta0
> 
> # ratio1 id for beta1
> if(current.b1[i]==0) {ratio1[i]=1} else
{ratio1[i]=proposal.b1[i]-current.b1[i]}
> if (ratio1[i] < log(runif(1))) {parameter.matrix[i,2] <-
parameter.matrix[i-1,2]}
> # for beta1
> cat("At Iteration ", i, "ratio0 and ratio1 are",
ratio0[i], ratio1[i], "\n" )
> }
> 
> x11()
> plot(parameter.matrix[,1], parameter.matrix[,2], type="b",
xlab="beta.0", ylab="beta.1")
> write.table(parameter.matrix, file="z:\\paramaters.txt", quote =
F, sep = " ")
> x11()
> 
> par(mfrow=c(3,1))
> plot(index, parameter.matrix[index,1], type="l",
xlab="Index", ylab="beta0")
> plot(index, parameter.matrix[index,2], type="l",
xlab="Index", ylab="beta1")
> plot(index, parameter.matrix[index,3], type="l",
xlab="Index", ylab="tau")
> }
> 
> p2(Nsim=1000)
> 
> 
> 
> 
> 
>      
____________________________________________________________________________________
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 
> ______________________________________________
> 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.