The Tukey-Lambda distribution doesn’t have a closed form. Is there a way to implement this distribution in Stan?

Hi, we will need a bit more information. Do you have any references we can follow?

1 Like

The first search result:

The quantile function (inverse CDF) has closed form so in theory it’s possible to implement this with Stan’s algebraic equation solver.

Here’s what I came up with. (caveat: this emits a lot a of warnings so a different parameterization may be needed for reliable inference)

```
functions {
real tukey_q(real logit_p, real lambda) {
real log_p = log_inv_logit(logit_p);
real log1m_p = log1m_inv_logit(logit_p);
real s = lambda*log_p;
real t = lambda*log1m_p;
if (abs(s) > 1e-6 || abs(t) > 1e-6)
return (exp(s) - exp(t))/lambda;
else
return (log_p - log1m_p) + (s*log_p - t*log1m_p);
}
real tukey_rng(real lambda) {
real q = uniform_rng(0, 1);
return tukey_q(logit(q), lambda);
}
vector tukey_system(vector logit_p, vector theta, real[] _y, int[] _z) {
return [tukey_q(logit_p[1], theta[1])-theta[2]]';
}
real tukey_lcdf(real x, real lambda) {
vector[1] logit_p = algebra_solver(tukey_system, [0.0]', [lambda, x]', {0.0},{0});
return log_inv_logit(logit_p[1]);
}
real tukey_lpdf(real x, real lambda) {
vector[1] logit_p = algebra_solver(tukey_system, [0.0]', [lambda, x]', {0.0},{0});
real log_p = log_inv_logit(logit_p[1]);
real log1m_p = log1m_inv_logit(logit_p[1]);
return -log_sum_exp((lambda-1)*log_p, (lambda-1)*log1m_p);
}
}
data {
real lambda;
}
parameters {
real x;
}
model {
x ~ tukey(lambda);
}
generated quantities {
real z = tukey_rng(lambda);
}
```

3 Likes