Helpful Stan Functions
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;
}
array[] matrix binomial_marginal(array[,] int num, array[,] int den, matrix mu_glm, matrix u_raw)
Definition: centered_gaussian_copula.stanfunctions:92
real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L)
Definition: centered_gaussian_copula.stanfunctions:164
array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma)
Definition: centered_gaussian_copula.stanfunctions:28
array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw)
Definition: centered_gaussian_copula.stanfunctions:61
array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw)
Definition: centered_gaussian_copula.stanfunctions:133
real multi_normal_cholesky_copula_lpdf(matrix U, matrix L)
Definition: normal_copula.stanfunctions:130
*copyright Sean *author Sean Pinkney *param x Vector *param mu Real *param sigma Real f$ in(0, \infty)\f$ scale parameter * @param lambda Real \f$-1< \lambda< 1\f$ * @param q Real \f$\in(0

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 )

Parameters
yint[,] 2d array of binary outcomes
mu_glmmatrix of regression means
u_rawmatrix 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)

Parameters
numint[,] 2D array of numerator integers
denint[,] 2D array of denominator integers
mu_glmmatrix of regression means
u_rawmatrix 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

Parameters
marginalsNested arrays of matrices from marginal calculations
LCholesky 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.

Parameters
ymatrix of normal outcomes
mu_glmrow_vector of regression means
matrixvector 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 \).

Parameters
yint[,] 2D array of integer counts
mu_glmmatrix of regression means
u_rawmatrix of nuisance latent variables
Returns
2D array of matrices containing the random variables and jacobian adjustments