Как переписать значение наблюдаемого узла? - PullRequest
Новая
       89

Как переписать значение наблюдаемого узла?

1 голос
/ 10 ноября 2019

Новичок в JAGS и работа с моделью, указанной в этой статье . Похоже, что авторы использовали предыдущую версию JAGS, которая была ближе к BUGS, поскольку в какой-то момент это появляется в блоке модели (строки 28-29): 

  z[i,1:(K-1)] ~ dmnorm(etas[i,], omega)
  z[i,K] <- 0

JAGS жалуется на следующую ошибку: 

Error in setParameters(init.values[[i]], i) : Error in node z[1,4]
Cannot overwrite value of observed node

Проверка руководства JAGS, ошибка очевидна. В разделе A.4 Преобразования данных говорится, что BUGS позволяет дважды поместить стохастический узел в левую часть отношения. В качестве обходного пути предлагается включить преобразование данных в отдельный блок data. Это все еще терпит неудачу.

Кто-нибудь пытался повторить эту работу и нашел успех? Есть подсказки?

Полная модель JAGS ниже. Подозреваемые назначения: z[i,K] <- 0 и beta[j,K] <- 0

data {
   for (j in 1:(K-1)) {
      for (i in 1:N) {
         ones[i,j] <- 1
      }
      for (k in 1:K) {
         C[j,k] <- equals(j,k) - equals(k,K)
      }   
   }
   R <- (K-1) * C %*% t(C)
   lower <- -1e+3
   upper <-  1e+3
}
model { 
   for (i in 1:N) {
      for (k in 1:(K-1)) {
         bounds[i,k,1] <- equals(y[i,k],K)*lower + inprod(z[i,], equals(y[i,],y[i,k] + 1))
         bounds[i,k,2] <- equals(y[i,k],1)*upper + inprod(z[i,], equals(y[i,],y[i,k] - 1))
         ones[i,k] ~ dinterval(z[i,k], bounds[i,k,])
         etas[i,k] <- inprod(x[i,], beta[,k])
      }
      z[i,1:(K-1)] ~ dmnorm(etas[i,], omega)
      z[i,K] <- 0
   }
   for (j in 1:(P+1)) {
      for (k in 1:(K-1)) {
         beta[j,k] ~ dnorm(0,1e-3)
      }
      beta[j,K] <- 0
      for (k in 1:K) {   
          for (t in 1:K) {
             beta.identified[j,k,t] <- (beta[j,k] - beta[j,t])*sqrt(const)
          }
      }
   }
   omega ~ dwish(R, K-1)
   sigma <- inverse(omega)
   const <- pow(K/exp(logdet(sigma)), 1/K)
   sigma.identified <- sigma*const
}

1 Ответ

1 голос
/ 13 ноября 2019

Я обратился к автору статьи и признал, что последующий выпуск JAGS нарушил оригинальный код. По его словам, решение состоит в том, что ...

"Но есть относительно простой обходной путь, который включает использование дополнительной переменной, которая эффективно представляет различия К-1 в скрытых ответах, а затем сопоставляет ее сПеременная Z. Я приложу файл, который содержит обновленные сценарии из этого документа. Это требует изменения способа создания начальных значений для алгоритма. Функция makeinits (в файле data.R) былатакже обновлено. "

Ниже представлен код новой модели: 

data {
   for (j in 1:(K-1)) {
      for (i in 1:N) {
         ones[i,j] <- 1
      }
      for (k in 1:K) {
         C[j,k] <- equals(j,k) - equals(k,K)
      }   
   }
   R <- (K-1) * C %*% t(C)
   lower <- -1e+3
   upper <-  1e+3
}
model { 
   for (i in 1:N) {
      for (k in 1:(K-1)) {
         bounds[i,k,1] <- equals(y[i,k],K)*lower + inprod(z[i,], equals(y[i,],y[i,k] + 1))
         bounds[i,k,2] <- equals(y[i,k],1)*upper + inprod(z[i,], equals(y[i,],y[i,k] - 1))
         ones[i,k] ~ dinterval(z[i,k], bounds[i,k,])
         etas[i,k] <- inprod(x[i,], beta[,k])
      }
      u[i,1:(K-1)] ~ dmnorm(etas[i,], omega)
      z[i,1:(K-1)] <- u[i,1:(K-1)]
      z[i,K] <- 0
   }
   for (j in 1:(P+1)) {
      for (k in 1:(K-1)) {
         beta[j,k] ~ dnorm(0,1e-3)
      }
      beta[j,K] <- 0
      for (k in 1:K) {   
          for (t in 1:K) {
             beta.identified[j,k,t] <- (beta[j,k] - beta[j,t])*sqrt(const)
          }
      }
   }
   omega ~ dwish(R, K-1)
   sigma <- inverse(omega)
   const <- pow(K/exp(logdet(sigma)), 1/K)
   sigma.identified <- sigma*const
}

с новой функцией makeinits: 

makeinits <- function (y) {
    n <- nrow(y)
    m <- ncol(y)
    z <- matrix(NA, n, m)
    for (i in 1:n) { 
        z[i,] <- sort(rnorm(m))[rank(-y[i,])]
        z[i,] <- z[i,] - z[i,m]
        }
    z[,-ncol(y)]
}
...