Loglogistic distribution

[olivier....@student.uclouvain.be]

Hello,

I read in a previous post that one can always do recommendations for adding new distributions. I would need the loglogistic distribution in order to fit a growth curve model, but it is not yet available in Stan.

Kind regards,

Olivier

[Bob Carpenter]

If it's this one:

https://en.wikipedia.org/wiki/Log-logistic_distribution

then you can define its log density yourself with a user-defined
function.


functions {

real loglogistic_log(real y, real a, real b) {
return log(b) - log(a) + (b - 1) * (log(y) - log(a))
- 2 * log1p_exp(b * (log(y) - log(a));
}

}

For more efficiency, you want to store and reuse repeated
subexpressions like log(b) and (log(y) - log(a)).
For example, I think that top term reduces to:

log(b) - log(y) + b * (log(y) - log(a))

and then you can share the product term with the log1p_exp.
Oh, and log1p_exp(u) is log(1 + exp(u)), which allows us to
work on the log scale with u.

Once you know it works, you can put the functions block
in front of a Stan program and use

y ~ loglogistic(a, b);

in the model block. It won't be vectorized.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[olivier....@student.uclouvain.be]

Thank you for your response! Normally, I would only need the CDF of the loglogistic distribution. This is for instance the model for a weibull growth curve:

mu[i] <- ult[origin[i]] * weibull_cdf(dev[i], omega, theta);

Now, I tried to add the following to the beginning of my code :

model <- "

functions { 

  

real loglogistic_cdf(real y, real b, real a) { 

return 1/(1+(y/a)^(-b)); 

} 

} 

data { ...

...

mu[i] <- ult[origin[i]] * loglogistic_cdf(dev[i], omega, theta);

...

"

Do you think this is an efficient implementation? Sorry for the potentially stupid questions, I am just getting started with Stan.

[Bob Carpenter]

You almost always want to work on the log scale to prevent
underflow or overflow, and be careful about precision, so that'd
be:

loglogistic_cdf_log(real y, real a, real b)
return -log1p_exp(-b * (log(y) - log(a));
}

and then you can use it with truncation notation, but you won't
need that if you only need the CDF. Then you can calculate
the log of mu:

log_mu[i] <- log(ult[origin[i]]) + loglogistic_cdf_log(dev[i], omega, theta);

You don't then indicate what you want to do with mu, but you
want to avoid going back to the linear scale exp(log_mu[i]) if you
can at all help it.

- Bob

[olivier....@student.uclouvain.be]

Thank you very much! I hadn't heard of these concepts before!

[Bob Carpenter]

log_sum_exp is critical for computational stats when
you want to avoid underflow to 0 or overflow to 1 in
probability calculations involving sums on the log scale.

log_sum_exp(a, b) = log(exp(a) + exp(b))

= max(a,b) + log(exp(a - max(a,b)), exp(b - max(a, b))

You can see that the exp() now can't overflow since it gets an
argument less than or equal to 0 and you can see that the max(a,b)
term is going to get the leading digits of the result.


So to work with addition and multiplication on the log scale,
we have this

log(a * b) = log(a) + log(b)

log(a + b) = log(exp(log(a)) + exp(log(b))

= log_sum_exp(log(a), log(b))

The first expression computes the log of a*b given the log of
a and the log of b, and the second expression does this for addition.

- Bob