## -*- mode: R; -*-
############################################################
## General K component mixture model, no covariates.
############################################################
## p-variant -- strong priors on level 2 parameters
##
## Fixed Structural Parameters
## I -- number of essays, i -- index for essays
## J[I] -- number of events per essay, j -- index for events,
## Jmax <- max(J)
## K -- number of mixture components
##
## Data
## Y[I,Jmax] -- log pause time
## Note that this is a ragged array, so that Y[i,j] for j>J[i] will be NA
##
## 1st Level parameters
## pi[I,K] -- mixture probabilities
## mu[I,K] -- component means
## tau[I,K] -- component precisions
## sigma[I,K] <- 1/sqrt(tau[I,K])
##
## 2nd Level hyperparameters
## alpha[K] -- mixture component proportions
## mu0[K] -- grand means for components
## beta0[K] -- precisions for components
## tau0[K] -- grand geometric mean for precisions
## gamma0[K] -- geometric SD for precisions
##
## 3rd level fixed hyperparameters
## Need to convert sds into precisions for JAGS/Stan compatability
data {
   for (k in 1:K) {
     mu0p[k] <- 1/(mu0s[k]*mu0s[k])     # mu0 prior precision
     lbeta0p[k] <- 1/(lbeta0s[k]*lbeta0s[k]) # log beta0 prior precision
     tau0p[k] <- 1/(tau0s[k]*tau0s[k])   # log tau0 prior precision
     lgamma0p[k] <- 1/(lgamma0s[k]*lgamma0s[k]) # log gamma0 prior precision
   }
}

model{

  ## Hyperpriors for mixtures parameters
  alpha0[1:K]~ddirch(alpha0mu[1:K])
  alphaN ~dchisqr(alphaNu)T(0,500)
  for (k in 1:K) {
    mu0[k]~dnorm(mu0m[k],mu0p[k])
    beta0[k]~dlnorm(lbeta0m[k],lbeta0p[k])
    tau0[k] ~ dnorm(tau0m[k],tau0p[k])
    gamma0[k] ~ dlnorm(lgamma0m[k],lgamma0p[k])
    alpha[k] <- alpha0[k]*alphaN
  }

  ## Priors for individual components
  for(i in 1:I){
    pi[i,1:K] ~ ddirch(alpha[]+.01)
    for (k in 1:K) {
      theta[i,k] ~dnorm(0,1)
      mu[i,k] <- mu0[k]+theta[i,k]*beta0[k]
      eta[i,k] ~dnorm(0,1)
      tau[i,k] <- exp(tau0[k]+eta[i,k]*gamma0[k])
    }
  }

  ## Likelihood for data
  for (i in 1:I) {
    for (j in 1:J[i]) {
###     Z[i,j] ~ dcat(pi[i,1:K])
###     Y[i,j] ~dnorm(mu[i,Z[i,j]], tau[i,Z[i,j]])
      Y[i,j] ~ dnormmix(mu[i,],tau[i,],pi[i,])
    }
  }

}