[R] problem that arises after using the new version of "BRugs"

[moumita chatterjee]

Respected Sir,
                       With reference to my mail to you and the reply mail
by you dated 9th and 16th January, 2013, I am sending the reproducible code
in the attached document named " MODIFIED ANS ". I am also attaching the
txt file named "hazModel", which is required to save in my documents folder
to run the program. The file also contains the error message along with the
code.
                                          I would be very much obliged if
you kindly give me a solution for this question.
Thanking you.
                    Moumita Chatterjee
                    Research Scholar
                    University Of Calcutta.
-------------- next part --------------
x<-c(13,5,8,0,11)
delta<-c(1,1,0,1,0)


ZOSull<-
function (x, range.x, intKnots, drv = 0) 
{
    if (drv > 2) 
        stop("splines not smooth enough for more than 2 derivatives")
    library(splines)
    if (missing(range.x)) 
        range.x <- c(1.05 * min(x) - 0.05 * max(x), 1.05 * max(x) - 
            0.05 * min(x))
    if (missing(intKnots)) {
        numIntKnots <- min(length(unique(x)), 35)
        intKnots <- quantile(unique(x), seq(0, 1, length = (numIntKnots + 
            2))[-c(1, (numIntKnots + 2))])
    }
    numIntKnots <- length(intKnots)
    allKnots <- c(rep(range.x[1], 4), intKnots, rep(range.x[2], 
        4))
    K <- length(intKnots)
    L <- 3 * (K + 8)
    xtilde <- (rep(allKnots, each = 3)[-c(1, (L - 1), L)] + rep(allKnots, 
        each = 3)[-c(1, 2, L)])/2
    wts <- rep(diff(allKnots), each = 3) * rep(c(1, 4, 1)/6, 
        K + 7)
    Bdd <- spline.des(allKnots, xtilde, derivs = rep(2, length(xtilde)), 
        outer.ok = TRUE)$design
    Omega <- t(Bdd * wts) %*% Bdd
    eigOmega <- eigen(Omega)
    indsZ <- 1:(numIntKnots + 2)
    UZ <- eigOmega$vectors[, indsZ]
    LZ <- t(t(UZ)/sqrt(eigOmega$values[indsZ]))
    indsX <- (numIntKnots + 3):(numIntKnots + 4)
    UX <- eigOmega$vectors[, indsX]
    L <- cbind(UX, LZ)
    stabCheck <- t(crossprod(L, t(crossprod(L, Omega))))
    if (sum(stabCheck^2) > 1.0001 * (numIntKnots + 2)) 
        print("WARNING: NUMERICAL INSTABILITY ARISING\\\n              FROM SPECTRAL DECOMPOSITION")
    B <- spline.des(allKnots, x, derivs = rep(drv, length(x)), 
        outer.ok = TRUE)$design
    Z <- B %*% LZ
    attr(Z, "range.x") <- range.x
    attr(Z, "intKnots") <- intKnots
    return(Z)
}
 




BRugsMCMC<-
function (data, inits, parametersToSave, nBurnin, nIter, nThin, 
    modelFile) 
{
   
    if ((nBurnin < 100) | (nIter < 100)) 
        stop("currently only working for chains longer than 100")
    if ((100 * round(nBurnin/100) != nBurnin) | (100 * round(nIter/100) != 
        nIter)) 
        warning("chain lengths not multiples of 100; truncation may occur.")
    
    
    modelCheck(modelFile)
    
    modelData(bugsData(data))
    
    modelCompile(numChains = 1)
    
    modelInits(bugsInits(inits))
    
    burninIncSize <- round(nBurnin/(nThin * 100))
    for (i in 1:100) {
        modelUpdate(burninIncSize, thin = nThin)
        
    }
    samplesSet(parametersToSave)
    iterIncSize <- round(nIter/(nThin * 100))
    for (i in 1:100) {
        modelUpdate(iterIncSize, thin = nThin)
        
    }

    sims <- samplesStats("*")
    return(list(Stats = sims))
}











    require(BRugs)
    
        numKnots <- 25
     
        range.x <- c(0,202)
    ox <- order(x)
    x <- x[ox]
    delta <- delta[ox]
    n <- length(x)
    xTil <- as.numeric(names(table(x)))
    cnts <- as.numeric(table(x))
    nTil <- length(xTil)
    deltaTil <- rep(NA, nTil)
    for (iTil in 1:nTil) deltaTil[iTil] <- sum(delta[(1:n)[x == 
        xTil[iTil]]])
    sd.xTil <- sd(xTil)
    xTil <- xTil/sd.xTil
    range.x <- range.x/sd.xTil
    d1xTil <- xTil[2:nTil] - xTil[1:(nTil - 1)]
    d2xTil <- xTil[3:nTil] - xTil[1:(nTil - 2)]
    rcsrCnts <- rev(cumsum(rev(cnts)))
    trapLef <- 2 * cnts[1] * xTil[1] + (sum(cnts) - cnts[1]) * 
        (xTil[2] + xTil[1])
    trapMid <- cnts[2:(nTil - 1)] * d1xTil[-length(d1xTil)] + 
        rcsrCnts[1:(nTil - 2)] * d2xTil
    trapRig <- cnts[nTil] * d1xTil[length(d1xTil)]
    oTil <- log(c(trapLef, trapMid, trapRig)/2)
    X <- cbind(rep(1, nTil), xTil)
    numIntKnots <- 25
    intKnots <- quantile(xTil, seq(0, 1, length = numIntKnots + 
        2)[-c(1, numIntKnots + 2)])
    

Z <- ZOSull(xTil, intKnots = intKnots, range.x = range.x)
    numKnots <- ncol(Z)
    allData <- list(n = nTil, numKnots = numKnots, x = xTil, 
        delta = deltaTil, offs = oTil, Z = Z)
    parInits <- list(list(beta0 = 0.1, beta1 = 0.1, u = rep(0.1, 
        numKnots), tauU = 1))



#################################################
# the file hazModel.txt must be in the
# My Documents folder for the next step to run, I am attaching this with my mail.
#################################################



    BRugsMCMC(data = allData, inits = parInits, parametersToSave = c("beta0", 
        "beta1", "u", "tauU"), nBurnin =150, nIter =150, 
        nThin = 15, modelFile = "hazModel.txt")




betaMCMC <- rbind(samplesSample("beta0"), samplesSample("beta1"))


# An error is coming in this line. THe error window shows that BUGSHE~1.EXE has stopped working.# 

# the error message is #
model has probably not yet been updated
model has probably not yet been updated
Error in handleRes(res) : NA
In addition: Warning message:
running command '"C:/Users/MOU/Documents/R/win-library/2.15/BRugs/exec/BugsHelper.exe" "C:/Users/MOU/AppData/Local/Temp/RtmpwzJ97S" "C:/Users/MOU/AppData/Local/Temp/RtmpwzJ97S/trash" "file7f819516f3f.bug" "C:/Users/MOU/AppData/Local/Temp/RtmpwzJ97S/cmds.txt" "2"' had status 255 



-------------- next part --------------
model


{


   for(i in 1:n) 
 

  {


      log(mu[i]) <- beta0 + beta1*x[i] + inprod(u[],Z[i,]) + offs[i]
     

 delta[i] ~ dpois(mu[i])	
   

   }	


   for (k in 1:numKnots)


   {


      u[k] ~  dnorm(0,tauU)
  

 }



   beta0 ~ dnorm(0,1.0E-8)


   beta1 ~ dnorm(0,1.0E-8)
   
   
   tauU   ~ dgamma(0.01,0.01)
 

   sigU <- 1/sqrt(tauU)


}