1

我目前正在研究有关二元正态分布的 OpenBUGS 代码。该分布使用 Wishart 先验来提高精度，在更新模型时我遇到了一些麻烦。

我的模型将加载，它会使用我的数据进行编译。但是，在更新时，我收到一条错误消息，上面写着“TRAP 126（尚未实现）”，它与 tau.org[1,1] 有关。

我认为问题出在 Wishart 之前，但我不知道如何解决它。我尝试了一个更简单的模型，之前的更新没有问题，我意识到问题可能是由 eps_theta 引起的。

我已经检查了大多数来源，但错误仍然存​​在。如果有人能对这个问题有所了解，那就太好了。提前致谢！ 

model { 
for (i in 1:10) {
    y[i,1:2] ~ dmnorm(loc[i,1:2] , tau.t[i, ,])

    mu[i,1] <- alpha[1] + beta[2]*x[i]   #Eliminate 'j' due to theta 
    mu[i,2] <- alpha[2] + beta[2]*x[i]


    skew[i,1] <- lambda[i]*eps_theta[1]
    skew[i,2] <- lambda[i]*eps_theta[2]


    loc[i,1] <- mu[i,1] + skew[i,1]
    loc[i,2] <- mu[i,2] + skew[i,2]
    }  #close i loop 

    eps_theta[1] <- cov.org[1,1]*theta[1] + cov.org[1,2]*theta[2]
    eps_theta[2] <- cov.org[2,1]*theta[1] + cov.org[2,2]*theta[2]

    # tau.t
    for (i in 1:10){
        tau.t[i,1,1] <- inv.lambda[i]*tau.org[1,1] 
        tau.t[i,1,2] <- inv.lambda[i]*tau.org[1,2] 
        tau.t[i,2,1] <- inv.lambda[i]*tau.org[2,1] 
        tau.t[i,2,2] <- inv.lambda[i]*tau.org[2,2] 
        inv.lambda[i] ~ dgamma(dfby2, dfby2)
        lambda[i] <- 1/inv.lambda[i]
    }

    cov.org[1:2,1:2] <- inverse(tau.org[ , ])

    for (i in 1:10){     
        cov.t[i,1:2,1:2] <- inverse(tau.t[i, ,]) 
    } #Think if this is necessary 

    #Priors 
    for (j in 1:2){
        alpha[j] ~ dnorm(0, 1)
        beta[j] ~ dnorm(1, 1)
        theta[j] ~ dunif(-30,30) 
    }   

    tau.org[1:2,1:2] ~ dwish(R[,],2)  #R value provided in data
    dfby2 <- df/2
    df ~ dexp(0.1)I(1.01,40)


}




# Initials  
list( 
alpha = c(0,0), 
beta = c(1,1), 
df = 15, 
theta = c(0,0),
tau.org = structure(
    .Data = c(1,0.1,0.1,1),
    .Dim=c(2,2)))   

list( 
alpha = c(-5,-5), 
beta = c(5,5), 
df = 50, 
theta = c(5,5),
tau.org = structure(
    .Data = c(0.1,0,0,0.1),
    .Dim=c(2,2)))   

list( 
alpha = c(2,2), 
beta = c(-1,-1), 
df = 30, 
theta = c(-5,-5),
tau.org = structure(
    .Data = c(10,1,1,10),
    .Dim=c(2,2)))   



#Provide data
list(x=c(-2.6378, 0.7446, 0.5019, -0.9829, 1.2866, -0.2914, 0.5272, -4.9235, -1.5481, 0.1563),
y = structure(
.Data = c(-3.2529,  -3.4632,
                -1.4764,    -0.2982,
                2.8630, 0.7316,
                -2.1607,    -0.6788,
                0.6701, 1.0510,
                -1.3155,    0.6880,
                1.0596, 0.0861,
                -3.6570,    -3.8032,
                -1.6598,    -2.7763,
                3.9920, 0.3544), 
    .Dim = c(10,2)),

    R = structure(
        .Data=c(0.001,0,0,0.001),
        .Dim=c(2,2)))
4

0 回答 0