model {
 
# Following loop defines the base MPT model
 
for (n in 1:subjs) {   # Loop iterates through, subjs, subjects
 
    # theta[s,n] defines the base parameters where s indexes the parameters and n the subjects
    # theta[1,n] <- c
    # theta[2,n] <- r
    # theta[2,n] <- u
    # p[n,k] defines the category probabilities where, n indexes the subject and k the category
 
    p[n,1] <- theta[1,n]*theta[2,n]
    p[n,2] <- (1-theta[1,n])*pow(theta[3,n],2)
    p[n,3] <- (1-theta[1,n])*(2*theta[3,n])*(1-theta[3,n])
    p[n,4] <- theta[1,n]*(1-theta[2,n])+(1-theta[1,n])*pow((1-theta[3,n]),2)
 
    # Designates the data follows a multinomial distribution given the category probabilities response[n,k] is the data matrix
    response[n,1:4] ~ dmulti(p[n,1:4],items)
}
 
# Following loop defines the hierarchical distribution 
 
for (s in 1:3) {
    for (n in 1:subjs) {

# Designates that the base MPT parameters are beta distributed and constraint so that it doesn't crash
    theta[s,n] ~ dbeta(alph[s],bet[s])I(.001,.999) 
    }
}


# Priors for each parameter
# Sets up the priors to be close to flat on the mean and standard deviation of the beta distribution
# The phi is used to keep the variance flat so it samples correctly.

for (s in 1:3) {
 
    zero[s] <- 0
    zero[s] ~ dpois(phi[s])
    phi[s] <- -log(1/pow(alph[s]+bet[s],5/2)) 
    
    alph[s] ~ dunif(.01,5000)
    bet[s] ~ dunif(.01,5000)
    
    mnb[s] <- alph[s]/(alph[s]+bet[s])          # Used to monitor the mean for c,r,u 
    varp[s] <- sqrt(alph[s]*bet[s]/(pow(alph[s]+bet[s],2)*(alph[s]+bet[s]+1)))   # Used to monitor the variance for c,r,u 
    }
}