Better way of modelling stick-breaking process?

Modeling
2

Why do you want to model a stick-breaking process? It’s built into our simplex data type when you do this:

simplex[K] theta;

That makes sure that theta[k] > 0 for all k in 1:K and makes sure sum(theta) = 1. It works under the hood as a stick-breaking process as outlined in the manual.

Sorry nobody answered this before, but you’re absolutely right. You do not want to be using this kind of hack.

As it should be. You have a non-linear transform on the left-hand side of a sampling statement, which is exactly what this warning is supposed to be doing!

You can’t assign to data. So I’m not even sure what you’re trying to do here.

You can vectorize your b sampling to just b ~ beta(1, a); and same for normal, but you don’t want to do that normal trick.

What you want is something like this to draw the breaks from a beta distribution:

parameters {
  vector<lower=0, upper=1>[N - 1] breaks;
  ...
transformed parameters {
  simplex[N] b;
  {
    real sum = 0;
    b[1] = breaks[1];
    sum = b[1];
    for (n in 2:(N - 1)) {
      b[n] = (1 - sum) * breaks[n];
      sum += b[n];
    }
    b[N] = 1 - sum;
  }
  ...
model {
  breaks ~ beta(1, a);

I didn’t test that, but it’s the right way to lay the code out for stick-breaking, with the result being written into the simplex b. You can just loop if you need to break more sticks.

1 Like