Je suis déconcerté par une simple question dans R JAGS. J'ai par exemple 10 paramètres: d [1], d [2], ..., d [10]. Il est intuitif d'après les données qu'ils devraient augmenter. Je veux donc leur imposer une contrainte.

Voici ce que j'ai essayé de faire mais cela donne des messages d'erreur disant "Nœud incompatible avec les parents":

    model{
    ...
    for (j in 1:10){
    d.star[j]~dnorm(0,0.0001)
    }
    d=sort(d.star)
    }

Ensuite, j'ai essayé ceci:

  d[1]~dnorm(0,0.0001)
  for (j in 2:10){
   d[j]~dnorm(0,0.0001)I(d[j-1],)
  }

Cela a fonctionné, mais je ne sais pas si c'est la bonne façon de procéder. Pouvez-vous partager vos réflexions?

Merci!

1
user3669725 5 avril 2017 à 16:42

2 réponses

Meilleure réponse

Si vous n'êtes jamais sûr de quelque chose comme ça, il est préférable de simplement simuler quelques données pour déterminer si la structure du modèle que vous suggérez fonctionne (alerte spoiler: c'est le cas).

Voici le modèle que j'ai utilisé:

cat('model{
  d[1] ~ dnorm(0, 0.0001) # intercept
  d[2] ~ dnorm(0, 0.0001)
  for(j in 3:11){
    d[j] ~ dnorm(0, 0.0001) I(d[j-1],)
  }
  for(i in 1:200){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- inprod(d, x[i,])
  }
  tau ~ dgamma(0.01,0.01)
  }',
file = "model_example.R")```

Et voici les données que j'ai simulées pour les utiliser avec ce modèle.

library(run.jags)
library(mcmcplots)

# intercept with sorted betas 
set.seed(161)
betas <- c(1,sort(runif(10, -5,5)))

# make covariates, 1 for intercept
x <- cbind(1,matrix(rnorm(2000), nrow = 200, ncol = 10))

# deterministic part of model
y_det <- x %*% betas

# add noise
y <- rnorm(length(y_det), y_det, 1)

data_list <- list(y = as.numeric(y), x = x)

# fit the model
mout <- run.jags('model_example.R',monitor = c("d", "tau"), data = data_list)

Ensuite, nous pouvons tracer les estimations et superposer les vraies valeurs des paramètres

caterplot(mout, "d", reorder = FALSE)
points(rev(c(1:11)) ~ betas, pch = 18,cex = 0.9)

Les points noirs sont les vraies valeurs des paramètres, les points et les lignes bleus sont les estimations. Il semble que cette configuration fonctionne bien tant qu'il y a suffisamment de données pour estimer tous ces paramètres. entrez la description de l'image ici

1
mfidino 6 avril 2017 à 14:47

Il semble qu'il y ait une erreur de syntaxe dans la première implémentation. Essayez simplement:

model{
  ...
  for (j in 1:10){
    d.star[j]~dnorm(0,0.0001)
  }
  d[1:10] <- sort(d.star)  # notice d is indexed.
}

Et comparez les résultats avec ceux de la deuxième implémentation. D'après la documentation, ils sont tous les deux corrects, mais il est conseillé d'utiliser la fonction sort.

1
altroware 6 avril 2017 à 16:47