To test some ideas I threw around in another thread, I’m trying to write a special-purpose sampler using some of those ideas and test it against what Stan gives me. But, my problem of interest has a lot of simplexes. I looked at how Stan unconstrains simplexes but was thinking to do something different.

As I understand it, if y_i ~ gamma(1,1) then s_i = y_i/sum(y_i) is dirichlet(1_vector), which is uniform on the simplex, so if p(s) is a density on the simplex then I get what I want by p(y_i/sum(y_i))

so if x is unconstrained on -inf,inf and I constrain to [0,inf] with y = exp(x), then

p(x) dx = gamma_pdf(exp(x),1,1)*dy/dx dx

gives a gamma distribution on y = exp(x)

dropping normalization constants, this is:

exp(x)^0 * exp(-exp(x)/1) exp(x) = exp(-exp(x) + x)

so for x in [-inf,inf] I can use lp of -exp(x) + x

Does anyone see anything wrong with this calculation? Is there a reason Stan doesn’t use this reparameterization? overflow as x goes positive?