Sampling from a multivariate logistic distribution

I am trying to implement a multivariate logistic function into Stan and
figured I will need the following things in order to make the desired
likelihood function.

v Degrees of freedom
p Dimension of y
T_{p,v} ~ Multivariate T
g_v() Inverse CDF T
\phi^{-1} ~ \Gamma^{-1} (v/2,v/2)
z_{j} = log{\frac{F(\phi^{-1}e_j)}{(1-F(\phi^{-1}e_j))}}
e_1..e_p ~ Multivariate Normal N(0,R)
F(e_j) ~ Uniform distribution (0,1)
L() ~ Logit
R Covariance matrix of e_1..e_p

Some of these are easily accessible by using predefined functions
in Stan. However I’m wondering if it is possible to obtain the
Inverse CDF of a T distribution in an efficient way. Also I do not quite
understand how to sample from a distribution that is undefined. The model
I want to use is build up like this for the p-dimensional case:

Multivariate_logit = T_{p,v} ( g_v(z_1-\mu_1),...,g_v(z_p-\mu_p) ) | 0, R) \times \prod_{j=1}^p [{L(z_j|\mu_j)}/ {T_{1,v}(g_v(z_j-\mu_j)|0,1)}]

So the first part is a multivariate T distribution multiplied by scaling the
univariate logit model by a univariate T distribution. You can see that if
the number of dimensions is p=1, only L(z_1|\mu_1) remains, hence the univariate
logit distribution.

To summarize, does the inverse CDF of the T distribution exist in Stan?
Are there any alternative approaches to solve this. Would it be possible
to define a model in the way defined above and sample from this?
Is it likely that this would lead to efficient sampling?


I don’t think so. The student-t functions are documented here: 15.5 Student-T Distribution | Stan Functions Reference

The steps would be:

  1. Can you write out the log density of the thing you want to sample?
  2. Are the variables you want to sample continuous?
  3. Can you express the constraints you need on those variables in the Stan language?

If so, ideally yes. If not, no.

It’s possible that practical things like inverting a CDF that is numerically finicky make this impossible to do in Stan.

This is a description of the algorithm from which all the Stan algos are derived: 15.1 Hamiltonian Monte Carlo | Stan Reference Manual . Keep in mind it’s all about incrementing a log density.

Here are a couple other sections in the manual that talk about this in terms of the Stan model syntax: 7.3 Increment Log Density | Stan Reference Manual, 7.4 Sampling Statements | Stan Reference Manual

That help at all? You can surround things with dollar signs to use \LaTeX math or double dollar signs to center. You can use tick marks to write preformatted code too.

Has anybody figured out how to do this?

Is the inverse multivariate T cdf needed or just the univariate? The univariate T inverse cdf, what will be called quantile functions with _qf suffix in the next release can be coded in Stan now using

student_t_qf(real u, real nu)  {
  real x = u < 0.5 ? 2 * u : 2 * (1 - u);
  return sign(u - 0.5) * sqrt( nu * (1 / (inv_inc_beta(nu / 2, 0.5, x))) - 1);
1 Like

I think only univariate (looking at this paper -

An easier way to define a multivariate logistic is as a standard location/scale family based on the univariate marginals. See, e.g., Multivariate logistic distribution - Cross Validated and the link therein, Elliptical distribution - Wikipedia

Would it be possible to define multivariate ordered logistic model in Stan? Like in e.g. this paper?

For non-ordinal (binary) response one could use the method described in the manual for multivariate probit here and sample from the PDF of a multivariate logistic rather than multivariate normal.

However, for ordinal responses one would need to generalise the transformed data code - which sorts the data into +'ve and -'ve components whilst keeping track of the indexes - to sort the data into components according to which category each response falls into. Furthermore in the parameters block the latent variable z would need to be broken into components constrained by the thresholds. Would this be possible to code in Stan?

The transformed data component is no problem. I do that kind of thing all the time to avoid branching in the model code.

Im not sure what you mean by latent z being broken into components. Can you be more specific bout what you’re trying to do? I’m afraid I’m not going to have time to digest the paper you linked.

In the multivariate probit code in the Stan manual, z is truncated at 0 so its broken into positive-constrained and negative-constrained components. For ordinal regression it would need to be broken into components defined by cutpoint parameters, so the transformed data and transformed parameters blocks would need to be generalised.