Possible to get a multivariate gamma density in STAN?

[Royce Anders]

Hello, 

Would it be possible to get a multivariate gamma density in STAN?  

I am aware that the Wishart is a special case of this distribution, though I'm interested to use the more general case (mvgamma) where the shape and rate are specifiable. 

Best wishes,

Royce

[Ben Goodrich]

On Wednesday, November 4, 2015 at 10:02:38 AM UTC-5, Royce Anders wrote:

Would it be possible to get a multivariate gamma density in STAN?  

Coding-wise, it would not be difficult at all to copy what we have for the Wishart distribution and modifying it slightly for the multivariate gamma. The problem is finding someone who is sufficiently motivated to do it because in order to put it into Stan it would also have to be tested at least as much as the Wishart distribution, documented, complimented with a function to draw from the multivariate gamma, etc. And that is all pretty tedious to do for a distribution that is not widely used and is not likely to lead to efficient sampling in Stan unless X were data.

If you are that motivated, then we can help you jump through the hoops involved when adding new functions to Stan. But if you are not so motivated, it is pretty easy to write that log density as a user-defined Stan function within a Stan program. And depending on what special case (diagonal scale matrix, etc.) you might be interested in, it may actually run faster than would a general implementation in the Stan Math Library.

Ben

[Bob Carpenter]

> On Nov 4, 2015, at 11:11 AM, Ben Goodrich <goodri...@gmail.com> wrote:
>
> On Wednesday, November 4, 2015 at 10:02:38 AM UTC-5, Royce Anders wrote:
> Would it be possible to get a multivariate gamma density in STAN?
>
> Coding-wise, it would not be difficult at all to copy what we have for the Wishart distribution and modifying it slightly for the multivariate gamma. The problem is finding someone who is sufficiently motivated to do it because in order to put it into Stan it would also have to be tested at least as much as the Wishart distribution, documented, complimented with a function to draw from the multivariate gamma, etc.

I think of the testing and API doc as part of the total "coding" process.

> And that is all pretty tedious to do for a distribution that is not widely used and is not likely to lead to efficient sampling in Stan unless X were data.
>
> If you are that motivated, then we can help you jump through the hoops involved when adding new functions to Stan. But if you are not so motivated, it is pretty easy to write that log density as a user-defined Stan function within a Stan program. And depending on what special case (diagonal scale matrix, etc.) you might be interested in, it may actually run faster than would a general implementation in the Stan Math Library.

We obviously need more specialized matrix operations in
the math library. We've been kicking around the idea of
having various specializations, but the syntax is going to
be ugly. Things like multiply_lower_triangular(matrix, real)
and multiple_lower_triangular(matrix, vector) and so on.
Plus all the positive-definite versions, and the symmetric
versions, multiply_symmetric(...), ... I don't see how we
can overload the basic operations the way things are written now,
because run-time typing is just int / real / matrix / vector / row vector
and we don't have basic types for type inference for the triangular
and symmetric and other types of matrices we might want.

We could also add new data types the same way we added Cholesky.
That's not too hard to do --- Marcus managed to add one without
much help, but then he's a very skilled coder, so others' mileage
may vary.

- Bob

[Bob Carpenter]

There are two QR-decomposition functions (qr_Q, qr_R) for the two matrices.

It would be more efficient if we had a way to return both qr_Q and qr_R as
a single object --- we should probably have a function that just concatenates
one on top of the other, which would be easy to decompose back in the Stan code.
But we don't have that yet.

- Bob


> On Nov 5, 2015, at 6:36 AM, Royce Anders <ander...@gmail.com> wrote:

>
>
>
> On Wed, Nov 4, 2015 at 7:03 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> > On Nov 4, 2015, at 11:11 AM, Ben Goodrich <goodri...@gmail.com> wrote:
> >
> > On Wednesday, November 4, 2015 at 10:02:38 AM UTC-5, Royce Anders wrote:
> > Would it be possible to get a multivariate gamma density in STAN?
> >
> > Coding-wise, it would not be difficult at all to copy what we have for the Wishart distribution and modifying it slightly for the multivariate gamma. The problem is finding someone who is sufficiently motivated to do it because in order to put it into Stan it would also have to be tested at least as much as the Wishart distribution, documented, complimented with a function to draw from the multivariate gamma, etc.
>
> I think of the testing and API doc as part of the total "coding" process.
>
> > And that is all pretty tedious to do for a distribution that is not widely used and is not likely to lead to efficient sampling in Stan unless X were data.
> >
> > If you are that motivated, then we can help you jump through the hoops involved when adding new functions to Stan. But if you are not so motivated, it is pretty easy to write that log density as a user-defined Stan function within a Stan program. And depending on what special case (diagonal scale matrix, etc.) you might be interested in, it may actually run faster than would a general implementation in the Stan Math Library.
>

> I think indeed an interesting first step could be to try the log density as a user-defined Stan function. I'd be happy to report on it, and if the experiences pove worthwhile to implement it in STAN formally, I'd be happy to also help with its implementation.
>
> I took a try at a log density user-defined Stan function for the multivariate gamma, and according to Bob's comment below, I believe I got stuck due to limited matrix operations.
>
> I followed a technique based on Song (2000), in which one uses normal copula's, involving QR decomposition and solving the linear inverse problem of the QR decomposition, to estimate the density. I believe with this third step, is where I couldn't find such a function for STAN. According to your note here, perhaps I'm highly overcomplicating the calculation of this density, and I would be happy to receive help. I'm looking for the aka matrix gamma distribution, parameterized by (shape, scale, E), where E is the symmetric correlation matrix.
>
> Otherwise, if there are simple STAN codes to easily get these R equivalents, I might be able to figure this out already:
> If E is a pxp symmetric correlation matrix, I would need the following STAN equivalents for this code in R:
> 1) The diagonal of qr(E)$qr ## here $qr is a pxp matrix, in which the upper triangle is R and the lower is Q; and
> 2) qr.solve(qr(E)) ## which solves the linear inverse problem, giving a pxp matrix.
>
> again maybe there is a much easier solution to calculating this density than what I've been looking at.

>
>
> We obviously need more specialized matrix operations in
> the math library. We've been kicking around the idea of
> having various specializations, but the syntax is going to
> be ugly. Things like multiply_lower_triangular(matrix, real)
> and multiple_lower_triangular(matrix, vector) and so on.
> Plus all the positive-definite versions, and the symmetric
> versions, multiply_symmetric(...), ... I don't see how we
> can overload the basic operations the way things are written now,
> because run-time typing is just int / real / matrix / vector / row vector
> and we don't have basic types for type inference for the triangular
> and symmetric and other types of matrices we might want.
>
> We could also add new data types the same way we added Cholesky.
> That's not too hard to do --- Marcus managed to add one without
> much help, but then he's a very skilled coder, so others' mileage
> may vary.
>
> - Bob
>

> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/BBn8U67hcbQ/unsubscribe.
> To unsubscribe from this group and all its topics, 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.
>
>
> --
> 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.

[Ben Goodrich]

According to your note here, perhaps I'm highly overcomplicating the calculation of this density, and I would be happy to receive help. I'm looking for the aka matrix gamma distribution, parameterized by  (shape, scale, E), where E is the symmetric correlation matrix.

This one, right?

https://en.wikipedia.org/wiki/Matrix_gamma_distribution

Anything involving a QR decomposition is not going to be very efficient in Stan. The best way to do something like this would be to parameterize it in terms of the Cholesky factor, L, of Sigma = L * L'. So you get a log-PDF

(-2 * alpha) * log(det(L)) - (p * alpha) * log(beta)
- 0.25 * p * (p - 1) * sum(lgamma(alpha), ..., lgamma(alpha + 0.5 * (1 - p)))
+ (alpha - 0.5 * (p - 1)) * log(det(X)) - 1 / beta * trace( L^{-1}' * L^{-1} * X )

• log(det(L)) is straightforward because it is the sum of the logarithms of the diagonal elements of L
• if X is fixed then you can omit the term involving log(det(X)) or calculate it once in transformed data; if X is a parameter utilize Stan's log_determinant(X) function, which is more efficient than log(determinant(X))
• Rewrite trace( L^{-1}' * L^{-1} * X ) as trace( L^{-1} * X * L^{-1}' )  by the cyclic permutation property and utilize Stan's trace_quad_form(L_inv_t, X) function for efficiency. This will require you compute L_inv_t = L^{-1}', which can be done efficiently utilizing Stan's transpose(mdivide_left_tri_low(L, I))
  
• If you think of Sigma as a correlation matrix, then L needs to a parameter of type cholesky_factor_corr[p] and a good prior on L is lkj_corr_cholesky(1.0). If you think of Sigma as a covariance matrix, you will need to fix its trace, in which case a good choice for the vector of variances is a simplex. You can utilize Stan's diag_pre_multipy() function to multiply a diagonal matrix of standard deviations by the Cholesky factor of a correlation matrix to get a Cholesky factor of a covariance matrix.

By all means, use the expose_stan_functions() function in rstan to verify that your user-defined log-PDF returns the same number as whatever R implementation you have at your disposal.

Ben

[Royce Anders]

Thank you for the very thorough and helpful responses all!! I will have a go at implementing the advice and return with news.  

One quick thing to clarify if you please: in the general case, if one uses a user-defined distribution as a parameter prior; then instead of using, for example:

for(i in 1:nparams){

y[i] ~ dcustom_prior(a,b);

}

One would instead use 

for(i in 1:nparams){

increment_log_prob(dcustom_prior(y[i],a,b));

}

?

[Bob Carpenter]

See the chapter in the manual on user-defined functions.
The short answer is that if the function is foo_log(y, ...)
then you can call it as y ~ foo(...) in a Stan program.

- 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.