4
$\begingroup$

I'm trying to model data $0 < Y_i < 1$ with a finite mixture of Beta components. To do this, I've adapted the code given in section 5.3 of the Stan manual. Instead of (log)normal priors, I am using $\mathrm{Exponential}(1)$ priors for the $\alpha$ and $\beta$ parameters. Thus, as I understand it, my model is as follows:

\begin{align*} \alpha_k, \beta_k &\overset{iid}{\sim} \mathrm{Exponential}(1) \\ Z_i &\sim \mathrm{Categorical}(1, \ldots, K) \\ Y_i \mid \left(Z_i = k\right) &\sim \mathrm{Beta}_{\alpha_k, \beta_k} \end{align*}


Now, for my implementation in stan, I have the following two code chunks:

# fit.R
y <- c(rbeta(100, 1, 5), rbeta(100, 2, 2))
stan(file = "mixture-beta.stan", data = list(y = y, K = 2, N = 200))

and

// mixture-beta.stan

data {
  int<lower=1> K;
  int<lower=1> N;
  real y[N];
}

parameters {
  simplex[K] theta;
  vector<lower=0>[K] alpha;
  vector<lower=0>[K] beta;
}

model {
  vector[K] log_theta = log(theta);

  // priors
  alpha ~ exponential(1);
  beta ~ exponential(1);
  
  for (n in 1:N) {
    vector[K] lps = log_theta;

    for (k in 1:K) {
      lps[k] += beta_lpdf(y[n] | alpha[k], beta[k]);
    }

    target += log_sum_exp(lps);
  }
}


After running the code above (defaults to 4 chains of 2000 iterations, with 1000 warmup) I find that all the posterior components are essentially the same:

> print(fit)
Inference for Stan model: mixture-beta.
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.01 0.13  0.26  0.42  0.50  0.58  0.75   259 1.01
theta[2]  0.50    0.01 0.13  0.25  0.42  0.50  0.58  0.74   259 1.01
alpha[1]  2.40    0.38 1.73  0.70  0.94  1.20  3.89  6.01    21 1.16
alpha[2]  2.57    0.37 1.74  0.70  0.96  2.29  4.01  6.05    22 1.16
beta[1]   3.54    0.11 1.10  1.84  2.66  3.46  4.26  5.81    93 1.04
beta[2]   3.58    0.12 1.07  1.88  2.77  3.49  4.26  5.89    82 1.05
lp__     30.80    0.05 1.74 26.47 29.92 31.21 32.08 33.02  1068 1.00

Samples were drawn using NUTS(diag_e) at Thu Sep 17 12:16:13 2020.
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).

I read the warning about label switching, but I can't see how to use the trick of ordered[K] alpha since I also need to integrate the constraint of $\alpha$ and $\beta$ being positive.

Could someone help explain what's going on here?

$\endgroup$
1
$\begingroup$

I haven't (and won't) checked what I'm saying in Stan (I haven't time to spend on compilation today!), so please try this out and let us know what happens.

First, I'm pretty sure you're right that the issue is label switching. You should plot the traceplots (traceplot(my_stan_fit)) to confirm this. Basically, on some chains, alpha[1] and beta[1] belong to the high-probability distribution, while in others they belong to the low probability distribution.

Second, I think you can set constaints on ordered vectors, e.g. ordered<lower=0>[K] alpha;.

Third, rather than enforcing alpha[1] < alpha[2] and beta[1] > beta[2], it's probably more effective to create a transformed parameter encoding the mean of each of your mixture distributions and enforce and order on this, e.g. something like (again, I haven't tried to compile this):

transformed parameters { 
    ordered<lower=0,upper=1> mu[K];
    for (k in 1:K) {
        mu[k] = alpha[k] / (alpha[k] + beta[k]);
    }
} 
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ +1. One of the simplest possible expedients would be to let the $\alpha_k$ be the cumulative sum of the exponential variates, thereby assuring they are always ordered from smallest to largest. That would require a change at only one place in the code, replacing alpha[k] by cumsum(alpha)[k] (although a more efficient implementation would compute the cumulative sum just once upon generating alpha). $\endgroup$ – whuber Sep 28 at 14:06
  • $\begingroup$ Oh, interesting! Would that complicate the priors though? $\endgroup$ – Eoin Sep 28 at 14:09
  • $\begingroup$ My understanding is that these are rather generic priors, intended to function "uninformatively," so it shouldn't make any difference provided they have sufficiently large SDs. Indeed, it would be nice to treat the parameters in a symmetrical fashion--but note that the proper meaning of symmetry here may be complex, because one would expect the components with low mean to have high positive skew and those with high mean to have high negative skew. $\endgroup$ – whuber Sep 28 at 14:14
  • $\begingroup$ Another option you could consider is post-processing the MCMC output to find consistent labels across the iterations. If you are using R, the label.switching package provides a couple of tools for this. $\endgroup$ – baruuum Sep 29 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.