Help Wanted: Custom Multivariate pdf For State Space Model

General
, , ,
1

Hello,

I’m trying to fit a state space model of the form;

x_t = Gx_{t-1} + \epsilon_s
y_t = Fx_t + \epsilon_o

I would like to specify a user-defined PDF (multivariate) on the errors, \epsilon_s, \epsilon_o. Currently I’ve been using the below Stan code to sample the state sequence, x_t, and state transition matrix, G ;

data{  
int<lower=1> t_max;                // Time series length
  array[t_max]  vector[5] y;       // Observations
  
  cov_matrix[5]      W;       // Variance of state noise
  cov_matrix[5]      V;       // Variance of observation noise
  row_vector[5]     m0;       // Mean of prior distribution--A 2D row vector
  cov_matrix[5]     C0;       // Variance of prior distribution

  matrix[5,5]     Fmat;       // Factor loadings
}

parameters{
  
  row_vector[5]    x0;            // State [0]
  array[t_max]  vector[5]   x;    // State [1:t_max]
  
  vector<lower=-1, upper=1>[5]     Gmat;       // Vector containing autocorrelations

}

 model{
 
 // Pre-compute array of mean vectors;
array[t_max]  vector[5] y_mu;
array[t_max]  vector[5] x_mu_prev;

for(t in 1:t_max){

    if(t == 1){
    
    x_mu_prev[t] = diag_matrix(Gmat) * x0';
    
    }else{
    
    x_mu_prev[t] = diag_matrix(Gmat) * x[t-1];
    
    }
    
    y_mu[t] = Fmat * x[t];
      
 }

  // Initial state
x0 ~ multi_normal(m0, sqrt(C0));
  
  target += multi_normal_lpdf(x | x_mu_prev, sqrt(W)) + multi_normal_lpdf( y | y_mu, sqrt(V));
  
  //Prior on the diagonal elements of Gmat;
  for(n in 1:5){
  
  Gmat[n] ~ uniform(-1, 1);
  
  }

Which performs satisfactorily. But when I define the PDF;

functions{

real MVPE_lpdf( vector  x, vector mu, matrix sigma, real beta){

int n = num_elements(x);

return log(n) + lgamma(n * 0.5) -0.5*log_determinant(sigma) - 0.5*n*log(pi()) -lgamma(1 + n/(2*beta)) -(1+n/(2*beta))*log(2) -0.5*((x - mu)' * inverse(sigma) * (x - mu))^beta;

  }

}

And run the model;

 model{
 
 // Pre-compute array of mean vectors;
array[t_max]  vector[5] y_mu;
array[t_max]  vector[5] x_mu_prev;

for(t in 1:t_max){

    if(t == 1){
    
    x_mu_prev[t] = diag_matrix(Gmat) * x0;
    
    }else{
    
    x_mu_prev[t] = diag_matrix(Gmat) * x[t-1];
    
    }
    
    y_mu[t] = Fmat * x[t];
      
 }

  
  // Initial state
x0 ~ multi_normal(m0, sqrt(C0));
  
  target += MVPE_lpdf(x | x_mu_prev, W, beta_s) + MVPE_lpdf(y | y_mu, V, beta_o);

  //Prior on the diagonal elements of Gmat;  

  for(n in 1:5){
  
  Gmat[n] ~ uniform(-1, 1);
  
  }

I get the following error message;

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  0

Semantic error in 'string', line 54, column 12 to column 47:

Ill-typed arguments supplied to function 'MVPE_lpdf'. Available signatures:
(vector, vector, matrix, real) => real
Instead supplied arguments of incompatible type: vector[], vector[], matrix, real.

So, I would like to have my MVPE_lpdf function behave like multi_normal_lpdf does in the above model block, but I don’t know how to accomplish this. My guess is that I want MVPE_lpdf to accept an array of vectors(?) for its first two arguments. I’ve dabbled with changing the argument types, but to no avail.

Thank you!

2

I think that the last line of the error message tells you that you are passing in an array of vectors, so you need extra brackets behind vector when declaring the function:

real MVPE_lpdf( vector[]  x, vector[] mu, matrix sigma, real beta){

but then I don’t think that the last term of the loglikelihood with (x - mu) will work correctly because you are sending in many vectors at once. For that last term, you might invert sigma, then loop through n times and perhaps also use quad_form_sym().

3

Hi,

After following your advice (and correcting some of my own errors), I arrived at the following function;

real MVPE_lpdf( vector[]  x, vector[] mu, matrix sigma, real beta){
  
  // Dimension of the state process
  int n = num_elements(x[1]);
  
  // Length of the time series
  int T = dims(x)[1];
  
  // Dimensions of sigma
  int p = cols(sigma);
  
 matrix[p,p] sigma_inverse = inverse(sigma);
 
 real quad_term = 0.0;
  
  for(t in 1:T){
    
    vector[n] B = x[t] - mu[t];
    
    quad_term += quad_form_sym(sigma_inverse, B)^beta;
    
  }
  
  return T*(log(n) + lgamma(n * 0.5) -0.5*log_determinant(sigma) - 0.5*n*log(pi()) -lgamma(1 + n/(2*beta)) -(1+n/(2*beta))*log(2)) -0.5*quad_term;

}

Thank you for your help–You were spot on about (x-mu) breaking. And thank you for suggesting the quad_form_sym() function.

To anyone who might find this post: This density is the Multivariate Power Exponential, and for beta = 1, you obtain the multi_normal_lpdf. After following @edm’s advice and correcting mistakes, I was able to obtain the exact same output from the above code block, when using multi_normal_lpdf and MVPE_lpdf with beta = 1.

Thanks again for your help, I really appreciate it!

1 Like