moumita chatterjee
2013-Jan-18 10:55 UTC
[R] problem that arises after using the new version of "BRugs"
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) }