(This article was originally published at Statistical Modeling, Causal Inference, and Social Science, and syndicated at StatsBlogs.)
There was some confusion on the Stan list that I wanted to clear up, having to do with fitting mixture models. Someone quoted this from John Kruschke’s book, Doing Bayesian Data Analysis:
The lack of discrete parameters in Stan means that we cannot do model comparison as a hierarchical model with an indexical parameter at the top level. There might be ways to work around this restriction by using clever programming contrivances, but presently there is nothing as straight forward as the model specification in JAGS.
Actually it’s easy to write hierarchical mixture models in Stan. Here’s the code (supplied by Bob Carpenter) for the specific model from Kruschke:
data {
int N;
int y[N];
}
parameters {
real<lower=0, upper=1> theta;
}
model {
y ~ bernoulli(theta);
target += log_mix(0.5, beta_lpdf(theta | 2.75, 8.25),
beta_lpdf(theta | 8.25, 2.75));
}
Here’s the fit in R:
N = 9
z = 6
y = c(rep(0, N - z), rep(1, z))
fit = stan("kruschke-db.stan")
print(fit)
Inference for Stan model: kruschke-db.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta 0.67 0.00 0.14 0.32 0.6 0.70 0.77 0.87 1409 1
lp__ -7.69 0.02 0.77 -9.69 -8.0 -7.38 -7.14 -7.09 2089 1
By comparison, here’s Krushke’s model in Jags:
model {
for (i in 1:N) {
y[i] ~ dbern(theta)
}
theta ~ dbeta(omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1)
omega[1] <- .25
omega[2] <- .75
kappa <- 12
m ~ dcat(mPriorProb[])
mPriorProb[1] <- .5
mPriorProb[2] <- .5
}
The Stan discussion thread is here.
Speaking more generally on mixture models, Bob adds:
There are huge advantages in efficiency from marginalization, and in some cases it allows us to do inference in models where discrete sampling gets stuck. In other cases, discrete sampling is challenging if not impossible, as in Ising models or the combinatorially similar variable selection models.
It's come up with mixtures, the Cormack-Jolly-Seber model in ecology, and change-point models, all of which are explained in the manual. All of Kery and Straub's book has been translated and also Wagenmakers and Lee, both of which contain quite a few BUGS/JAGS models with discrete parameters.
Bob also points out that what's "very easy" for one person is a "contrivance" for another. And this goes in both directions. If you're already familiar with some software, whether it's Excel or Emacs or Mathematica or whatever, and it does the job for you, then you might as well continue using it. But if the software doesn't do everything you need, then it's time to learn something new.
The point of this post is that if Stan is working for you, or it could be working for you, but you heard that "we cannot do model comparison as a hierarchical model with an indexical parameter at the top level" then, no, this is not actually something to worry about. You actually can fit these models in Stan, it's completely straightforward. If you don't want to use Stan, that's fine too. We wrote Stan because existing programs couldn't fit the models we want to fit. Also we are planning to write a mixture model function in Stan that does exactly what's shown in the example above but doesn't have the "lpdf" things in case that upsets people.
(28 October: Stan model code fixed by Aki to show the constraint for theta)
The post Yes, despite what you may have heard, you can easily fit hierarchical mixture models in Stan appeared first on Statistical Modeling, Causal Inference, and Social Science.
Please comment on the article here: Statistical Modeling, Causal Inference, and Social Science
