1

I am currently working on an OpenBUGS code regarding bivariate normal distribution. The distribution uses a Wishart prior for precision, and I am having some troubling when updating the model.

My model will load, and it does compile with my data. However, when updating, I receive an error message saying 'TRAP 126 (not yet implemented)', and it is regarding tau.org[1,1].

I think the problem lies with Wishart prior but I am not sure how to fix it. I have tried a simpler model with the same prior and there was no problem updating, and I realised the problem may be caused by eps_theta.

I have checked most sources I could, but the error persist. It would be great if someone could shade some light on this problem. Thanking in advance! 

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)))
Oliver
  • 11
  • 1

0 Answers0