Tukey-Lambda Distribution

Modeling
1

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

2

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

1 Like
3

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