How do I specify for a variable to be uniformly sampled from the discrete values $1,2,\ldots,10$ in JAGS? I guess I want to use the dcat() distribution, but I have yet to find good documentation of the parameters for this.
Solved – Categorical distribution in JAGS
jags
Related Solutions
You can avoid this problem altogether by sampling from the untruncated distribution of Y[k], then (in R) discarding all samples for which Y[k] doesn't lie within the constraint bounds. This is a perfectly valid operation, however, if you have few posterior observations in the feasible region, you'll naturally have a large simulation error associated with your posterior distribution(s).
As a side note, you might want to avoid the Gamma(0.01,0.01) prior on the variance; see for example this presentation by Andrew Gelman and this paper, also by Gelman, for reasons why and alternative suggestions.
Here is one example of implementing a basic macro substitution system for JAGS scripts.
Explanation of the system
- Define a function that takes as arguments any optional elements of the script.
- For any aspects of the script that vary across argument values, record a macro token. This should be some unique text. Starting and ending with some symbols may assist to make this unambiguous.
- For each macro token include code for what the replacement string should be under all possible combinations of the function arguments.
- Replace the macro tokens with appropriate placement text based on arguments.
- The code below provides one example of how macros could be specified and how to apply the macros to the raw script (suggestions for doing this more elegantly are welcome).
- The function returns a replaced model string that could if necessary be passed to a
textConnectionfunction for use with rjags.
I like this system for a few reasons:
- the raw script is easy to read
- the resulting script has proper indentation
- you don't have to worry about errors related to commands appearing on the same line
Example
Specifically, this example, aims to allow the user to fit a particular type of multilevel nonlinear model. It is designed to allow for three functional forms: a two parameter power, three parameter power, and a three parameter exponential.
The macro section structures macros as a list of lists. The top level list contains one element for each macro token. For each macro token, there is the macro token text, and the conditional replacement text.
Finally, a for loop applies all the macro replacements to the raw script.
See below (scrolling is required):
jags_model <- function (f=c('power2', 'power3', 'exp3')) {
f <- match.arg(f)
# raw script
script <-
"model {
# Model
for (i in 1:length(y)) {
$FUNCTION
y[i] ~ dnorm(mu[i], tau[subject[i]])
}
# Random coefficients
for (i in 1:N) {
theta1[i] ~ dnorm(theta1.mu, theta1.tau)T(0, 1000)
theta2[i] ~ dnorm(theta2.mu, theta2.tau)T(-10, 0)
$THETA3DISTRIBUTION
sigma[i] ~ dnorm(sigma.mu, sigma.tau)T(0, 100);
tau[i] <- 1/(sigma[i]^2)
}
theta1.mu ~ dunif(0, 100)
theta2.mu ~ dunif(-2, 0)
$THETA3PRIOR.MU
sigma.mu ~ dunif(0, 20)
theta1.sigma ~ dunif(0, 100)
theta2.sigma ~ dunif(0, 2)
$THETA3PRIOR.SIGMA
sigma.sigma ~ dunif(0, 10)
# Transformations
theta1.tau <- 1/(theta1.sigma^2)
theta2.tau <- 1/(theta2.sigma^2)
$THETA3.TAU
sigma.tau <- 1/(sigma.sigma^2)
}"
# define macros
macros <- list(list("$FUNCTION",
switch(f,
power2="mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]])",
power3="mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]]) + theta3[subject[i]];",
exp3="mu[i] <- theta1[subject[i]] * exp(theta2[subject[i]] * (trial[i] - 1)) + theta3[subject[i]];") ),
list("$THETA3DISTRIBUTION",
switch(f,
power3=, exp3= "theta3[i] ~ dnorm(theta3.mu, theta3.tau)T(0, 1000)",
power2="") ),
list("$THETA3PRIOR.MU",
switch(f,
power3=, exp3= "theta3.mu ~ dunif(0, 100)",
power2="") ),
list("$THETA3PRIOR.SIGMA",
switch(f,
power3=, exp3= "theta3.sigma ~ dunif(0, 100)",
power2="") ),
list("$THETA3.TAU",
switch(f,
power3=, exp3= "theta3.tau <- 1/(theta3.sigma^2)",
power2="") )
)
# apply macros
for (m in seq(macros)) {
script <- gsub(macros[[m]][1], macros[[m]][2], script, fixed=TRUE)
}
script
}
Demonstration
Thus, we can produce the processed JAGS model with
x <- jags_model(f='power3')
And if we want to view the resulting model, we can do
cat(x)
which results in
model {
# Model
for (i in 1:length(y)) {
mu[i] <- theta1[subject[i]] * pow(trial[i], theta2[subject[i]]) + theta3[subject[i]];
y[i] ~ dnorm(mu[i], tau[subject[i]])
}
# Random coefficients
for (i in 1:N) {
theta1[i] ~ dnorm(theta1.mu, theta1.tau)T(0, 1000)
theta2[i] ~ dnorm(theta2.mu, theta2.tau)T(-10, 0)
theta3[i] ~ dnorm(theta3.mu, theta3.tau)T(0, 1000)
sigma[i] ~ dnorm(sigma.mu, sigma.tau)T(0, 100);
tau[i] <- 1/(sigma[i]^2)
}
theta1.mu ~ dunif(0, 100)
theta2.mu ~ dunif(-2, 0)
theta3.mu ~ dunif(0, 100)
sigma.mu ~ dunif(0, 20)
theta1.sigma ~ dunif(0, 100)
theta2.sigma ~ dunif(0, 2)
theta3.sigma ~ dunif(0, 100)
sigma.sigma ~ dunif(0, 10)
# Transformations
theta1.tau <- 1/(theta1.sigma^2)
theta2.tau <- 1/(theta2.sigma^2)
theta3.tau <- 1/(theta3.sigma^2)
sigma.tau <- 1/(sigma.sigma^2)
}
Best Answer
Have a look at the JAGS user manual.
dcat(pi)is defined with a density of $\pi_x/\sum_i{\pi_i}$.