Error in the custom multivariate normal density function

General
1

Hi, I’m working on the custom multivariate normal density function, which is:

functions {
  real mvn_lpdf(vector X, vector mu, cov_matrix sigma) {
  
    real e_part = 0.0;
    real logl = 0.0;
 
    e_part = (X - mu)' * inverse(sigma) * (X - mu);
  
    logl = log(2*pi()) + log(determinant(sigma))
    logl = -.5*logl - .5*e_part

    return logl
  }
}

When I run this code with the proper data, it gave me an error looking like a simple one, but I don’t have a clue where it comes.

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
 error in 'model41f44c6a75b0_example_model' at line 2, column 35
  -------------------------------------------------
     1: functions {
     2:   real mvn_lpdf(vector X, vector mu, cov_matrix sigma) {
                                         ^
     3:   
  -------------------------------------------------

PARSER EXPECTED: <argument declaration or close paren ) to end argument declarations>
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'example_model' due to the above error.

Any solution to address this error?

Best,

1 Like
2

For user defined functions you can’t input a constrained type. Declare sigma as a matrix instead of a cov_matrix.

Is there a reason why you’re writing your own? Your implementation above will be much slower than the inbuilt. Or even writing it in a more computationally efficient way. A few notes,

  • we try to avoid explicitly calling inverse on a matrix
  • the cholesky parameterization will be faster and more numerically stable
  • there is a log_determinant function

If you have a matrix X where the rows equal the number of observations then a more efficient version of the mvn using the cholesky parameterization in the functions block will look like

 real mvn_cholesky_lpdf(matrix X, matrix mu, matrix L) {
    int K = rows(L);
    int N = rows(X);
    real sqrt_det = -N * sum(log(diagonal(L)));
    real norm_const =  -K * N * 0.5 * log(2 * pi());
    real mahab = sum(columns_dot_self(mdivide_left_tri_low(L, (X' - mu))));
    
    return norm_const + sqrt_det - 0.5 * mahab;
  }
1 Like
3

Thank you for your answer.

I’m trying to do some simulation about Bayesian analysis, using Texas Advanced Computing Center (TACC). For some reason, multi_norm() function didn’t work on TACC. I guess something related to compliers. So, I’m testing my own function based on the very basic stan functions to avoid such errors.

Thank you again. I’ll try your code on TACC.

1 Like