Skip to content

Commit

Permalink
Formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
actions-user committed Mar 12, 2022
1 parent 2c06a47 commit e9cac29
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 36 deletions.
51 changes: 25 additions & 26 deletions functions/copula/centered_gaussian_copula.stanfunctions
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,12 @@ array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) {
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]);
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;
}

Expand All @@ -63,11 +63,11 @@ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
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) };
return {inv_Phi(Lbound + UmL .* u_raw), log(UmL)};
}

/**
Expand All @@ -89,15 +89,15 @@ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
* @return 2D array of matrices containing the random variables
* and jacobian adjustments
*/
array[] matrix binomial_marginal(array[,] int num, array[,] int den,
matrix mu_glm, matrix u_raw) {
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) {
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) {
Expand All @@ -108,7 +108,7 @@ array[] matrix binomial_marginal(array[,] int num, array[,] int den,
rtn[2][n, j] = log(UmL);
}
}

return rtn;
}

Expand All @@ -135,9 +135,9 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
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) {
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) {
Expand All @@ -148,11 +148,10 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
rtn[2][n, j] = log(UmL);
}
}

return rtn;
}


/**
* Mixed Copula Log-Probability Function
*
Expand All @@ -167,23 +166,23 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L)
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)) {
for (m in 1 : size(marginals)) {
int curr_cols = cols(marginals[m][1]);

U[ : , pos:(pos + curr_cols - 1)] = 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;
}

/** @} */
/** @} */
19 changes: 9 additions & 10 deletions functions/copula/normal_copula.stanfunctions
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
*/
real normal_copula_lpdf(real u, real v, real rho) {
real rho_sq = square(rho);

return (0.5 * rho * (-2. * u * v + square(u) * rho + square(v) * rho)) / (-1. + rho_sq)
- 0.5 * log1m(rho_sq);
}
Expand All @@ -52,12 +52,12 @@ real normal_copula_lpdf(real u, real v, real rho) {
real normal_copula_vector_lpdf(vector u, vector v, real rho) {
int N = num_elements(u);
real rho_sq = square(rho);

real a1 = 0.5 * rho;
real a2 = rho_sq - 1;
real a3 = 0.5 * log1m(rho_sq);
real x = -2 * u' * v + rho * (dot_self(u) + dot_self(v));

return a1 * x / a2 - N * a3;
}

Expand Down Expand Up @@ -97,7 +97,7 @@ real multi_normal_copula_lpdf(matrix u, matrix L) {
int N = cols(u);
real inv_sqrt_det_log = sum(log(diagonal(L)));
matrix[K, N] x = mdivide_left_tri_low(L, u);

return -N * inv_sqrt_det_log - 0.5 * sum(columns_dot_self(x) - columns_dot_self(u));
}

Expand All @@ -122,10 +122,10 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
real alpha_u = a * (pv / pu - rho);
real alpha_v = a * (pu / pv - rho);
real d = 0;

if (u[1] < 0.5 && u[2] >= 0.5 || u[1] >= 0.5 && u[2] < 0.5)
if (u[1] < 0.5 && u[2] >= 0.5 || u[1] >= 0.5 && u[2] < 0.5)
d = 0.5;

return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d;
}

Expand Down Expand Up @@ -165,7 +165,6 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
* @return log density of the distribution
*/


/**
* Multivariate Normal Copula lpdf (Cholesky parameterisation)
*
Expand All @@ -176,7 +175,7 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) {
int N = rows(U);
int J = cols(U);
matrix[J,J] Gammainv = chol2inv(L);
matrix[J, J] Gammainv = chol2inv(L);
return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U));
}
/** @} */
/** @} */

0 comments on commit e9cac29

Please sign in to comment.