Ordered simplex constraint transform

Modeling
18

Hi,
just returning to this topic as I was looking at this for another reason. I noticed that the way to make the softmax approach suggested by @Bob_Carpenter work is to note that I actually want to compute the Jacobian for transforming the positive_ordered vector to the first/last K-1 elements of the simplex - the remaining element of the simplex is then simply 1 - sum(rest) which is a transformation with constant Jacobian which can be dropped.

The log determinant of the K-1 \times K - 1 Jacobian matrix then has an analytical form (I used Mathematica to help me with it). The full code is then:

data {
  int K;
  array[K] int<lower=0> observed;
  vector<lower=0>[K] prior_alpha;
}


parameters {
  positive_ordered[K - 1] y;
}

transformed parameters {
  simplex[K] x = softmax(append_row(0, y));
}

model {
  x ~ dirichlet(prior_alpha);
  // The required log Jacobian determinant
  target += sum(y) - K * log_sum_exp(append_row(0, y));
  // Likelihood, because why not
  observed ~ multinomial(x);
}

I verified that it works with SBC (code almost identical as the above cases)
In quick tests it appears to be slightly slower, less stable and inducing worse geometry than the variant I proposed in a previous post (Ordered simplex constraint transform - #14 by martinmodrak) but I thought someone may find this interesting.

The problem in the code tested by @spinkney is that the full K \times K Jacobian is not full rank and has determinant 0, which explains the numerical issues.

1 Like