[R] missing value where TRUE/FALSE needed ERROR

[Maram SAlem]

Hi all,

I'm trying to write a function to implement a Metropolis-within-Gibbs
algorithm for two parameters.I'm including a naive version here so as to be
able to spot the error I got. So I first generate the vectors, X and R,
that will help to start the algorithm using (for example):

n=8; m=5; p=0.1; t=0.9 ; JH=10;

R <- numeric(m)

W <- numeric(m)

V <- numeric(m)

U <- numeric(m)

X <- numeric(m)

Bay.alpha<- numeric (JH)

Bay.beta<- numeric (JH)

Bay.Surv <- numeric (JH)

hyp=c(3,15,6,22.5)

theta<-c(0.2,2)


alpha.curr<-theta[1]

beta.curr<- theta[2]



R[1]<-rbinom(1, n-m, p)

 for (i in 2:m-1) {

 R[i]<-rbinom(1,n-m-sum(R[1:i-1]),p)

 }

 R[m]<-n-m-sum(R[1:m-1])

W<-runif(m, min = 0, max = 1)

for (i in 1:m){

 V[i]<-W[i]^(1/(i+sum(R[(m-i+1):m])))

}

for (i in 1:m){

 U[i]<- 1- prod(V[(m-i+1):m])

}

for (i in 1:m){

X[i]<- ((-1/theta[1])*log(1-U[i]))^(1/theta[2])

}


Then, I defined three functions 1- alpha.update() for updating alpha (Gibbs
step) 2- bettarg(), for the target distribution of beta 3- beta.update()
for updating beta using a Metropolis Hastings technique.


   alpha.update=function(X, R, alpha.curr, beta.curr ,m, hyp)

  {

  o<-numeric(m)

     for (i in 1:m) {

       o[i]<- (1+R[i])*((X[i])^(beta.curr))

        }

   sh<-sum(o) + hyp[2] + (hyp[4]* beta.curr)

   rg<-rgamma(1, shape= m+hyp[1]+hyp[3] , rate = sh )

   return(rg)

   }



   bettarg<- function(X, R, alpha.curr, beta.curr ,m, hyp)

       {

           o<-numeric(m)

           for (i in 1:m) {

       o[i]<- (1+R[i])*((X[i])^( beta.curr))

        }

      bt<- beta.curr ^(m+hyp[3]-1) * prod((X)^( beta.curr -1)) *
exp(-1*alpha.curr *(sum(o) +      (hyp[4]* beta.curr)))

    return(bt)

       }

       beta.update<- function (X, R, alpha.curr, beta.curr ,m, hyp, cand.sd)

       {

       beta.cand<- rnorm(1, mean = beta.curr, sd = cand.sd)

      AB<- bettarg (X, R, alpha.curr, beta.curr = beta.cand ,m, hyp)

      CD<- bettarg (X, R, alpha.curr, beta.curr ,m, hyp)

      accept.prob <- AB/CD

     if (runif(1) <= accept.prob)

     beta.cand

    else  beta.curr

     }


Then I started a tiny chain using ten iterations but got this ERROR:


for (k in 1:JH)

{

alpha.curr<- alpha.update(X, R, alpha.curr, beta.curr ,m, hyp)

 Bay.alpha[k]<-alpha.curr

beta.curr<- beta.update (X, R, alpha.curr, beta.curr ,m, hyp, cand.sd=2)

     Bay.beta[k]<- beta.curr

     Bay.Surv [k]<- exp (-1*Bay.alpha[k] * (t)^Bay.beta[k])

     }


Error in if (runif(1) <= accept.prob) beta.cand else beta.curr :

  missing value where TRUE/FALSE needed

In addition: Warning message:

In rgamma(1, shape = m + hyp[1] + hyp[3], rate = sh) : NAs produced


I ran the code several times for only one iteration but everything was fine
with no errors or warnings, so I don't know from where does the missing
value/ NA come from?

In addition, I want to calculate the acceptance rate for the
Metropolis-Hastings step, is this possible?

Any help would be appreciated.
Thanks,

Maram Salem

	[[alternative HTML version deleted]]