Normal copula with Rstan

Modeling
1

Implementing a Gaussian Copula with Gamma Margins in Stan - Numerical Stability Issues

Model Description

I’m implementing a model in Stan where the latent parameters lambda, mu, and nu have Gamma marginal priors and are correlated via a Gaussian copula.

Here’s the implementation of the multi_normal_cholesky_copula function and the model block:

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

model {
  matrix[N,3] u; // To hold the three uniforms
  
  // 1) Marginal Gamma priors
  lambda ~ gamma(r, a);
  mu ~ gamma(s, b);
  nu ~ gamma(q, c);
  
  // 2) Build uniforms via their marginal CDFs
  for (i in 1:N) {
    u[i, 1] = gamma_cdf(lambda[i] | r, a);
    u[i, 2] = gamma_cdf(mu[i] | s, b);
    u[i, 3] = gamma_cdf(nu[i] | q, c);
  }
  
  // 3) Copula ties the three margins together
  target += multi_normal_cholesky_copula_lpdf(u | L);

Is this the correct and stable way to specify Gamma marginals with a Gaussian copula in Stan?
Also I got these warning errors:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: gamma_lpdf: Random variable[65] is inf, but must be positive finite! (in 'C:/Users/riana/AppData/Local/Temp/RtmpA3j600/model-73807ba7106.stan', line 58, column 2 to column 23)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Thank you for your help.
2

Hi, the quoted function doesn’t take parameter u (a unit vector). Instead, you should feed it a standard (multivariate) normally distributed vector, which is constructed from what you now call parameter u. What is missing from your code is a translation of u from the unit scale to the unbounded scale, using the normal copula. This can be achieved simply by the inv_Phi() function (approximate but fast):

should be

Instead of inv_Phi(), you could also use the inverse CDF of the normal distribution, which is available in Stan. But in my experience, it is slower (more computation involved) and can suffer from numerical instability. My advice would be to stick with the polynomial approximation via inv_Phi().

EDIT: I did not check your definition of the multi_normal_cholesky_copula_lpdf() function, but I’m assuming you adopted @spinkney’s implementation as described on their github here.

1 Like
3

@bgautijonsson did a presentation on copulas with tutorial material at StanCon last year. Here’s his presentation on YouTube!

and he’s working on adding a chapter to the Stan User’s Guide on copulas.

2 Likes
4

Hey, Valera.

I have a small example with Exponential marginals in the following blog post: A Gentle Introduction: The Gaussian Copula – bggj

It would be great to see the rest of the model if you’d like more comments.

1 Like
5

Thank you for your response. yes I adopted @spinkney’s implementation.

6

Thank you, I will watch it.

7

Thank you very much. After adjustment, here is my complete stan code.

Define the Stan model code as a string

stan_code ← "

functions {
  // Define the normal copula log-likelihood for multivariate data
  real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) {
    int N = rows(U);
    int J = cols(U);
    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=1> N;           // Number of samples
  vector[N] x;              // Data vector for x
  vector[N] k;              // Data vector for t_x
  vector[N] T;              // Data vector for T
  vector[N] m_x;            // Data vector for m_x
  
  // Hyperparameters for priors
  real<lower=0> r;          // Shape parameter for lambda's prior
  real<lower=0> a;          // Rate parameter for lambda's prior
  real<lower=0> s;          // Shape parameter for mu's prior
  real<lower=0> b;          // Rate parameter for mu's prior
  real<lower=0> p_;         // Parameter for nu's likelihood contribution
  real<lower=0> q;          // Shape parameter for nu's prior
  real<lower=0> c;          // Rate parameter for nu's prior
  
  // Correlation matrix entries for copula
  real<lower=-1,upper=1> rho_lm;
  real<lower=-1,upper=1> rho_ln;
  real<lower=-1,upper=1> rho_mn;
}

transformed data {
  matrix[3,3] Omega;
  Omega = [ 
    [1,      rho_lm, rho_ln],
    [rho_lm, 1,      rho_mn],
    [rho_ln, rho_mn, 1     ]
  ];
  
  // Cholesky decomposition of the correlation matrix
  matrix[3,3] L = cholesky_decompose(Omega);
}

parameters {
  vector<lower=1e-8>[N] lambda;
  vector<lower=1e-8>[N] mu;
  vector<lower=1e-8>[N] nu;
}

model {
  matrix[N, 3] z;
  real eps = 1e-10;

  // Marginal priors
  lambda ~ gamma(r, a);
  mu ~ gamma(s, b);
  nu ~ gamma(q, c);

  // Transform to Gaussian copula space
  for (i in 1:N) {
    z[i, 1] = inv_Phi(fmin(fmax(gamma_cdf(lambda[i] | r, a), eps), 1 - eps));
    z[i, 2] = inv_Phi(fmin(fmax(gamma_cdf(mu[i] | s, b), eps), 1 - eps));
    z[i, 3] = inv_Phi(fmin(fmax(gamma_cdf(nu[i] | q, c), eps), 1 - eps));

  }

  // Apply copula likelihood
  target += multi_normal_cholesky_copula_lpdf(z | L);

  // 4) Likelihood model: per-chain likelihood
  for (i in 1:N) {
    real x_i = x[i];
    real k_i = k[i];
    real T_i = T[i];
    real m_i = m_x[i];

    real log_lik1 = log(fmax(eps, x_i * log(lambda[i]) + log(mu[i]) 
                  - log(lambda[i] + mu[i]) 
                  - k_i * (lambda[i] + mu[i])));
    real log_lik2 = log(fmax(eps, (x_i + 1) * log(lambda[i]) 
                  - log(lambda[i] + mu[i]) 
                  - T_i * (lambda[i] + mu[i])));
    real log_Aik = log_sum_exp(log_lik1, log_lik2);


    if (x_i == 0) {
      target += log_Aik;
    } else {
      real log_ABik = log_Aik
                    + (x_i * p_ - 1) * log(m_i) 
                    - m_i * x_i * nu[i]
                    - lgamma(x_i * p_)
                    + p_ * x_i * log(x_i * nu[i]);
      target += log_ABik;
    }
  }
}

[edit: escaped Stan code]

8

As the warning error says, if it’s sporadic (typically at the start of sampling) then it should be fine and you won’t have to worry about it.

Regarding your code:

  1. You could create L as a cholesky_factor_corr instead of manually creating all the parts of the correlation matrix.
  2. Is the fmin/fmax code because of the warnings?
1 Like
9

Yes, the fmin/fmax is set because of the warnings. After making adjustments - such as using cholesky_factor_corr and applying inv_Phi() - the warning still appears at the beginning, but as you mentioned, MCMC sampling runs successfully. Many thanks for your help!

1 Like