What is the best package for fitting a model where both the mean and the variance are latent variables? I'm trying to model skew that arises from correlation, i.e., a positive correlation between mean and variance causes positive skew. For example, m ~ a_m + b_m . x + e_m s ~ a_s + b_s . x + e_s y ~ N(m, s) cov(e_m, e_s) > 0, and m and s are not directly observable. [[alternative HTML version deleted]]