How do I specify a logit-normal distribution

Modeling
1

I want to fit a logit-normal distribution to a data set.

logit(y/ymax) ~ normal(mu,sigma)

where mu, sigma and ymax are the parameters. (This distribution is also called an upper limit log normal distribution.)

I have implemented it using a Jacobian adjustment below. Michael Betancourt adviced me against this approach on twitter, and suggested:

Can anybody help me interpret this advice?

data {
  int<lower=0> N; // number of events
  real<lower=0> y[N]; // observed value
  real<lower=0> yobsmax; //highest observed value 
}
parameters {
  real<lower=0> sigma;
  real mu;
  real<lower=yobsmax> ymax;
}
model {
    target += normal_lpdf(mu | 0, 5);
    target += lognormal_lpdf(sigma | 0, 5);
    target += -ymax;
    for (j in 1:N) {
        target += lognormal_lpdf( (y[j]*ymax)/(ymax-y[j]) | mu, sigma);
        target += log(y[j]) + log(ymax) - log(ymax - y[j]);
    }
}
  • Surrogate Data *
    this data was generated with Ī¼=0, Ļƒ=3, ymax=300.
y <- c(0.126200, 0.118300, 7.778100, 0.036400, 0.096600, 83.073800, 7.079400, 0.294000, 1.949300, 4.726800, 3.918100, 0.005800, 5.822200, 52.167300, 0.747500, 0.997900, 56.475600, 2.597600, 0.003100, 2.578700, 46.248000, 4.251200, 0.107100, 204.197900, 29.969100, 0.118900, 3.834800, 2.496200, 3.833400, 11.984300, 36.084200, 39.096500, 0.371500, 25.997400, 0.002300, 4.884200, 0.546900, 0.795900, 2.336200, 1.005800, 9.580700, 19.440900, 0.180800, 14.342300, 0.523200, 137.215500, 55.234000, 0.666800, 10.892300, 0.088300, 3.370900, 20.666900, 0.845300, 9.211700, 5.667300, 0.474600, 3.149000, 0.463000, 0.119900, 1.643800, 0.066500, 2.036400, 2.272600, 8.318800, 61.534800, 0.005200, 0.013600, 2.572100, 0.761400, 0.029100, 0.060200, 0.121700, 0.216400, 6.413700, 0.771200, 0.056700, 0.066400, 15.507000, 0.497900, 0.031600, 9.732200, 0.442500, 2.119500, 3.947000, 0.676800, 2.086400, 5.764600, 3.004400, 13.824000, 0.574900, 3.246000, 1.723700, 0.242900, 1.225500, 0.032600, 0.007900, 0.831700, 2.563500, 0.017600, 0.197600, 0.086300, 0.020000, 19.560600, 0.025700, 0.551300, 0.046000, 0.429000, 0.010600, 0.032000, 0.005900, 0.107900, 6.545200, 0.277000, 136.810000, 3.865500, 0.272500, 0.077200, 0.010700, 74.743900, 14.806200, 3.562500, 0.086000, 0.033400, 0.004500, 0.376200, 7.004900, 1.629100, 1.542000, 28.362300, 0.809300, 0.092400, 0.015200, 0.825600, 0.041800, 7.593200, 9.530400, 6.149100, 0.274000, 0.265900, 0.030100, 3.061100, 13.143200, 7.629600, 36.954100, 0.815500, 0.507200, 0.385300, 0.591200, 0.734600, 0.064500, 280.928900, 0.009400, 1.182500, 0.648600, 4.428400, 0.738300, 7.708000, 0.004200, 7.156400, 0.691000, 9.833700, 0.123900, 0.078900, 3.676300, 0.863700, 0.009500, 1.277500, 14.559000, 0.010400, 76.919400, 0.652900, 0.796800, 0.049800, 4.488900, 7.441300, 1.565300, 0.016300, 60.066200, 0.068900, 0.068000, 0.949300, 7.018900, 16.992700, 1.052800, 11.426600, 0.377700, 29.089100, 0.126200, 0.022500, 0.006400, 1.825400, 0.000300, 93.251800, 155.191100, 0.009300, 10.605800, 1.708700, 44.859400, 71.884900, 15.400800, 0.663200, 0.001400, 17.293200, 3.056700, 29.001300, 0.132400, 98.233800, 0.718400, 2.989800, 17.986700, 0.005700, 0.006400, 5.533600, 16.106500, 3.164500, 0.741200, 1.551600, 62.235400, 0.719400, 28.484000, 4.372600, 0.000700, 0.035000, 0.013400, 218.007500, 127.601500, 0.002000, 141.031200, 0.722200, 16.013700, 3.598500, 23.820800, 5.317500, 1.556500, 0.163800, 0.573800, 0.008700, 1.620100, 3.308600, 1.329900, 0.210400, 1.416500, 20.737000, 0.251200, 3.843800, 26.738100, 0.120200, 2.535800, 0.768000, 3.150000, 73.607800, 0.680400, 37.693600, 9.588500, 0.349500, 8.256300, 0.011200, 53.254400, 39.776500, 1.741900, 2.736000, 23.374000, 0.079000, 0.133200, 0.481400, 7.267300, 0.023700, 49.941700, 0.084100, 0.079400, 8.241300, 11.766100, 91.412500, 13.292900, 1.411700, 0.210600, 16.738800, 0.034400, 6.167000, 0.148700, 13.635400, 6.065600, 0.083400, 0.263900, 0.022700, 2.386100, 0.138900, 6.630300, 2.640200, 0.053800, 0.215000, 0.107200, 0.241400, 1.471300, 2.988400, 3.448000, 14.083800, 0.031400, 27.396700, 0.139500, 0.041000, 0.140300, 0.001200, 85.915600, 95.288000, 6.191700, 0.061700, 2.237000, 3.462400, 10.992700, 0.063700, 0.000100, 3.941600, 8.039300, 0.001100, 3.785300, 0.088800, 10.573700, 39.616000, 2.797700, 0.463300, 0.196900, 0.074800, 6.791800, 2.524800, 260.086800, 0.236000, 2.392200, 0.001600, 0.252700, 0.272200, 0.089300, 1.071400, 47.590100, 23.691800, 0.779000, 0.683200, 8.491100, 0.771200, 0.084800, 0.134800, 0.851400, 0.034400, 0.575800, 0.042200, 0.357800, 0.233600, 0.203300, 0.035700, 0.131300, 66.652300, 0.024600, 1.092700, 0.049900, 298.905200, 0.328000, 0.430300, 5.383100, 0.081200, 0.069000, 61.124700, 0.569700, 0.007500, 0.151700, 0.309200, 99.214900, 0.908400, 0.046900, 0.043200, 0.209500, 5.141900, 76.921600, 0.177700, 0.211800, 0.262600, 10.779900, 2.080100, 4.621700, 0.009100, 55.759000, 0.029000, 0.127400, 0.101500, 0.170100, 5.244100, 7.910100, 5.758200, 0.014500, 0.056700, 1.571000, 0.376900, 6.089500, 0.399300, 0.286500, 0.075100, 8.230300, 4.589900, 31.676700, 0.013000, 0.605100, 0.297400, 0.035000, 3.436300, 5.233400, 0.003400, 0.070200, 35.251000, 0.541000, 1.357800, 42.871000, 0.028300, 0.520600, 0.037000, 0.209400, 8.395700, 0.045000, 0.353800, 1.281200, 0.685500, 3.175600, 0.448200, 0.049000, 0.447700, 49.674900, 0.362800, 0.886700, 2.431500, 0.356000, 1.848200, 0.001700, 3.884600, 20.572400, 0.100700, 1.151200, 0.100800, 0.029000, 0.053900, 0.047500, 6.431100, 99.930800, 0.145400, 1.748800, 0.058600, 20.170900, 199.048000, 3.980400, 0.842900, 32.573800, 1.969300, 0.072800, 0.551700, 1.394600, 8.853300, 0.031100, 0.985600, 5.201900, 0.103700, 0.689900, 0.813800, 0.455000, 2.666900, 0.031400, 0.178700, 0.413500, 229.050500, 78.518500, 0.095000, 101.611800, 7.094600, 0.648100, 3.830900, 10.764500, 0.821500, 9.856800, 0.143400, 0.015300, 0.062400, 99.875100, 3.095300, 0.637900, 5.025700, 0.470600, 0.057500, 0.540200, 0.132000, 5.708600, 0.009300, 4.499800, 9.615600, 0.345500, 0.021100, 14.496500, 0.074400, 0.119700, 8.814100, 3.169000, 0.003800, 43.916300, 0.133900, 4.038500)
N <- 500
yobsmax <- 299
2

Thereā€™s a discussion in the userā€™s guide part of the manual in the ā€œcorrelated topic modelā€ example.

The basic idea is that to generate theta from a logistic normal, use a latent unconstrained parameter X and then run through softmax.

X ~ normal(mu, Sigma);
theta = softmax(X);
2 Likes
3

I encountered a very similar problem as yours, so I wondered if you have solved the problem. If so could you please share the solution? It will be very much appreciated! :)

1 Like
4

I mistakenly elaborated with the multivariate version. For the univariate case where you want

\mathrm{logit}(\alpha) \sim \mathsf{Normal}(\mu, \sigma),

it looks like this:

parameters {
  real<lower = 0, upper = 1> alpha;
...
model {
  logit(alpha) ~ normal(mu, sigma);
  target += log_inv_logit(alpha) + log_1m_inv_logit(alpha);  // log Jacobian

The better way to code this is to make logit_alpha the parameter and define alpha as a transform; that way you donā€™t need the Jacobian adjustment because thereā€™s no change of variables in a distribution.

parameters {
  real logit_alpha;
  ...
transformed parameters {
  real alpha = inv_logit(logit_alpha);
  ...
model {
  logit_alpha ~ normal(mu, sigma);

@Aslak_Grinstedā€™s case is confusing, because itā€™s essentially putting a lot of distributions on ymaxā€”one for each N. I donā€™t see how you could define a parameter logit_y_div_ymax and hope to recover ymax from it.

parameters {
  vector[N] logit_y_div_ymax;
  ...
transformed parameters {
  vector[N] ymax = ... some function of y and logit_y_div_ymax ... ???

P.S. Also, Iā€™d suggest not taking any chances and defining yobsmax in transformed data as:

transformed data {
  real<lower = 0> yobsmax = max(y);

or even just defining the bound directly in the declaration:

real<lower = max(y)> ymax;

P.P.S. You can also vectorize the loop:

(y .* ymax) ./ (ymax - y) ~ lognormal(mu, sigma);  // drops constant terms, but is more efficient
target += log(y) + N * log(ymax) - log(ymax - y);

The .* and ./ are the elementwise vector multiplication and division. And of course, given that ymax is a parameter, you need a Jacobian here (though it may simplify given that y is data).

1 Like
5

I agree with Bob when ymax = 1, which is the same as logit(alpha). When ymax is not 1, but instead a parameter value (which I believe is @Aslak_Grinstedā€™s case), one can use the following function I wrote with specification of the mean in the model statement.

    function{
        real scaled_logit_normal_lpdf(real y, real mu, real sigma, real ymax){
            real logit_ymax = logit(y/ymax);
            return normal_lpdf(logit_ymax | mu, sigma) + log(ymax) - log(y) - log(ymax-y);
          }
        }

In the model block, specify mu according to oneā€™s model and pass it to the function above.

for(i in 1:n1){
    mu[i] = beta1[1] + beta1[2]*logd1[i];
    y[i] ~ scaled_logit_normal_lpdf(mu[i], std, ymax);
  }

This way, ymax can be estimated nicely.

1 Like
6

@Bob_Carpenter thereā€™s a typo above. I think it should be log1m_inv_logit.