Mixed Discrete-Continuous Gaussian Copula Functions

Functions
array[] matrix	normal_marginal (matrix y, matrix mu_glm, vector sigma)
array[] matrix	bernoulli_marginal (array[,] int y, matrix mu_glm, matrix u_raw)
array[] matrix	binomial_marginal (array[,] int num, array[,] int den, matrix mu_glm, matrix u_raw)
array[] matrix	poisson_marginal (array[,] int y, matrix mu_glm, matrix u_raw)
real	centered_gaussian_copula_cholesky_lpdf (array[,] matrix marginals, matrix L)

Detailed Description

The mixed-discrete normal copula. The method is from Smith, Michael & Khaled, Mohamad. (2011). Estimation of Copula Models With Discrete Margins via Bayesian Data Augmentation. Journal of the American Statistical Association. 107. 10.2139/ssrn.1937983.

 
#include copula/normal_copula.stanfunctions
 
array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) {
  int N = rows(mu_glm);
  int J = cols(mu_glm);
  array[2] matrix[N, J] rtn;
  // Initialise the jacobian adjustments to zero, as vectorised lpdf will be used
  rtn[2] = rep_matrix(0, N, J);
  
  for (j in 1 : J) {
    rtn[1][ : , j] = (y[ : , j] - mu_glm[ : , j]) / sigma[j];
    rtn[2][1, j] = normal_lpdf(y[ : , j] | mu_glm[ : , j], sigma[j]);
  }
  
  return rtn;
}
 
array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
  int N = rows(mu_glm);
  int J = cols(mu_glm);
  matrix[N, J] matrix_y = to_matrix(y);
  matrix[N, J] mu_glm_logit = 1 - inv_logit(mu_glm);
  
  matrix[N, J] Lbound = matrix_y .* mu_glm_logit;
  matrix[N, J] UmL = fabs(matrix_y - mu_glm_logit);
  
  return {inv_Phi(Lbound + UmL .* u_raw), log(UmL)};
}
 
array[] matrix binomial_marginal(array[,] int num, array[,] int den, matrix mu_glm,
                                 matrix u_raw) {
  int N = rows(mu_glm);
  int J = cols(mu_glm);
  matrix[N, J] mu_glm_logit = inv_logit(mu_glm);
  array[2] matrix[N, J] rtn;
  
  for (j in 1 : J) {
    for (n in 1 : N) {
      real Ubound = binomial_cdf(num[n, j] | den[n, j], mu_glm_logit[n, j]);
      real Lbound = 0;
      if (num[n, j] > 0) {
        Lbound = binomial_cdf(num[n, j] - 1 | den[n, j], mu_glm_logit[n, j]);
      }
      real UmL = Ubound - Lbound;
      rtn[1][n, j] = inv_Phi(Lbound + UmL * u_raw[n, j]);
      rtn[2][n, j] = log(UmL);
    }
  }
  
  return rtn;
}
 
array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
  int N = rows(mu_glm);
  int J = cols(mu_glm);
  matrix[N, J] mu_glm_exp = exp(mu_glm);
  array[2] matrix[N, J] rtn;
  
  for (j in 1 : J) {
    for (n in 1 : N) {
      real Ubound = poisson_cdf(y[n, j] | mu_glm_exp[n, j]);
      real Lbound = 0;
      if (y[n, j] > 0) {
        Lbound = poisson_cdf(y[n, j] - 1 | mu_glm_exp[n, j]);
      }
      real UmL = Ubound - Lbound;
      rtn[1][n, j] = inv_Phi(Lbound + UmL * u_raw[n, j]);
      rtn[2][n, j] = log(UmL);
    }
  }
  
  return rtn;
}
 
real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) {
  // Extract dimensions of final outcome matrix
  int N = rows(marginals[1][1]);
  int J = rows(L);
  matrix[N, J] U;
  
  // Iterate through marginal arrays, concatenating the outcome matrices by column
  // and aggregating the log-likelihoods (from continuous marginals) and jacobian
  // adjustments (from discrete marginals)
  real adj = 0;
  int pos = 1;
  for (m in 1 : size(marginals)) {
    int curr_cols = cols(marginals[m][1]);
    
    U[ : , pos : (pos + curr_cols - 1)] = marginals[m][1];
    
    adj += sum(marginals[m][2]);
    pos += curr_cols;
  }
  
  // Return the sum of the log-probability for copula outcomes and jacobian adjustments
  return multi_normal_cholesky_copula_lpdf(U | L) + adj;
}
 

Function Documentation

bernoulli_marginal()

array[] matrix bernoulli_marginal	(	array int	y[,],
		matrix	mu_glm,
		matrix	u_raw
	)

Bernoulli marginal

Bernoulli marginal for mixed continuous-discrete Gaussian copula.

The lb and ub: if y == 0, upper bound at inv_Phi( y, mu ) if y == 1, lower bound at inv_Phi( y - 1, mu )

Copyright

Ethan Alt, Sean Pinkney 2021

Parameters

y	int[,] 2d array of binary outcomes
mu_glm	matrix of regression means
u_raw	matrix of nuisance latent variables

Returns

2D array of matrices containing the random variables and jacobian adjustments

binomial_marginal()

array[] matrix binomial_marginal	(	array int	num[,],
		array int	den[,],
		matrix	mu_glm,
		matrix	u_raw
	)

Binomial marginal

Binomial marginal for mixed continuous-discrete Gaussian copula.

The lb and ub: Always upper bound at inv_Phi( y, mu ) If n != 0, lower bound at inv_Phi(n-1, mu)

Copyright

Ethan Alt, Sean Pinkney 2021

Parameters

num	int[,] 2D array of numerator integers
den	int[,] 2D array of denominator integers
mu_glm	matrix of regression means
u_raw	matrix of nuisance latent variables

Returns

2D array of matrices containing the random variables and jacobian adjustments

centered_gaussian_copula_cholesky_lpdf()

real centered_gaussian_copula_cholesky_lpdf	(	array matrix	marginals[,],
		matrix	L
	)

Mixed Copula Log-Probability Function

Copyright

Andrew Johnson 2022

Parameters

marginals	Nested arrays of matrices from marginal calculations
L	Cholesky Factor Correlation

Returns

Real log-probability

normal_marginal()

array[] matrix normal_marginal	(	matrix	y,
		matrix	mu_glm,
		vector	sigma
	)

Normal marginal

Standardized normal marginal for mixed continuous-discrete Gaussian copula.

Copyright

Ethan Alt, Sean Pinkney 2021

Parameters

y	matrix of normal outcomes
mu_glm	row_vector of regression means
matrix	vector of outcome SD's

Returns

2D array of matrices containing the random variables and jacobian adjustments

poisson_marginal()

array[] matrix poisson_marginal	(	array int	y[,],
		matrix	mu_glm,
		matrix	u_raw
	)

Poisson marginal

Poisson marginal for mixed continuous-discrete Gaussian copula.

The lower-bound and upper-bound:
The upper-bound is always at

inv_Phi( y, mu ) 

If \(y \ne 0\), lower-bound at

inv_Phi(y-1, mu) 

.
At \(y = 0\) the lower-bound is \( 0 \).

Copyright

Ethan Alt, Sean Pinkney 2021

Parameters

y	int[,] 2D array of integer counts
mu_glm	matrix of regression means
u_raw	matrix of nuisance latent variables

Returns

2D array of matrices containing the random variables and jacobian adjustments