Gaussian copula for spherical data

Modeling
6

True, it’s not complex, it’s the unit circle. I assumed they use complex valued pdfs.

This may be interesting too:

But can you try this function:

real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) {
    int N = rows(U); // total number of observations
    int J = cols(U); // total number of responses (J+K)
    matrix[J, J] Gammainv = chol2inv(L);
    return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U));
  }

functions {  
real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) {
    int N = rows(U); // total number of observations
    int J = cols(U); // total number of responses (J+K)
    matrix[J, J] Gammainv = chol2inv(L);
    return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U));
  }
}

data {
  int<lower=0> N;
  int<lower=0> K;
  matrix[N, K] Y;
}

parameters {
  vector<lower=0>[K] kappa;
  vector<lower=-pi(), upper=pi()>[K] mu;
  cholesky_factor_corr[K] rho_chol;
}

model {

  rho_chol ~ lkj_corr_cholesky(2.0);
  kappa ~ exponential(0.05);
  
  mu ~ von_mises(0.0, 0.0);
  
  matrix[N, K] u;

  for (n in 1:N) {
    for (k in 1:K){
       u[n, k] = std_normal_log_qf(von_mises_lcdf(Y[n, k] | mu[k], kappa[k]));
    }
  }
  
  for (k in 1:K){
    Y[,k] ~ von_mises(mu[k], kappa[k]);
  }

  u ~ multi_normal_cholesky_copula(rho_chol);
}

generated quantities {
  matrix[K,K] rho = multiply_lower_tri_self_transpose(rho_chol);
}