RNG function for induced dirichlet distribution

For my R package I had to implement an RNG function for @betanalpha induced dirichlet distribution for ordinal cutpoints–see link. The RNG function is in the attached Stan code along with the PDF to show parameter recovery.

Not sure to what extent this will be more widely useful but wanted to post it here in case it is. Thanks again to @betanalpha for creating this useful distribution.
induced_dirichlet_rng.stan (2.0 KB)

2 Likes

There is a rng function in that case study, is it not the same?

Reproduced below

functions {
  vector induced_dirichlet_rng(int K, real phi) {
    vector[K - 1] c;
    vector[K] p = dirichlet_rng(rep_vector(1, K));
    
    c[1] = phi - logit(1 - p[1]);
    for (k in 2:(K - 1))
      c[k] = phi - logit(inv_logit(phi - c[k - 1]) - p[k]);
      
    return c;
  }
}

transformed data {
  int<lower=1> N = 50; // Number of observations
  int<lower=1> K = 5;  // Number of ordinal categories
}

generated quantities {
  real gamma = normal_rng(0, 1);                  // Latent effect
  ordered[K - 1] c = induced_dirichlet_rng(K, 0); // (Internal) cut points
  int<lower=1, upper=K> y[N];                     // Simulated ordinals
  
  for (n in 1:N)
    y[n] = ordered_logistic_rng(gamma, c);
}
2 Likes

Haha that’s awesome. I literally just reinvented the wheel.

On the other hand, at least I did it correctly!

1 Like