### Generate artificial data for a HierMixModel.
library(rje)

hmmgen <- function (I,alpha0,mu0,beta0,tau0,gamma0,alphaN,
                    lambda=50) {
  K <- length(alpha0)
  if (length(mu0)!=K || length(beta0)!=K ||
        length(tau0)!=K || length(gamma0)!=K ||
        length(alphaN)!=1) {
    stop("Not all parameters have length ",K)
  }
  ## Figure out size of sample
  J <- rpois(I,lambda)
  Jmax <- max(J)
  Y <- matrix(NA,I,Jmax)
  
  ## draw level 1 parameters
  pi <- rdirichlet(I,alphaN*alpha0)
  mu <- matrix(NA,I,K)
  tau <- matrix(NA,I,K)
  for (k in 1:K) {
    mu[,k] <- rnorm(I,mu0[k],beta0[k])
    tau[,k] <- exp(rnorm(I,tau0[k],gamma0[k]))
  }
  
  ## Draw data
  for (i in 1:I) {
    Z <- sample(1:K,J[i],replace=TRUE,prob=pi[i,])
    Y[i,1:J[i]] <- rnorm(J[i],mu[i,Z],1/sqrt(tau[i,Z]))
  }
  
  list(Y=Y,J=J,I=I,pi=pi,mu=mu,tau=tau,
       alpha0=alpha0,mu0=mu0,beta0=beta0,
       tau0=tau0,gamma0=gamma0,alphaN=alphaN)
}