Stan Reference Manual - an easy question on a mixture model

[Stefano]

Hello Stan community,

I approached the deconvolution of a mixed distribution with an example that is present in the reference manual.

The problem is that although the model make sense, it is not able to deconvolve really simple data. There must be something I am missing.

Here I post the R code executable

my.data <- 

  c(
    rnorm(100, 0, 1), 

    rnorm(50,6, 1), 

    rnorm(50, 9, 1)

  )

my.m <- "

    data {

      int<lower=1> K; // number of mixture components

      int<lower=1> N; // number of data points

      real y[N]; // observations

    }

    parameters {

      simplex[K] theta; // mixing proportions

      real mu[K]; // locations of mixture components

    }

    model {

      real ps[K]; // temp for log component densities

      mu ~ normal(0,50);

      for (n in 1:N) {

        for (k in 1:K) ps[k] <- log(theta[k]) + normal_log(y[n],mu[k],1);

        increment_log_prob(log_sum_exp(ps));

      }

    }

  "

my.in <- list(K = 3, N = length(my.data), y = my.data)

fit <- stan(model_code=my.m, data=my.in) 

Thanks a lot!

[Bob Carpenter]

Mixtures can be hard to fit, but that one has enough data
and is well separated enough.

The problem you're likely running into is running with the default four
chains. With mixtures, there's an index-swapping problem. See
the chapter on problematic posteriors.

If you fit with one chain, the fit looks OK:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta[1] 0.27 0.00 0.04 0.20 0.24 0.27 0.30 0.34 656 1
theta[2] 0.50 0.00 0.03 0.43 0.48 0.50 0.52 0.56 1000 1
theta[3] 0.23 0.00 0.04 0.17 0.21 0.23 0.26 0.31 629 1
mu[1] 8.70 0.01 0.17 8.37 8.59 8.70 8.81 9.02 528 1
mu[2] 0.00 0.00 0.10 -0.21 -0.07 0.00 0.06 0.19 526 1
mu[3] 5.90 0.01 0.19 5.51 5.77 5.91 6.04 6.28 616 1
lp__ -477.71 0.08 1.50 -481.62 -478.55 -477.38 -476.57 -475.76 393 1


mu[1] just squeaked into the 97.5% interval, but everything else was in its
50% interval.

- Bob

P.S. You want to define log_theta as a temp and reuse rather than computing
log(theta) N times.

...
model {
real log_ps[K];
vector[K] log_theta;
log_theta <- log(theta);
mu ~ normal(0, 50);

for (n in 1:N) {
for (k in 1:K)

log_ps[k] <- log_theta[k] + normal_log(y[n], mu[k], 1);
increment_log_prob(log_sum_exp(log_ps));

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Andrew Gelman]

An example of a mixture model in Stan is here: http://andrewgelman.com/2015/12/09/why-i-decided-not-to-enter-the-100000-global-warming-time-series-challenge/

[Eric Novik]

That is a really cool example, but there you are passing locations (mu) as data, which makes it a lot easier to estimate the model.

[Bob Carpenter]

The suggestions we've had to identify the model are to order
the theta or order the mu, but it can bias estimates (different
answers than the unordered case) if there's not a big trough
between them.

For instance, changing the declaration of parameter mu from

vector[K] mu;

to

ordered[K] mu;

you get good identification in this artificial data case. Here's
a run with four chains and default parameters:

Inference for Stan model: stefano-mix.
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[1] 0.50 0.00 0.03 0.43 0.48 0.50 0.52 0.57 2417 1
theta[2] 0.23 0.00 0.04 0.17 0.21 0.23 0.26 0.31 2271 1
theta[3] 0.27 0.00 0.04 0.20 0.24 0.27 0.29 0.34 2473 1
mu[1] 0.00 0.00 0.10 -0.20 -0.06 0.00 0.07 0.19 2421 1
mu[2] 5.91 0.00 0.20 5.51 5.78 5.91 6.05 6.30 2091 1
mu[3] 8.71 0.00 0.18 8.37 8.59 8.70 8.83 9.06 4000 1
lp__ -475.03 0.04 1.67 -479.23 -475.92 -474.64 -473.80 -472.92 1382 1

Samples were drawn using NUTS(diag_e) at Fri Jan 8 20:17:48 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).

- Bob

[Stefano]

Thanks Guys,

I wasn't aware of the "index-swapping problem". I think it would be helpful to reference it in the paragraph where the simple mixture deconvolution model is presented, as it is something not intuitive to pick for first Bayes modellers like me.

Bw.

[Bob Carpenter]

OK, will do:

https://github.com/stan-dev/stan/issues/1709#issuecomment-170430146

- Bob