Hi,

I tried putting your implementation to a test via SBC and it appears that both the original version you posted and the one produced with adding `target += log1p(stick_len);`

is somwhat biased, so maybe the Jacobian is off?

Here’s my code:

`ordered_simplex.stan`

```
functions {
vector ordered_simplex_constrain_lp(vector y) {
int Km1 = rows(y);
vector[Km1] x;
real stick_len = 1.0;
// use only negative y's
// start from the smallest -y, which is the last value in the vector
for (k in 1:Km1) {
real adj_y_k = -y[Km1 - k + 1] - log(Km1 - k + 1);
real z_k = inv_logit(adj_y_k);
x[k] = stick_len * z_k;
target += log(stick_len) - log1p_exp(-adj_y_k) - log1p_exp(adj_y_k);
stick_len -= x[k];
}
// this is new,
// instead of adding the stick_len to the last element
// distribute it to all the K - 1 values
// could also just distribute evenly
// comment out the x* = 1 + stick_len and just
// return x + stick_len / Km1;
x *= 1 + stick_len;
return x + (1 - sum(x) ) / Km1;
}
vector ordered_simplex_constrain_v2_lp(vector y) {
int Km1 = rows(y);
vector[Km1] x;
real stick_len = 1.0;
// use only negative y's
// start from the smallest -y, which is the last value in the vector
for (k in 1:Km1) {
real adj_y_k = -y[Km1 - k + 1] - log(Km1 - k + 1);
real z_k = inv_logit(adj_y_k);
x[k] = stick_len * z_k;
target += log(stick_len) - log1p_exp(-adj_y_k) - log1p_exp(adj_y_k);
stick_len -= x[k];
}
// this is new,
// instead of adding the stick_len to the last element
// distribute it to all the K - 1 values
// could also just distribute evenly
// comment out the x* = 1 + stick_len and just
// return x + stick_len / Km1;
x *= 1 + stick_len;
target += log1p(stick_len);
return x + (1 - sum(x) ) / Km1;
}
}
data {
int K;
int<lower=1,upper=2> version;
int<lower=0> observed[K];
}
parameters {
positive_ordered[K] y;
}
transformed parameters {
simplex[K] x;
if(version == 1 ) x = ordered_simplex_constrain_lp(y);
else x = ordered_simplex_constrain_v2_lp(y);
}
model {
x ~ dirichlet(rep_vector(2, K));
observed ~ multinomial(x);
}
```

and the SBC code:

```
library(SBC) # remotes::install_github("hyunjimoon/SBC")
library(cmdstanr)
library(MCMCpack)
library(ggplot2)
library(future)
plan(multisession)
options(SBC.min_chunk_size = 5)
m <- cmdstan_model("ordered_simplex.stan")
backend <- SBC_backend_cmdstan_sample(m, chains = 2)
generate_one_dataset <- function(N, K, version) {
x_raw <- rdirichlet(1, alpha = rep(2, K))
x <- sort(x_raw)
observed <- as.integer(rmultinom(1, size = N, prob = x))
list(
parameters = list(x = x),
generated = list(K = K, observed = observed, version = version)
)
}
datasets_v1 <- generate_datasets(
SBC_generator_function(generate_one_dataset, N = 30, K = 6, version = 1),
n_datasets = 500)
res_v1 <- compute_results(datasets_v1, backend)
plot_ecdf_diff(res_v1) + ggtitle("V1")
plot_rank_hist(res_v1) + ggtitle("V1")
plot_sim_estimated(res_v1, alpha = 0.3) + ggtitle("V1")
datasets_v2 <- generate_datasets(
SBC_generator_function(generate_one_dataset, N = 30, K = 6, version = 2),
n_datasets = 500)
res_v2 <- compute_results(datasets_v2, backend)
plot_ecdf_diff(res_v2) + ggtitle("V2")
plot_rank_hist(res_v2) + ggtitle("V2")
plot_sim_estimated(res_v2) + ggtitle("V2")
```

The ranks for “v1”:

The ranks for “v2”:

So either I misunderstood how to use the code, or how to simulate an ordered simplex or there is some problem with my code or there is a problem with your code :-)

Hope that’s at least slightly helpful :-)