4

我正在尝试使用 R2WinBugs 和 BRugs 通过 R 运行贝叶斯模型，但遇到了我无法解决的错误。我检查了我的模型在语法上是否正确，这是我尝试运行时的输出（下面的示例）：

model is syntactically correct
data loaded
model compiled
Initializing chain 1: 
model is initialized
model is already initialized
Sampling has been started ...
Error in handleRes(res) : NA
In addition: Warning message:
running command '"C:/Users/user/Documents/R/win-library/3.0/BRugs/exec/BugsHelper.exe" "C:/Users/user/AppData/Local/Temp/RtmpWweqIM" "C:/Users/user/AppData/Local/Temp/RtmpWweqIM/trash" "file32f44e4c5ab7.bug" "C:/Users/user/AppData/Local/Temp/RtmpWweqIM/cmds.txt" "1"' had status 144

会话信息：

> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-w64-mingw32/x64 (64-bit)

RStudio 0.99.441 Windows 7

这是一个示例，虽然我不确定如何附加数据，如果您想要完整的数据集，我可以在这里复制它，请告诉我。

library(BRugs)
library(R2WinBUGS)
library(data.table)
> dat
       date mth inc ind  response
  1: Jan-10   1   1   0  100.0000
  2: Jan-10   1   2   0  103.1304
  3: Jan-10   1   3   0  106.2609
  4: Jan-10   1   4   0  109.3913
  5: Jan-10   1   5   0  112.5217
 ---                             
572: Dec-10  12  44   1 1887.4783
573: Dec-10  12  45   1 1890.6087
574: Dec-10  12  46   1 1893.7391
575: Dec-10  12  47   1 1896.8696
576: Dec-10  12  48   1 1900.0000

# set data values as individual vectors
N <- nrow(dat)
periods <- length(unique(dat$date))
real_response <- dat$response
TT <- dat$mth
P <- dat$inc

ind <- dat$ind

x_dat <- list("real_response"=real_response,"P"=P,"N"=N,"TT"=TT,"periods"=periods, "ind"=ind)

CreateModel <- function(){

    for(k in 1:N) {
        P2[k] <- P[k]*30+ind[k]*speed
    }

    for(i in 1:N) {

        real_response[i] ~ dnorm(real_responseHat[i],t_response[i])
        t_response[i]  <- pow(var_response[i],-1)
        var_response[i] <- pow(sigma,2)* mu_target * (1 - exp(-pow((P[i]/theta),
                                                                 omega[TT[i]])))
        real_responseHat[i] <- TARGET[TT[i]] * (1 - exp(-pow((P2[i]/theta),
                                                             omega[TT[i]])))
    }

    sigma ~ dgamma(0.5,0.001)
    theta ~ dgamma(0.5,0.001)

    for(j in 1:periods){
        newmu[j] <- mu_target*pow((1+adj_factor), j)
        TARGET[j] ~ dnorm(newmu[j],t_target)
        omega[j] ~ dgamma(s1,s2)
    }

    speed ~ dcat(p[])
    for(j in 1:15){
        p[j] <- 1/15
    }

    adj_factor ~ dnorm(2, 1.2)
    t_target <- pow(s_target,-2)
    s_target ~ dgamma(0.5,0.0001)
    s1 ~ dgamma(1.5,.5)
    s2 ~ dgamma(1.5,.5)
    mu_target ~ dnorm(1000,0.0001)
}

mod <- file.path(paste(getwd(), "\\stackex_model.txt", sep=""))
writeModel(CreateModel, mod)
modelCheck(mod)
# file.show(mod)

##

inits <- function() {
    list("TARGET" = rep(1000,periods), "sigma" = 5, "theta" = 15, "omega" = rep(0.5,periods),"mu_target"=1000,"s_target"=60, "s1"=1.5, "s2"=1, "speed"=3, "adj_factor"=2)
}

parms <- c("TARGET", "sigma", "theta", "omega","mu_target","s_target","s1","s2", "speed", "adj_factor")

system.time(bmod <- bugs(data=x_dat, 
                         inits=inits, parameters.to.save=parms,
                         n.iter=100, n.chains=1, n.burnin=50, 
                         n.thin=1, model.file=mod, 
                         debug=TRUE, codaPkg=FALSE,
                         bugs.directory="C:/Program Files (x86)/OpenBUGS/OpenBUGS323",
                         program="OpenBUGS"))

我正在继续尝试调试这个问题，但如果有人以前见过这个或者对正在发生的事情有线索，我会非常感激。

此外，迭代和细化设置非常低，因此调试更容易。一旦我知道模型运行，我计划增加这些（并希望并行化）。

4

0 回答 0