efficiency of stan code

[spinkney]

This is about twice as fast for me. Of course when you first wrote the code Stan didn't have tuples nor did it have the cholesky version of the inverse Wishart.

functions {
  matrix square_root (matrix A) {
    // assume A is SPD
    int R = rows(A);
    tuple(matrix[R, R], vector[R]) eig = eigendecompose_sym(A);
    return diag_post_multiply(eig.1, sqrt(eig.2)) * eig.1';
  }
  
  real log_norm_const(int d, real rho) {
    real num = (d - 1) * log2() + log1p(inv_sqrt(1 - square(rho)));
    real den = d * log(sqrt(1 - rho) + sqrt(1 + rho));
    return num - den;
  }
  
  real multi_expectile_lpdf (array[] vector x, vector m, matrix L, vector nu, real r) {
    int d = rows(L);
    int n = dims(x)[1];
    // real lpdf = multi_normal_cholesky_lpdf(x | m, L);
   // lpdf -= n * log_norm_const(d, r);
    
    matrix[d, d] sigma_inv = chol2inv(L);
    matrix[d, d] sigma_inv_sqrt = square_root(sigma_inv);
    real lpdf = multi_normal_prec_lpdf(x | m, sigma_inv);
    lpdf -= n * log_norm_const(d, r);
    
    real cache = 0;
    vector[d] sigma_inv_sqrt_nu = sigma_inv_sqrt * nu;
    for (i in 1:n) {
      vector[d] x_m = x[i] - m; 
      cache += dot_product(x_m, sigma_inv_sqrt_nu) * sqrt(quad_form(sigma_inv, x_m));
    }
    
    return lpdf - 0.5 * r * cache;
    
  }
}
data {
  int<lower=1> d;         // data dimension
  int<lower=0> n;         // number of observations
  array[n] vector[d] Y;         // observations
}
parameters {
  vector[d] mu;                     // mean
  real<lower=0,upper=1> rho;        // asymmetry parameter
  cholesky_factor_cov[d] L;          
  vector<lower=0>[d] diag_sigma_sq; // mu hyperprior cov mat
  vector[d] nu_tilde_raw;          // direction (tilde version)
}
transformed parameters {
  real r = sqrt(dot_self(nu_tilde_raw));
  vector[d] nu_tilde = nu_tilde_raw / r;
}
model {
  target += -0.5 * r^2 + 10 * log(r);
  mu ~ normal(0, diag_sigma_sq);
  L ~ inv_wishart_cholesky(d, diag_matrix(rep_vector(1, d)));
  diag_sigma_sq ~ inv_gamma(1, 1);
  Y ~ multi_expectile(mu, L, nu_tilde, rho);
}