Rstan Stan model for a simple mixture of normals

This blog post from Rbloggers describes how to code a simple three-part normal mixture model with known mixing coefficients, means and standard deviations. While it describes the procedures in some detail for many of the Bayesian sampling tools available in R, the author reported that s/he was unable to sample from this model using rstan. Is it possible to construct this model in rstan? And if so, how?

The model is a mixture of normals with known means and standard deviations: $$p(y)= 0.3\mathcal{N}(-3,2)+0.4\mathcal{N}(2,1)+0.3\mathcal{N}(10,4)$$

Because all of the parameters of this distribution are known, and we merely want to draw samples from this distribution, coding the model in rstan is straightforward. Note that this is, by far, one of the least efficient paths to sampling from this particular model, in terms of time that I spent coding it (15 minutes). The author of the original post is correct when he notes that the easiest way to sample from this model is to use the sample function creatively.

library(rstan)
mix_model   <- "
data{
int                 J;
vector<lower=0>[J]  weights;
vector[J]           means;
vector<lower=0>[J]  sdevs;
}
transformed data{
vector[J]       ln_weights;
ln_weights      <- log(weights);
}
parameters{
real y;
}
model{
vector[J]   probs;
for(j in 1:J){
probs[j]    <- exp(ln_weights[j]+normal_log(y,means[j],sdevs[j]));
}
increment_log_prob(log(sum(probs)));
}
"

mixdata <- list(J=3, weights=c(0.3,0.4,0.3),means=c(-3,2,10),sdevs=c(2,1,4))
testfit <- stan(model_code=mix_model, data=mixdata, iter=10)
fit     <- stan(fit=testfit, data=mixdata, iter=25000, chains=5)


I took the step of reading in each of the parameters of the mixture as data so that the "sum of several known normals" model is easily extended to cases of arbitrary numbers of mixture components.

Transforming the mixture weights to the log scale is done in transformed data because, in this model, it is known. Transforming it there, rather than in the model block means we just read off the stored value, rather than recomputing the log at every iteration.

The only part of this model that I'm unsatisfied with is the loop over the log-probabilities of each mixture component. In general, one prefers to use the native composed functions of rstan because they already have the derivatives worked out, so you don't have to use the slower autodiff routine. On the other hand, the composed function in this case only accepts two arguments, not 3 or more...

y       <- extract(fit, "y")[[1]]
plot(density(y))
x       <- seq(-10,25,by=0.01)
y1      <- 0.3*dnorm(x, mean=-3,sd=2)
y2      <- 0.4*dnorm(x, mean=2,sd=1)
y3      <- 0.3*dnorm(x, mean=10,sd=4)
lines(x,y1, col="red", lty="dashed")
lines(x,y2, col="red", lty="dashed")
lines(x,y3, col="red", lty="dashed")


Visually, the results appear to be a reasonable approximation of the mixture density.